diff --git a/analysis/src/main/java/org/hps/analysis/MC/SvtEventForensicsDriver.java b/analysis/src/main/java/org/hps/analysis/MC/SvtEventForensicsDriver.java new file mode 100644 index 000000000..aa03472cc --- /dev/null +++ b/analysis/src/main/java/org/hps/analysis/MC/SvtEventForensicsDriver.java @@ -0,0 +1,397 @@ +package org.hps.analysis.MC; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; + +import hep.physics.vec.Hep3Vector; + +import org.lcsim.event.EventHeader; +import org.lcsim.event.LCRelation; +import org.lcsim.event.MCParticle; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.ReconstructedParticle; +import org.lcsim.event.SimTrackerHit; +import org.lcsim.event.Track; +import org.lcsim.event.TrackState; +import org.lcsim.event.TrackerHit; +import org.lcsim.event.Vertex; +import org.lcsim.util.Driver; + +/** + * Full per-hit dump of a handful of named events, for debugging a specific case where a + * real track lost its truth link under pulser overlay. + * + * SvtHitProvenanceDriver answers the aggregate question. This one answers the forensic + * question for one event, and in particular separates the two ways a real particle's + * track can end up with zero truth relations: + * + * A. The particle's SimTrackerHits never produced a truthed raw hit. The digitization + * dropped the relation -- e.g. the MC contribution failed the G1 threshold on a + * channel that a pulser hit kept alive (MERGED_SUBTHRESH). + * + * B. The particle's SimTrackerHits DID produce truthed raw hits, but pattern + * recognition built the track out of different raw hits instead. + * + * The distinguishing measurement is per SimTrackerHit: walk SVTTrueHitRelations backwards + * to the raw hits that carry it, then ask whether those raw hits are on the track. A is + * "no raw hit carries this SimTrackerHit"; B is "one does, but the track did not use it". + * + * Requires the readout to have been run with writeHitOriginCollections=true. + * + * Configure with one or more entries; with none, every event is dumped, so + * always set at least one. + */ +public class SvtEventForensicsDriver extends Driver { + + private String trackCollectionName = "KalmanFullTracks"; + private String truthRelationCollectionName = "SVTTrueHitRelations"; + private String pulserOriginCollectionName = "SVTHitOriginPulser"; + private String mcContribCollectionName = "SVTHitOriginMCContrib"; + private String simHitCollectionName = "TrackerHits"; + private String rawHitCollectionName = "SVTRawTrackerHits"; + private String clusterCollectionName = "StripClusterer_SiTrackerHitStrip1D"; + private String[] vertexCollectionNames = { "UnconstrainedV0Vertices_KF" }; + + /** Half-width, in strips, of the neighbourhood listed around a lost hit. */ + private int neighbourWindow = 12; + + /** Emit one "RAWHIT event sensor channel category" line per raw hit, for diffing arms. */ + private boolean dumpAllRawHits = false; + + /** Only dump MCParticles above this momentum, to keep the listing readable. [GeV] */ + private double mcpMinMomentum = 0.05; + + private final Set eventNumbers = new HashSet(); + + private static final int NOISE = 0, MC_PURE = 1, MC_PURE_SUBTHRESH = 2, + PULSER_PURE = 3, MERGED = 4, MERGED_SUBTHRESH = 5; + private static final String[] CAT_NAME = { + "NOISE", "MC_PURE", "MC_PURE_SUBTHRESH", "PULSER_PURE", "MERGED", "MERGED_SUBTHRESH" + }; + + public void setTrackCollectionName(String val) { this.trackCollectionName = val; } + public void setTruthRelationCollectionName(String val) { this.truthRelationCollectionName = val; } + public void setPulserOriginCollectionName(String val) { this.pulserOriginCollectionName = val; } + public void setMcContribCollectionName(String val) { this.mcContribCollectionName = val; } + public void setSimHitCollectionName(String val) { this.simHitCollectionName = val; } + public void setRawHitCollectionName(String val) { this.rawHitCollectionName = val; } + public void setClusterCollectionName(String val) { this.clusterCollectionName = val; } + public void setNeighbourWindow(int val) { this.neighbourWindow = val; } + public void setDumpAllRawHits(boolean val) { this.dumpAllRawHits = val; } + public void setVertexCollectionNames(String val) { this.vertexCollectionNames = val.split(" +"); } + public void setMcpMinMomentum(double val) { this.mcpMinMomentum = val; } + + /** Add one event number to dump. Repeat the element to dump several. */ + public void setEventNumber(int val) { this.eventNumbers.add(val); } + + private Set fromSide(EventHeader event, String name) { + Set out = new HashSet(); + if(!event.hasCollection(LCRelation.class, name)) { return out; } + for(LCRelation rel : event.get(LCRelation.class, name)) { + if(rel.getFrom() instanceof RawTrackerHit) { out.add((RawTrackerHit) rel.getFrom()); } + } + return out; + } + + private static String sensorOf(RawTrackerHit hit) { + return hit.getDetectorElement().getName(); + } + + private static int channelOf(RawTrackerHit hit) { + return hit.getIdentifierFieldValue("strip"); + } + + private static int categorise(RawTrackerHit hit, Set pulser, + Set mcContrib, Set truthed) { + boolean p = pulser.contains(hit); + boolean m = mcContrib.contains(hit); + boolean t = truthed.contains(hit); + if(!m) { return p ? PULSER_PURE : NOISE; } + if(p) { return t ? MERGED : MERGED_SUBTHRESH; } + return t ? MC_PURE : MC_PURE_SUBTHRESH; + } + + private static double momentum(MCParticle p) { + return p.getMomentum().magnitude(); + } + + /** tan(lambda) for an MCParticle, matching the track-parameter convention. */ + private static double tanLambda(MCParticle p) { + Hep3Vector m = p.getMomentum(); + double pt = Math.hypot(m.x(), m.z()); + return pt > 0 ? m.y() / pt : Double.NaN; + } + + @Override + protected void process(EventHeader event) { + if(!eventNumbers.isEmpty() && !eventNumbers.contains(event.getEventNumber())) { return; } + + Set pulser = fromSide(event, pulserOriginCollectionName); + Set mcContrib = fromSide(event, mcContribCollectionName); + Set truthed = fromSide(event, truthRelationCollectionName); + + // raw hit -> the SimTrackerHits it carries, and the reverse + Map> rawToSim = new HashMap>(); + Map> simToRaw = new HashMap>(); + if(event.hasCollection(LCRelation.class, truthRelationCollectionName)) { + for(LCRelation rel : event.get(LCRelation.class, truthRelationCollectionName)) { + if(!(rel.getFrom() instanceof RawTrackerHit)) { continue; } + if(!(rel.getTo() instanceof SimTrackerHit)) { continue; } + RawTrackerHit r = (RawTrackerHit) rel.getFrom(); + SimTrackerHit s = (SimTrackerHit) rel.getTo(); + if(!rawToSim.containsKey(r)) { rawToSim.put(r, new ArrayList()); } + rawToSim.get(r).add(s); + if(!simToRaw.containsKey(s)) { simToRaw.put(s, new ArrayList()); } + simToRaw.get(s).add(r); + } + } + + int nRaw = event.hasCollection(RawTrackerHit.class, rawHitCollectionName) + ? event.get(RawTrackerHit.class, rawHitCollectionName).size() : -1; + + System.out.println(); + System.out.println("################ FORENSICS run " + event.getRunNumber() + + " event " + event.getEventNumber() + " ################"); + System.out.println(" raw hits in event : " + nRaw); + System.out.println(" with truth relation : " + truthed.size()); + System.out.println(" with pulser origin : " + pulser.size()); + System.out.println(" with MC contribution : " + mcContrib.size()); + + // ---- machine-readable dump of every raw hit ---- + // Emitted so the two arms can be diffed channel by channel. The question it + // answers: are the hits labelled PULSER_PURE really independent data, or are + // they a second copy of the MC deposit at shifted channels? A PULSER_PURE + // channel that also fires in the no-pulser arm cannot have come from the + // pulser file, since that arm has no pulser file. + if(dumpAllRawHits && event.hasCollection(RawTrackerHit.class, rawHitCollectionName)) { + for(RawTrackerHit r : event.get(RawTrackerHit.class, rawHitCollectionName)) { + System.out.printf("RAWHIT %d %s %d %s%n", event.getEventNumber(), sensorOf(r), + channelOf(r), CAT_NAME[categorise(r, pulser, mcContrib, truthed)]); + } + } + + // ---- MCParticles of interest ---- + System.out.println(" -- MCParticles (|pdg|=11, p > " + mcpMinMomentum + ") --"); + List mcps = event.getMCParticles(); + for(MCParticle p : mcps) { + if(Math.abs(p.getPDGID()) != 11 || momentum(p) < mcpMinMomentum) { continue; } + System.out.printf(" pdg=%+3d p=%.4f tanL=%+.5f vtx_z=%+.3f%n", + p.getPDGID(), momentum(p), tanLambda(p), p.getOriginZ()); + } + + // ---- SimTrackerHits grouped by particle ---- + Map> simByParticle = new HashMap>(); + if(event.hasCollection(SimTrackerHit.class, simHitCollectionName)) { + for(SimTrackerHit s : event.get(SimTrackerHit.class, simHitCollectionName)) { + MCParticle p = s.getMCParticle(); + if(p == null) { continue; } + if(!simByParticle.containsKey(p)) { simByParticle.put(p, new ArrayList()); } + simByParticle.get(p).add(s); + } + } + + // ---- raw hit -> the cluster containing it ---- + // Tracking never sees raw hits; it sees the clusters TrackerHitDriver builds. A + // raw hit that is present and truthed but is in no cluster, or in a cluster no + // track used, has been removed from tracking's view before any pattern + // recognition happens -- which is a different failure from losing a + // combinatorial competition. + Map rawToCluster = new HashMap(); + List clusters = event.hasCollection(TrackerHit.class, clusterCollectionName) + ? event.get(TrackerHit.class, clusterCollectionName) + : new ArrayList(); + for(TrackerHit c : clusters) { + for(Object o : c.getRawHits()) { + if(o instanceof RawTrackerHit) { rawToCluster.put((RawTrackerHit) o, c); } + } + } + + // sensor -> every raw hit on it, so the neighbourhood of a lost hit can be shown + Map> bySensor = new HashMap>(); + if(event.hasCollection(RawTrackerHit.class, rawHitCollectionName)) { + for(RawTrackerHit r : event.get(RawTrackerHit.class, rawHitCollectionName)) { + String s = sensorOf(r); + if(!bySensor.containsKey(s)) { bySensor.put(s, new ArrayList()); } + bySensor.get(s).add(r); + } + } + + // ---- raw hits used by each track ---- + Map hitToTrack = new HashMap(); + List tracks = event.hasCollection(Track.class, trackCollectionName) + ? event.get(Track.class, trackCollectionName) : new ArrayList(); + for(int it = 0; it < tracks.size(); it++) { + for(TrackerHit th : tracks.get(it).getTrackerHits()) { + for(Object o : th.getRawHits()) { + if(o instanceof RawTrackerHit) { hitToTrack.put((RawTrackerHit) o, it); } + } + } + } + + // ---- per-track hit listing ---- + System.out.println(" -- tracks in " + trackCollectionName + ": " + tracks.size() + " --"); + for(int it = 0; it < tracks.size(); it++) { + Track t = tracks.get(it); + TrackState ts = t.getTrackStates().isEmpty() ? null : t.getTrackStates().get(0); + double[] mom = ts != null ? ts.getMomentum() : new double[]{0, 0, 0}; + double pmag = Math.sqrt(mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2]); + double tanL = ts != null ? ts.getTanLambda() : Double.NaN; + System.out.printf(" track %d: charge=%+d p=%.4f tanL=%+.5f chi2=%.2f ndf=%d nClusters=%d%n", + it, (int) t.getCharge(), pmag, tanL, t.getChi2(), t.getNDF(), + t.getTrackerHits().size()); + int nTruthed = 0, nHits = 0; + int[] byCat = new int[6]; + for(TrackerHit th : t.getTrackerHits()) { + for(Object o : th.getRawHits()) { + if(!(o instanceof RawTrackerHit)) { continue; } + RawTrackerHit r = (RawTrackerHit) o; + nHits++; + int cat = categorise(r, pulser, mcContrib, truthed); + byCat[cat]++; + StringBuilder who = new StringBuilder(); + if(rawToSim.containsKey(r)) { + nTruthed++; + for(SimTrackerHit s : rawToSim.get(r)) { + MCParticle p = s.getMCParticle(); + who.append(" <- pdg=").append(p == null ? "null" : p.getPDGID()) + .append(" p=").append(p == null ? "?" : String.format("%.3f", momentum(p))); + } + } + System.out.printf(" %-28s ch=%4d %-17s%s%n", + sensorOf(r), channelOf(r), CAT_NAME[cat], who); + } + } + System.out.printf(" -> %d raw hits, %d truthed;", nHits, nTruthed); + for(int c = 0; c < 6; c++) { if(byCat[c] > 0) { System.out.print(" " + CAT_NAME[c] + "=" + byCat[c]); } } + System.out.println(); + } + + // ---- vertices, with the provenance of the tracks they actually reference ---- + // Identifying a vertex's track by matching momenta against the track collection + // is unreliable: the momentum a vertex reports is the fitted one, and hpstr reads + // a different track state again. Walk the references instead. + for(String vtxColl : vertexCollectionNames) { + if(!event.hasCollection(Vertex.class, vtxColl)) { continue; } + List vertices = event.get(Vertex.class, vtxColl); + System.out.println(" -- " + vtxColl + ": " + vertices.size() + " --"); + int iv = 0; + for(Vertex v : vertices) { + System.out.printf(" vertex %d: chi2=%.2f pos=(%.3f, %.3f, %.3f) invM=%.4f%n", + iv++, v.getChi2(), v.getPosition().x(), v.getPosition().y(), + v.getPosition().z(), + v.getAssociatedParticle() == null ? Double.NaN + : v.getAssociatedParticle().getMass()); + if(v.getAssociatedParticle() == null) { continue; } + for(ReconstructedParticle rp : v.getAssociatedParticle().getParticles()) { + for(Track t : rp.getTracks()) { + int nHits = 0, nTruthed = 0; + int[] byCat = new int[6]; + Set carried = new HashSet(); + for(TrackerHit th : t.getTrackerHits()) { + for(Object o : th.getRawHits()) { + if(!(o instanceof RawTrackerHit)) { continue; } + RawTrackerHit r = (RawTrackerHit) o; + nHits++; + byCat[categorise(r, pulser, mcContrib, truthed)]++; + if(rawToSim.containsKey(r)) { + nTruthed++; + for(SimTrackerHit s : rawToSim.get(r)) { + if(s.getMCParticle() != null) { carried.add(s.getMCParticle()); } + } + } + } + } + TrackState ts = t.getTrackStates().isEmpty() ? null : t.getTrackStates().get(0); + double[] m = ts != null ? ts.getMomentum() : new double[]{0, 0, 0}; + Integer which = null; + for(int it = 0; it < tracks.size(); it++) { + if(tracks.get(it) == t) { which = it; } + } + System.out.printf(" rp charge=%+d rp_p=%.4f | track#%s p=%.4f tanL=%+.5f " + + "chi2=%.2f nRawHits=%d nTruthed=%d", + (int) rp.getCharge(), rp.getMomentum().magnitude(), + which == null ? "NOT-IN-" + trackCollectionName : which.toString(), + Math.sqrt(m[0]*m[0] + m[1]*m[1] + m[2]*m[2]), + ts != null ? ts.getTanLambda() : Double.NaN, + t.getChi2(), nHits, nTruthed); + for(int c = 0; c < 6; c++) { + if(byCat[c] > 0) { System.out.print(" " + CAT_NAME[c] + "=" + byCat[c]); } + } + for(MCParticle p : carried) { + System.out.printf(" <- pdg=%+d p=%.3f", p.getPDGID(), momentum(p)); + } + System.out.println(); + } + } + } + } + + // ---- the decisive part: what became of each MC hit ---- + System.out.println(" -- fate of every SimTrackerHit of each e+- --"); + List ordered = new ArrayList(simByParticle.keySet()); + Collections.sort(ordered, (a, b) -> Double.compare(momentum(b), momentum(a))); + for(MCParticle p : ordered) { + if(Math.abs(p.getPDGID()) != 11 || momentum(p) < mcpMinMomentum) { continue; } + List sims = simByParticle.get(p); + System.out.printf(" particle pdg=%+d p=%.4f tanL=%+.5f : %d SimTrackerHits%n", + p.getPDGID(), momentum(p), tanLambda(p), sims.size()); + int lost = 0, onTrack = 0, orphan = 0, noCluster = 0; + for(SimTrackerHit s : sims) { + List raws = simToRaw.get(s); + if(raws == null || raws.isEmpty()) { + lost++; + System.out.printf(" %-28s NO RAW HIT CARRIES THIS -- truth dropped in digitization%n", + s.getDetectorElement().getName()); + continue; + } + for(RawTrackerHit r : raws) { + Integer trk = hitToTrack.get(r); + if(trk == null) { orphan++; } else { onTrack++; } + TrackerHit c = rawToCluster.get(r); + String clus; + if(c == null) { + clus = "NOT IN ANY CLUSTER"; + noCluster++; + } else { + boolean used = false; + for(Object o : c.getRawHits()) { + if(o instanceof RawTrackerHit && hitToTrack.containsKey(o)) { used = true; } + } + clus = String.format("cluster size=%d t=%.2f %s", + c.getRawHits().size(), c.getTime(), + used ? "USED by a track" : "used by NO track"); + } + System.out.printf(" %-28s ch=%4d %-17s %-16s %s%n", + sensorOf(r), channelOf(r), + CAT_NAME[categorise(r, pulser, mcContrib, truthed)], + trk == null ? "NOT on any track" : ("on track " + trk), + clus); + // What else is on this sensor nearby? A pulser hit on an adjacent + // channel merges into the same cluster and can move or spoil it. + List near = bySensor.get(sensorOf(r)); + if(near != null) { + StringBuilder nb = new StringBuilder(); + for(RawTrackerHit o : near) { + int d = Math.abs(channelOf(o) - channelOf(r)); + if(o != r && d <= neighbourWindow) { + nb.append(String.format(" ch=%d(%s,d=%d)", channelOf(o), + CAT_NAME[categorise(o, pulser, mcContrib, truthed)], d)); + } + } + if(nb.length() > 0) { System.out.println(" neighbours:" + nb); } + } + } + } + System.out.printf(" -> %d SimTrackerHits with no raw hit (case A), " + + "%d raw hits on a track, %d truthed raw hits used by no track (case B), " + + "%d of those in NO cluster at all (case C)%n", + lost, onTrack, orphan, noCluster); + } + System.out.println("################ end event " + event.getEventNumber() + " ################"); + } +} diff --git a/analysis/src/main/java/org/hps/analysis/MC/SvtHitProvenanceDriver.java b/analysis/src/main/java/org/hps/analysis/MC/SvtHitProvenanceDriver.java new file mode 100644 index 000000000..feccedad1 --- /dev/null +++ b/analysis/src/main/java/org/hps/analysis/MC/SvtHitProvenanceDriver.java @@ -0,0 +1,297 @@ +package org.hps.analysis.MC; + +import java.util.ArrayList; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.TreeSet; + +import org.lcsim.event.EventHeader; +import org.lcsim.event.LCRelation; +import org.lcsim.event.RawTrackerHit; +import org.lcsim.event.Track; +import org.lcsim.event.TrackerHit; +import org.lcsim.util.Driver; + +/** + * Attributes the hits on each track to their origin, using the provenance relations + * written by SvtDigitizationWithPulserDataMergingReadoutDriver. + * + * Under pulser overlay a large fraction of tracks carry no truth relation at all. Two + * explanations survive the aggregate hit counts: + * + * 1. The tracks are fakes assembled from pulser data hits. + * 2. The tracks follow a real trajectory but were built from pulser hits on strips + * neighbouring the ones the MC particle actually hit. Adjacent strips are ~55 um + * apart, which over a metre of lever arm is ~55 urad, so such a track would still + * point at the true particle to well inside a milliradian while carrying zero + * truth relations. + * + * These differ observably in whether the untruthed hits on a track sit next to channels + * that did receive MC charge. That is what this driver measures, against the baseline + * rate at which any pulser hit in the event happens to be near an MC channel -- with + * high occupancy, adjacency alone proves nothing, so the comparison to the baseline is + * the whole point. + * + * Requires writeHitOriginCollections=true on the digitization driver. + */ +public class SvtHitProvenanceDriver extends Driver { + + private String trackCollectionName = "KalmanFullTracks"; + private String rawHitCollectionName = "SVTRawTrackerHits"; + private String truthRelationCollectionName = "SVTTrueHitRelations"; + private String pulserOriginCollectionName = "SVTHitOriginPulser"; + private String mcContribCollectionName = "SVTHitOriginMCContrib"; + + /** Neighbour distances, in strips, at which adjacency is reported. */ + private static final int[] DISTANCES = { 1, 2, 3, 5, 10 }; + /** Bins for the minimum strip distance to an MC-contributing channel. */ + private static final int MAX_DIST_BIN = 12; + + private boolean debug = false; + private int debugMaxPrint = 20; + private int debugPrinted = 0; + + // Hit categories, matching the digitization driver. + private static final int NOISE = 0; + private static final int MC_PURE = 1; + private static final int MC_PURE_SUBTHRESH = 2; + private static final int PULSER_PURE = 3; + private static final int MERGED = 4; + private static final int MERGED_SUBTHRESH = 5; + private static final String[] CAT_NAME = { + "NOISE ", "MC_PURE ", "MC_PURE_SUBTHRESH", + "PULSER_PURE ", "MERGED ", "MERGED_SUBTHRESH " + }; + + private long nEvents = 0; + private long nTracks = 0; + private long nZeroTruthTracks = 0; + private boolean warnedMissing = false; + + // [0] = tracks with at least one truthed hit, [1] = tracks with none + private final long[] nTracksByClass = new long[2]; + private final long[] nHitsByClass = new long[2]; + private final long[][] nHitsByClassAndCat = new long[2][6]; + + // Minimum strip distance from a PULSER_PURE hit to an MC-contributing channel on the + // same sensor. Index MAX_DIST_BIN is the overflow, MAX_DIST_BIN+1 means the sensor had + // no MC-contributing channel at all. + private final long[][] distHistByClass = new long[2][MAX_DIST_BIN + 2]; + private final long[] distHistBaseline = new long[MAX_DIST_BIN + 2]; + private long nPulserPureBaseline = 0; + + public void setTrackCollectionName(String val) { this.trackCollectionName = val; } + public void setRawHitCollectionName(String val) { this.rawHitCollectionName = val; } + public void setTruthRelationCollectionName(String val) { this.truthRelationCollectionName = val; } + public void setPulserOriginCollectionName(String val) { this.pulserOriginCollectionName = val; } + public void setMcContribCollectionName(String val) { this.mcContribCollectionName = val; } + public void setDebug(boolean val) { this.debug = val; } + public void setDebugMaxPrint(int val) { this.debugMaxPrint = val; } + + /** Collects the "from" side of a relation collection, tolerating a missing collection. */ + private Set fromSide(EventHeader event, String name) { + Set out = new HashSet(); + if(!event.hasCollection(LCRelation.class, name)) { return out; } + for(LCRelation rel : event.get(LCRelation.class, name)) { + if(rel.getFrom() instanceof RawTrackerHit) { + out.add((RawTrackerHit) rel.getFrom()); + } + } + return out; + } + + private static String sensorOf(RawTrackerHit hit) { + return hit.getDetectorElement().getName(); + } + + private static int channelOf(RawTrackerHit hit) { + return hit.getIdentifierFieldValue("strip"); + } + + private static int categorise(RawTrackerHit hit, Set pulser, + Set mcContrib, Set truthed) { + boolean p = pulser.contains(hit); + boolean m = mcContrib.contains(hit); + boolean t = truthed.contains(hit); + if(!m) { return p ? PULSER_PURE : NOISE; } + if(p) { return t ? MERGED : MERGED_SUBTHRESH; } + return t ? MC_PURE : MC_PURE_SUBTHRESH; + } + + /** + * Minimum distance in strips from this hit to a channel that received MC charge on the + * same sensor. Returns MAX_DIST_BIN+1 if the sensor had no MC contribution anywhere. + */ + private int minDistanceToMC(RawTrackerHit hit, Map> mcChannels) { + TreeSet chans = mcChannels.get(sensorOf(hit)); + if(chans == null || chans.isEmpty()) { return MAX_DIST_BIN + 1; } + int ch = channelOf(hit); + Integer lo = chans.floor(ch); + Integer hi = chans.ceiling(ch); + int best = Integer.MAX_VALUE; + if(lo != null) { best = Math.min(best, ch - lo); } + if(hi != null) { best = Math.min(best, hi - ch); } + if(best == Integer.MAX_VALUE) { return MAX_DIST_BIN + 1; } + return Math.min(best, MAX_DIST_BIN); + } + + @Override + public void process(EventHeader event) { + if(!event.hasCollection(Track.class, trackCollectionName)) { return; } + + Set truthed = fromSide(event, truthRelationCollectionName); + Set pulser = fromSide(event, pulserOriginCollectionName); + Set mcContrib = fromSide(event, mcContribCollectionName); + + if(!warnedMissing && !event.hasCollection(LCRelation.class, mcContribCollectionName)) { + warnedMissing = true; + System.out.println("SvtHitProvenanceDriver: WARNING collection '" + mcContribCollectionName + + "' not found. Was writeHitOriginCollections set on the digitization driver?"); + } + + nEvents++; + + // Channels that received MC charge, by sensor. Built from the MC-contribution + // relation so it is independent of whether the truth gate kept the relation. + Map> mcChannels = new HashMap>(); + for(RawTrackerHit hit : mcContrib) { + String s = sensorOf(hit); + TreeSet set = mcChannels.get(s); + if(set == null) { set = new TreeSet(); mcChannels.put(s, set); } + set.add(channelOf(hit)); + } + + // Baseline: how near an MC channel does an arbitrary pulser hit in this event sit? + // Tracks are compared against this, not against zero. + if(event.hasCollection(RawTrackerHit.class, rawHitCollectionName)) { + for(RawTrackerHit hit : event.get(RawTrackerHit.class, rawHitCollectionName)) { + if(categorise(hit, pulser, mcContrib, truthed) != PULSER_PURE) { continue; } + nPulserPureBaseline++; + distHistBaseline[minDistanceToMC(hit, mcChannels)]++; + } + } + + for(Track track : event.get(Track.class, trackCollectionName)) { + nTracks++; + + List rawHits = new ArrayList(); + for(TrackerHit th : track.getTrackerHits()) { + for(Object o : th.getRawHits()) { + if(o instanceof RawTrackerHit) { rawHits.add((RawTrackerHit) o); } + } + } + if(rawHits.isEmpty()) { continue; } + + int nTruthedOnTrack = 0; + for(RawTrackerHit hit : rawHits) { + if(truthed.contains(hit)) { nTruthedOnTrack++; } + } + final int cls = (nTruthedOnTrack == 0) ? 1 : 0; + if(cls == 1) { nZeroTruthTracks++; } + nTracksByClass[cls]++; + nHitsByClass[cls] += rawHits.size(); + + for(RawTrackerHit hit : rawHits) { + int cat = categorise(hit, pulser, mcContrib, truthed); + nHitsByClassAndCat[cls][cat]++; + if(cat == PULSER_PURE) { + distHistByClass[cls][minDistanceToMC(hit, mcChannels)]++; + } + } + + if(debug && cls == 1 && debugPrinted < debugMaxPrint) { + debugPrinted++; + StringBuilder sb = new StringBuilder("[SvtProv] zero-truth track, nRawHits=" + + rawHits.size() + " hits:"); + for(RawTrackerHit hit : rawHits) { + int cat = categorise(hit, pulser, mcContrib, truthed); + sb.append(" ").append(CAT_NAME[cat].trim()) + .append("(").append(sensorOf(hit)).append(":").append(channelOf(hit)); + if(cat == PULSER_PURE) { + int d = minDistanceToMC(hit, mcChannels); + sb.append(", dMC=").append(d > MAX_DIST_BIN ? "none" : Integer.toString(d)); + } + sb.append(")"); + } + System.out.println(sb.toString()); + } + } + } + + /** Fraction of entries in a distance histogram at or below d strips. */ + private static double fracWithin(long[] hist, int d) { + long num = 0, den = 0; + for(int i = 0; i < hist.length; i++) { + den += hist[i]; + if(i <= d) { num += hist[i]; } + } + return den > 0 ? (double) num / den : 0.0; + } + + @Override + public void endOfData() { + String[] clsName = { "tracks with truth ", "ZERO-truth tracks " }; + System.out.println(); + System.out.println("============== SVT hit provenance by track =============="); + System.out.println(" events : " + nEvents); + System.out.println(" tracks (" + trackCollectionName + ") : " + nTracks); + System.out.println(" zero-truth tracks : " + nZeroTruthTracks + + (nTracks > 0 ? String.format(" (%.4f)", (double) nZeroTruthTracks / nTracks) : "")); + System.out.println(); + + for(int cls = 0; cls < 2; cls++) { + System.out.println(" ---- " + clsName[cls] + " ----"); + System.out.println(" tracks : " + nTracksByClass[cls] + + " raw hits/track : " + + (nTracksByClass[cls] > 0 + ? String.format("%.2f", (double) nHitsByClass[cls] / nTracksByClass[cls]) : "-")); + for(int cat = 0; cat < 6; cat++) { + long n = nHitsByClassAndCat[cls][cat]; + System.out.println(" " + CAT_NAME[cat] + " : " + n + + (nHitsByClass[cls] > 0 + ? String.format(" (%.4f)", (double) n / nHitsByClass[cls]) : "")); + } + System.out.println(); + } + + System.out.println(" ---- are PULSER_PURE hits next to channels that saw MC charge? ----"); + System.out.println(" The baseline is every PULSER_PURE hit in the event, so it already"); + System.out.println(" folds in the occupancy. Only an excess over baseline is meaningful."); + StringBuilder hdr = new StringBuilder(String.format(" %-26s %10s", "population", "nHits")); + for(int d : DISTANCES) { hdr.append(String.format(" within%3d", d)); } + System.out.println(hdr.toString()); + + long nZ = 0, nT = 0; + for(int i = 0; i < distHistByClass[1].length; i++) { nZ += distHistByClass[1][i]; } + for(int i = 0; i < distHistByClass[0].length; i++) { nT += distHistByClass[0][i]; } + + Object[][] rows = { + { "baseline (all in event)", distHistBaseline, nPulserPureBaseline }, + { "on tracks with truth", distHistByClass[0], nT }, + { "on ZERO-truth tracks", distHistByClass[1], nZ }, + }; + for(Object[] row : rows) { + StringBuilder sb = new StringBuilder(String.format(" %-26s %10d", row[0], (Long) row[2])); + for(int d : DISTANCES) { + sb.append(String.format(" %8.4f", fracWithin((long[]) row[1], d))); + } + System.out.println(sb.toString()); + } + + System.out.println(); + System.out.println(" minimum strip distance to an MC channel, ZERO-truth tracks:"); + StringBuilder sb = new StringBuilder(" "); + for(int i = 0; i <= MAX_DIST_BIN; i++) { + sb.append(i == MAX_DIST_BIN ? ">=" + MAX_DIST_BIN : Integer.toString(i)) + .append("=").append(distHistByClass[1][i]).append(" "); + } + sb.append("noMCOnSensor=").append(distHistByClass[1][MAX_DIST_BIN + 1]); + System.out.println(sb.toString()); + System.out.println("========================================================="); + System.out.println(); + super.endOfData(); + } +} diff --git a/detector-model/dependency-reduced-pom.xml b/detector-model/dependency-reduced-pom.xml index a57ace5df..6be2e5533 100644 --- a/detector-model/dependency-reduced-pom.xml +++ b/detector-model/dependency-reduced-pom.xml @@ -51,7 +51,7 @@ junit junit - 4.13.1 + 4.13.2 test diff --git a/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java b/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java index 41a69282e..cef66c006 100755 --- a/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java +++ b/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java @@ -73,10 +73,111 @@ public class SvtDigitizationWithPulserDataMergingReadoutDriver extends ReadoutDr private boolean enablePileupCut = true; private boolean dropBadChannels = true; + // ------------------------------------------------------------------ + // Truth-relation diagnostics. + // + // An MC strip hit only reaches SVTTrueHitRelations if it passes two + // independent gates in getOnTriggerData(): + // G1 totalContrib > 4.0 * meanNoise -- attaches hit.simHits to the channel + // G2 readoutCuts(hit) -- if this fails the raw hit is dropped + // and the relations are never built + // G1 sees only the MC contribution, so it cannot depend on the pulser overlay. + // G2 is evaluated on the combined waveform, and on a channel carrying a pulser + // hit that waveform is real data whose baseline need not agree with the pedestal + // that samplesAboveThreshold() subtracts. Every counter below is therefore + // indexed [0] = no pulser hit on this channel, [1] = pulser hit on this channel, + // so the two populations can be compared directly. + // ------------------------------------------------------------------ + private boolean debug = false; + private int debugMaxPrint = 200; + private int debugPrinted = 0; + + private long nChannelsSeen = 0; + private long nChannelsWithPulser = 0; + private long nPulserQueueEmptyNonNull = 0; // would make poll() return null at L668 + + private final long[] nMCStripHits = new long[2]; + private final long[] nTruthGatePass = new long[2]; + private final long[] nTruthGateFail = new long[2]; + private final long[] nSimHitsLostTruthGate = new long[2]; + + private final long[] nHitsReadoutPass = new long[2]; + private final long[] nHitsReadoutFail = new long[2]; + private final long[] nHitsReadoutFailWithTruth = new long[2]; + private final long[] nSimHitsLostReadoutCut = new long[2]; + private final long[] nRelationsWritten = new long[2]; + + // totalContrib / meanNoise, binned: <1, 1-2, 2-4, 4-8, 8-16, >=16 + private final long[][] truthRatioBins = new long[2][6]; + + // Mean offset of the pulser waveform from the conditions-DB pedestal, sampled + // before any MC contribution is added. A negative value means the overlaid data + // sits below the pedestal that the threshold cut subtracts, which eats into the + // MC signal's headroom over threshold. + private double sumPulserBaselineOffset = 0; + private double sumPulserBaselineOffsetSq = 0; + private long nPulserBaselineChannels = 0; + + //-------------------------// + //--- Hit origin labels ---// + //-------------------------// + + /** + * Opt-in provenance labelling for the output raw hits. + * + * The absence of a truth relation on a raw hit is ambiguous: it can mean the + * channel carried only overlaid pulser data, or that an MC particle did deposit + * charge there but fell below the truth gate G1 and so had its relation dropped. + * Those two cases need different fixes, so they are separated here by recording + * what actually contributed to each channel, independent of G1. + * + * Two extra relation collections are written: + * + * SVTHitOriginPulser (RawTrackerHit out) -> (RawTrackerHit pulser source) + * SVTHitOriginMCContrib (RawTrackerHit out) -> (MCParticle depositing charge) + * + * The MC-contribution relation is written whether or not G1 passed. Combined with + * the existing SVTTrueHitRelations this gives a complete category per hit: + * + * pulser mcContrib truthRel category + * - - - NOISE pedestal and noise only + * - Y Y MC_PURE MC only, truth kept + * - Y - MC_PURE_SUBTHRESH MC only, truth lost at G1 + * Y - - PULSER_PURE overlaid data only + * Y Y Y MERGED MC + data, truth kept + * Y Y - MERGED_SUBTHRESH MC + data, truth lost at G1 + * + * MERGED_SUBTHRESH is the category that distinguishes a genuine relation-propagation + * bug in the merging from tracks simply picking up real pulser hits. + * + * Off by default: when enabled the MCParticle collection is widened to cover the + * sub-threshold contributors, so the production output is only bit-identical with + * this disabled. + */ + private boolean writeHitOriginCollections = false; + private String hitOriginPulserCollection = "SVTHitOriginPulser"; + private String hitOriginMCContribCollection = "SVTHitOriginMCContrib"; + + private LCIOCollection hitOriginPulserCollectionParams; + private LCIOCollection hitOriginMCContribCollectionParams; + + private static final int ORIGIN_NOISE = 0; + private static final int ORIGIN_MC_PURE = 1; + private static final int ORIGIN_MC_PURE_SUBTHRESH = 2; + private static final int ORIGIN_PULSER_PURE = 3; + private static final int ORIGIN_MERGED = 4; + private static final int ORIGIN_MERGED_SUBTHRESH = 5; + private static final String[] ORIGIN_NAME = { + "NOISE ", "MC_PURE ", "MC_PURE_SUBTHRESH", + "PULSER_PURE ", "MERGED ", "MERGED_SUBTHRESH " + }; + private final long[] nHitsByOrigin = new long[6]; + private final long[] nSimHitsByOrigin = new long[6]; + // Collection Names private String outputCollection = "SVTRawTrackerHits"; private String relationCollection = "SVTTrueHitRelations"; - + private LCIOCollection trackerHitCollectionParams; private LCIOCollection truthRelationsCollectionParams; private LCIOCollection truthHitsCollectionParams; @@ -97,6 +198,41 @@ public SvtDigitizationWithPulserDataMergingReadoutDriver() { * analog hits, while false uses only contributions * from pulses generated from truth data. */ + /** + * Enables the truth-relation diagnostics. Counters are accumulated regardless; + * this additionally prints per-hit lines for the first {@link #debugMaxPrint} + * MC strip hits that lose their truth relation. + * @param debug - true to enable diagnostic printout. + */ + public void setDebug(boolean debug) { + this.debug = debug; + } + + /** + * Sets how many individual truth-loss records are printed when debug is enabled. + * @param debugMaxPrint - Maximum number of per-hit lines. + */ + public void setDebugMaxPrint(int debugMaxPrint) { + this.debugMaxPrint = debugMaxPrint; + } + + /** + * Write the per-hit provenance relation collections described above. Off by + * default; enabling it also widens the output MCParticle collection to cover + * sub-threshold contributors. + */ + public void setWriteHitOriginCollections(boolean writeHitOriginCollections) { + this.writeHitOriginCollections = writeHitOriginCollections; + } + + public void setHitOriginPulserCollection(String name) { + this.hitOriginPulserCollection = name; + } + + public void setHitOriginMCContribCollection(String name) { + this.hitOriginMCContribCollection = name; + } + public void setAddNoise(boolean addNoise) { this.addNoise = addNoise; } @@ -425,6 +561,16 @@ public void startOfData() { LCIOCollectionFactory.setCollectionName(relationCollection); LCIOCollectionFactory.setProductionDriver(this); truthRelationsCollectionParams = LCIOCollectionFactory.produceLCIOCollection(LCRelation.class); + + if(writeHitOriginCollections) { + LCIOCollectionFactory.setCollectionName(hitOriginPulserCollection); + LCIOCollectionFactory.setProductionDriver(this); + hitOriginPulserCollectionParams = LCIOCollectionFactory.produceLCIOCollection(LCRelation.class); + + LCIOCollectionFactory.setCollectionName(hitOriginMCContribCollection); + LCIOCollectionFactory.setProductionDriver(this); + hitOriginMCContribCollectionParams = LCIOCollectionFactory.produceLCIOCollection(LCRelation.class); + } LCIOCollectionFactory.setCollectionName("TrackerHits"); LCIOCollectionFactory.setFlags(0xc0000000); @@ -637,7 +783,14 @@ protected Collection> getOnTriggerData(double triggerTime) // Create a list to hold the analog data List hits = new ArrayList(); List truthHits = new ArrayList(); - List trueHitRelations = new ArrayList(); + List trueHitRelations = new ArrayList(); + + // Provenance relations, only populated when writeHitOriginCollections is set. + List hitOriginPulserRelations = new ArrayList(); + List hitOriginMCContribRelations = new ArrayList(); + // Contributors that failed G1 and so appear in no truth hit. They still need to + // reach the output MCParticle collection or the relations above would dangle. + Set extraContribParticles = new java.util.HashSet(); // Calculate time of first sample double firstSample = Math.floor(((triggerTime + 256) - readoutLatency - readoutOffset) / HPSSVTConstants.SAMPLING_INTERVAL) @@ -662,11 +815,22 @@ protected Collection> getOnTriggerData(double triggerTime) // the channel. double[] signal = new double[6]; + nChannelsSeen++; + if(pulserHitQueues[channel] != null && pulserHitQueues[channel].isEmpty()) { + // poll() below returns null on an empty queue, so this counts how often + // the != null test at the next line is not sufficient on its own. + nPulserQueueEmptyNonNull++; + } + //do the pulser hit first...if there is a pulser hit, don't add pedestal or noise to mc hit boolean hasPulserHit=false; // flag if this channel has a pulser hit + // Kept in scope so the output hit can be related back to the data hit + // it was merged with. + RawTrackerHit pulserSourceHit = null; if(pulserHitQueues[channel] != null){ StripHit ph=pulserHitQueues[channel].poll(); RawTrackerHit rth=ph.getRawTrackerHit(); + pulserSourceHit = rth; hasPulserHit=true; short[] samples =rth.getADCValues(); for(int sampleN = 0; sampleN < 6; sampleN++) { @@ -674,6 +838,23 @@ protected Collection> getOnTriggerData(double triggerTime) } } + // Index every counter by whether this channel carries pulser data. + final int pi = hasPulserHit ? 1 : 0; + if(hasPulserHit) { + nChannelsWithPulser++; + // Record where the overlaid data baseline sits relative to the pedestal + // that samplesAboveThreshold() will subtract. Sampled here, before the + // MC contribution is added below. + double offset = 0; + for(int sampleN = 0; sampleN < 6; sampleN++) { + offset += signal[sampleN] - ((HpsSiSensor) sensor).getPedestal(channel, sampleN); + } + offset /= 6; + sumPulserBaselineOffset += offset; + sumPulserBaselineOffsetSq += offset * offset; + nPulserBaselineChannels++; + } + if(!hasPulserHit){ // Create a buffer to hold the extracted signal for // the channel. Populate it with the appropriate @@ -690,12 +871,16 @@ protected Collection> getOnTriggerData(double triggerTime) // Create a list to store truth SVT hits. List simHits = new ArrayList(); - + // Every sim hit that deposited charge on this channel, whether or not it + // passed G1. simHits above is the G1-surviving subset. + List allContribSimHits = new ArrayList(); + // If there is data in the mc hit queues, process it. if(hitQueues[channel] != null) { for(StripHit hit : hitQueues[channel]) { processedHits.add(hit); - + allContribSimHits.addAll(hit.simHits); + // Track the noise and contribution to the // signal from the current hit. double meanNoise = 0; @@ -727,8 +912,35 @@ protected Collection> getOnTriggerData(double triggerTime) // from the hit. If it exceeds a the noise // threshold, store it as a truth hit. //meanNoise /= 6; + // ---- G1 accounting (no behaviour change) ---- + nMCStripHits[pi]++; + double truthRatio = (meanNoise > 0) ? (totalContrib / meanNoise) : -1; + if(truthRatio >= 0) { + int b; + if(truthRatio < 1) { b = 0; } + else if(truthRatio < 2) { b = 1; } + else if(truthRatio < 4) { b = 2; } + else if(truthRatio < 8) { b = 3; } + else if(truthRatio < 16) { b = 4; } + else { b = 5; } + truthRatioBins[pi][b]++; + } + if(totalContrib > 4.0 * meanNoise) { simHits.addAll(hit.simHits); + nTruthGatePass[pi]++; + } else { + nTruthGateFail[pi]++; + nSimHitsLostTruthGate[pi] += hit.simHits.size(); + if(debug && debugPrinted < debugMaxPrint) { + debugPrinted++; + System.out.println("[SvtDigiTruth] G1 FAIL pulser=" + hasPulserHit + + " sensor=" + sensor.getName() + " ch=" + channel + + " totalContrib=" + totalContrib + + " meanNoise=" + meanNoise + + " ratio=" + truthRatio + + " nSimHits=" + hit.simHits.size()); + } } } } @@ -747,6 +959,7 @@ protected Collection> getOnTriggerData(double triggerTime) // Only tracker hits that pass the readout cuts may // be passed through to readout. if(readoutCuts(hit)) { + nHitsReadoutPass[pi]++; // Add the hit to the readout hits collection. hits.add(hit); // Associate the truth hits with the raw hit and @@ -755,6 +968,64 @@ protected Collection> getOnTriggerData(double triggerTime) LCRelation hitRelation = new BaseLCRelation(hit, simHit); trueHitRelations.add(hitRelation); truthHits.add(simHit); + nRelationsWritten[pi]++; + } + + // ---- provenance labelling ---- + // Classify by what actually contributed charge, independent of G1, + // so that "no truth relation" resolves into a specific cause. + final boolean mcContrib = !allContribSimHits.isEmpty(); + final boolean truthKept = !simHits.isEmpty(); + final int origin; + if(!mcContrib) { + origin = hasPulserHit ? ORIGIN_PULSER_PURE : ORIGIN_NOISE; + } else if(hasPulserHit) { + origin = truthKept ? ORIGIN_MERGED : ORIGIN_MERGED_SUBTHRESH; + } else { + origin = truthKept ? ORIGIN_MC_PURE : ORIGIN_MC_PURE_SUBTHRESH; + } + nHitsByOrigin[origin]++; + nSimHitsByOrigin[origin] += allContribSimHits.size(); + + if(writeHitOriginCollections) { + if(pulserSourceHit != null) { + hitOriginPulserRelations.add(new BaseLCRelation(hit, pulserSourceHit)); + } + // One relation per distinct contributing particle. Duplicates are + // dropped so the collection size counts particles, not sim hits. + Set seen = new java.util.HashSet(); + for(SimTrackerHit simHit : allContribSimHits) { + MCParticle p = simHit.getMCParticle(); + if(p == null || !seen.add(p)) { continue; } + hitOriginMCContribRelations.add(new BaseLCRelation(hit, p)); + ReadoutDataManager.addParticleParents(p, extraContribParticles); + } + } + } else { + // ---- G2 accounting. Truth attached here is discarded with the hit. ---- + nHitsReadoutFail[pi]++; + if(!simHits.isEmpty()) { + nHitsReadoutFailWithTruth[pi]++; + nSimHitsLostReadoutCut[pi] += simHits.size(); + if(debug && debugPrinted < debugMaxPrint) { + debugPrinted++; + double baseline = 0; + int nAbove = 0; + for(int sampleN = 0; sampleN < 6; sampleN++) { + double ped = ((HpsSiSensor) sensor).getPedestal(channel, sampleN); + double nse = ((HpsSiSensor) sensor).getNoise(channel, sampleN); + baseline += samples[sampleN] - ped; + if(samples[sampleN] - ped > nse * noiseThreshold) { nAbove++; } + } + baseline /= 6; + System.out.println("[SvtDigiTruth] G2 FAIL pulser=" + hasPulserHit + + " sensor=" + sensor.getName() + " ch=" + channel + + " meanSampleMinusPed=" + baseline + + " nSamplesAboveThresh=" + nAbove + + "/" + samplesAboveThreshold + + " badChannel=" + !badChannelCut(hit) + + " nSimHitsLost=" + simHits.size()); + } } } } @@ -776,7 +1047,12 @@ protected Collection> getOnTriggerData(double triggerTime) for(SimTrackerHit simHit : truthHits) { ReadoutDataManager.addParticleParents(simHit.getMCParticle(), truthParticles); } - + // Sub-threshold contributors are referenced by the provenance relations but have + // no truth hit, so they would otherwise be missing from the written collection. + if(writeHitOriginCollections) { + truthParticles.addAll(extraContribParticles); + } + // Create the truth MC particle collection. LCIOCollectionFactory.setCollectionName("MCParticle"); LCIOCollectionFactory.setProductionDriver(this); @@ -792,17 +1068,97 @@ protected Collection> getOnTriggerData(double triggerTime) timestampData.getData().add(timestamp); // Store them in a single collection. - Collection> eventOutput = new ArrayList>(5); + Collection> eventOutput = new ArrayList>(7); eventOutput.add(hitCollection); eventOutput.add(truthParticleData); eventOutput.add(truthHitCollection); eventOutput.add(truthRelationCollection); eventOutput.add(timestampData); - + + if(writeHitOriginCollections) { + TriggeredLCIOData originPulserData = + new TriggeredLCIOData(hitOriginPulserCollectionParams); + originPulserData.getData().addAll(hitOriginPulserRelations); + eventOutput.add(originPulserData); + + TriggeredLCIOData originMCContribData = + new TriggeredLCIOData(hitOriginMCContribCollectionParams); + originMCContribData.getData().addAll(hitOriginMCContribRelations); + eventOutput.add(originMCContribData); + } + // Return the event output. return eventOutput; } + @Override + public void endOfData() { + String[] tag = { "no-pulser", "pulser " }; + System.out.println(); + System.out.println("================ SvtDigitization truth-relation summary ================"); + System.out.println(" channels processed : " + nChannelsSeen); + System.out.println(" channels with pulser hit : " + nChannelsWithPulser + + fraction(nChannelsWithPulser, nChannelsSeen)); + System.out.println(" pulser queues non-null but empty (poll() -> null) : " + nPulserQueueEmptyNonNull); + if(nPulserBaselineChannels > 0) { + double mean = sumPulserBaselineOffset / nPulserBaselineChannels; + double var = sumPulserBaselineOffsetSq / nPulserBaselineChannels - mean * mean; + System.out.println(" pulser baseline - DB pedestal [ADC] : mean=" + mean + + " rms=" + (var > 0 ? Math.sqrt(var) : 0.0) + + " (negative eats MC headroom over threshold)"); + } + System.out.println(); + System.out.println(" G1 = truth gate (totalContrib > 4*meanNoise)"); + System.out.println(" G2 = readoutCuts on the combined waveform"); + for(int pi = 0; pi < 2; pi++) { + System.out.println(" ---- " + tag[pi] + " ----"); + System.out.println(" MC strip hits : " + nMCStripHits[pi]); + System.out.println(" G1 pass / fail : " + nTruthGatePass[pi] + " / " + nTruthGateFail[pi] + + fraction(nTruthGateFail[pi], nMCStripHits[pi]) + " fail"); + System.out.println(" sim hits lost at G1 : " + nSimHitsLostTruthGate[pi]); + System.out.println(" raw hits G2 pass / fail : " + nHitsReadoutPass[pi] + " / " + nHitsReadoutFail[pi]); + System.out.println(" G2 failures carrying truth : " + nHitsReadoutFailWithTruth[pi]); + System.out.println(" sim hits lost at G2 : " + nSimHitsLostReadoutCut[pi]); + System.out.println(" relations written : " + nRelationsWritten[pi]); + long lost = nSimHitsLostTruthGate[pi] + nSimHitsLostReadoutCut[pi]; + System.out.println(" total sim hits lost : " + lost + + fraction(lost, lost + nRelationsWritten[pi])); + StringBuilder sb = new StringBuilder(" totalContrib/meanNoise : "); + String[] edges = { "<1", "1-2", "2-4", "4-8", "8-16", ">=16" }; + for(int b = 0; b < 6; b++) { + sb.append(edges[b]).append("=").append(truthRatioBins[pi][b]).append(" "); + } + System.out.println(sb.toString()); + } + System.out.println(); + System.out.println(" ---- provenance of the raw hits that reached readout ----"); + long totOrigin = 0; + for(int o = 0; o < nHitsByOrigin.length; o++) { totOrigin += nHitsByOrigin[o]; } + for(int o = 0; o < nHitsByOrigin.length; o++) { + System.out.println(" " + ORIGIN_NAME[o] + " : " + nHitsByOrigin[o] + + fraction(nHitsByOrigin[o], totOrigin) + + " contributing sim hits = " + nSimHitsByOrigin[o]); + } + long untruthed = nHitsByOrigin[ORIGIN_NOISE] + nHitsByOrigin[ORIGIN_PULSER_PURE] + + nHitsByOrigin[ORIGIN_MC_PURE_SUBTHRESH] + nHitsByOrigin[ORIGIN_MERGED_SUBTHRESH]; + long mcLostTruth = nHitsByOrigin[ORIGIN_MC_PURE_SUBTHRESH] + nHitsByOrigin[ORIGIN_MERGED_SUBTHRESH]; + System.out.println(" hits with no truth relation : " + untruthed + + fraction(untruthed, totOrigin)); + System.out.println(" of which MC did contribute : " + mcLostTruth + + fraction(mcLostTruth, untruthed) + + " <- relation lost, not absent"); + System.out.println(" origin collections written : " + writeHitOriginCollections); + System.out.println("======================================================================="); + System.out.println(); + super.endOfData(); + } + + /** Formats a count as a parenthesised fraction of a total, or "" if the total is zero. */ + private static String fraction(long num, long den) { + if(den <= 0) { return ""; } + return String.format(" (%.4f)", (double) num / den); + } + /** * Class StripHit is responsible for storing several * parameters defining a simulated hit object. diff --git a/recon/src/main/java/org/hps/recon/filtering/UnbiasedTriggerFilterDriver.java b/recon/src/main/java/org/hps/recon/filtering/UnbiasedTriggerFilterDriver.java index 2511ec681..1b95fcf9d 100644 --- a/recon/src/main/java/org/hps/recon/filtering/UnbiasedTriggerFilterDriver.java +++ b/recon/src/main/java/org/hps/recon/filtering/UnbiasedTriggerFilterDriver.java @@ -1,5 +1,6 @@ package org.hps.recon.filtering; +import java.util.logging.Logger; import org.lcsim.event.EventHeader; import org.lcsim.event.GenericObject; import org.lcsim.util.Driver; @@ -10,33 +11,99 @@ /** * Keep pulser triggered events. Also keep EPICS events, and Scaler events. Drop all other events. + * + * Set debug="true" in the steering file to enable per-event logging. + * A summary of accepted/rejected counts is always printed at end-of-data. */ public class UnbiasedTriggerFilterDriver extends Driver { + private static final Logger LOGGER = Logger.getLogger(UnbiasedTriggerFilterDriver.class.getName()); + + // Configurable via steering file: UnbiasedTriggerFilterDrivertrue + private boolean debug = false; + + // Counters + private long nTotal = 0; + private long nEpics = 0; + private long nScaler = 0; + private long nNoTSBank = 0; + private long nPulser = 0; + private long nFaradayCup = 0; + private long nDropped = 0; + + public void setDebug(boolean debug) { + this.debug = debug; + } + + @Override + public void startOfData() { + System.out.println("[UnbiasedTriggerFilterDriver] Driver started. debug=" + debug); + } + public void process(EventHeader event) { + nTotal++; + int evNum = event.getEventNumber(); + // 1. keep all events with EPICS data (could also use event tag = 31): - if (EpicsData.read(event) != null) + if (EpicsData.read(event) != null) { + nEpics++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": ACCEPTED (EPICS data)"); return; + } // 2. keep all events with Scaler data: - if (ScalerData.read(event) != null) + if (ScalerData.read(event) != null) { + nScaler++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": ACCEPTED (Scaler data)"); return; + } // 3. drop event if it doesn't have a TriggerBank - if (!event.hasCollection(GenericObject.class, "TSBank")) + if (!event.hasCollection(GenericObject.class, "TSBank")) { + nNoTSBank++; + nDropped++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": REJECTED (no TSBank)"); throw new Driver.NextEventException(); + } - // 4. keep event if it was from a Pulser trigger: + // 4. keep event if it was from a Pulser or FaradayCup trigger: for (GenericObject gob : event.get(GenericObject.class, "TSBank")) { if (!(AbstractIntData.getTag(gob) == TSData2019.BANK_TAG)) continue; TSData2019 tsd = new TSData2019(gob); - if (tsd.isPulserTrigger() || tsd.isFaradayCupTrigger()) + if (tsd.isPulserTrigger()) { + nPulser++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": ACCEPTED (Pulser trigger)"); return; + } + if (tsd.isFaradayCupTrigger()) { + nFaradayCup++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": ACCEPTED (FaradayCup trigger)"); + return; + } } // 5. Else, drop event: + nDropped++; + if (debug) System.out.println("[UnbiasedTriggerFilterDriver] Event " + evNum + ": REJECTED (no Pulser/FaradayCup trigger in TSBank)"); throw new Driver.NextEventException(); } + + @Override + public void endOfData() { + long nAccepted = nTotal - nDropped; + double efficiency = (nTotal > 0) ? 100.0 * nAccepted / nTotal : 0.0; + System.out.println("========== UnbiasedTriggerFilterDriver Summary =========="); + System.out.println(" Total events processed : " + nTotal); + System.out.println(String.format(" Accepted : %d (%.2f%%)", nAccepted, efficiency)); + System.out.println(" -> EPICS : " + nEpics); + System.out.println(" -> Scaler : " + nScaler); + System.out.println(" -> Pulser trigger : " + nPulser); + System.out.println(" -> FaradayCup trigger: " + nFaradayCup); + System.out.println(String.format(" Rejected : %d (%.2f%%)", nDropped, 100.0 - efficiency)); + System.out.println(" -> No TSBank : " + nNoTSBank); + System.out.println(" -> Other triggers : " + (nDropped - nNoTSBank)); + System.out.println("========================================================="); + } } diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021MCRecon_KF_WithSpacing_pass6dev_lowp_physics_HitSmear_HitKill_14272.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021MCRecon_KF_WithSpacing_pass6dev_lowp_physics_HitSmear_HitKill_14272.lcsim new file mode 100644 index 000000000..e96d2cb73 --- /dev/null +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021MCRecon_KF_WithSpacing_pass6dev_lowp_physics_HitSmear_HitKill_14272.lcsim @@ -0,0 +1,285 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 14 + false + + + false + + + pass5v9run14272_L1b_axial_hole.txt pass5v9run14272_L1b_stereo_hole.txt pass5v9run14272_L1t_axial_hole.txt pass5v9run14272_L1t_stereo_hole.txt pass5v9run14272_L2b_axial_hole.txt pass5v9run14272_L2b_stereo_hole.txt pass5v9run14272_L2t_axial_hole.txt pass5v9run14272_L2t_stereo_hole.txt + false + + + timeSmearing-2nsL4L7-4nsL1L3.txt + false + + + false + 0.8 + 0.8 + + + posSmearing_pass5_v9.txt + false + + + 1000 + + + HodoscopeReadoutHits + CONFIG + true + + + HodoscopeReadoutHits + 8 + CONFIG + true + + + + + + + 13.3 + + + WARNING + EcalClusters + + + EcalClusters + EcalClustersCorr + + + + SVTRawTrackerHits + + + .5 + 1 + Pileup + Migrad + true + 165 + true + true + false + true + false + false + true + false + + + 24.0 + 3.0 + false + 400 + 4.0 + 1.0 + 3.0 + 3.0 + true + true + false + + + true + true + 2 + 1 + + + 12.0 + + 55.0 + 20.0 + 11.546843987796496 + + + + + 9 + 9 + + 6 + 6 + + 3 + 3 + 40.0 + 5.0 + 12.320066328390354 + 9.206482863412027 + + + 8.0 + + 5.508828061070076 + 5 + false + + + 400 + + 0.3473319986601534 + true + 0.055 + 0.045 + true + -1.1 + + + + + 000BBS0 + 00BBS00 + 00ASBS0 + 00ABSS0 + 0A0SBS0 + 00B0BS0 + 00SBSA0 + 00SBB00 + 00SBAS0 + 0SA0BS0 + 0000SBB + 000SABS + 0BBS000 + 0SBB000 + 0000BBS + 000SSBA + 000BSSA + ABSS000 + SBB0000 + SABS000 + + BBS0000 + BSB0000 + 0BSB000 + + + 17 + + false + + + KalmanFullTracks + true + false + + + EcalClustersCorr + KalmanFullTracks + KalmanFullTracks + TrackClusterMatcherMinDistance + UnconstrainedV0Candidates_KF + UnconstrainedV0Vertices_KF + BeamspotConstrainedV0Candidates_KF + BeamspotConstrainedV0Vertices_KF + TargetConstrainedV0Candidates_KF + TargetConstrainedV0Vertices_KF + FinalStateParticles_KF + OtherElectrons_KF + true + false + true + false + 0.055 + 0.045 + -1.1 + 7.0 + 7.0 + 0.0 + 40.0 + 40 + 40 + false + true + true + false + true + true + UnconstrainedMollerCandidates_KF + UnconstrainedMollerVertices_KF + BeamspotConstrainedMollerCandidates_KF + BeamspotConstrainedMollerVertices_KF + TargetConstrainedMollerCandidates_KF + TargetConstrainedMollerVertices_KF + + + ${outputFile}_hit_eff.root + + + ${outputFile}.slcio + + + + + + diff --git a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass5_recon_skimmed_dataqual_physics.lcsim b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass5_recon_skimmed_dataqual_physics.lcsim index 464458884..e344df17c 100644 --- a/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass5_recon_skimmed_dataqual_physics.lcsim +++ b/steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021_pass5_recon_skimmed_dataqual_physics.lcsim @@ -187,7 +187,6 @@ ${outputFile}_v0skim.slcio FPGAData HelicalTrackHitRelations HelicalTrackHits HelicalTrackMCRelations KFGBLStripClusterData KFGBLStripClusterDataRelations ReadoutTimestamps RotatedHelicalTrackHitRelations RotatedHelicalTrackHits RotatedHelicalTrackMCRelations SVTFittedRawTrackerHits SVTShapeFitParameters SVTTrueHitRelations StripClusterer_SiTrackerHitStrip1D SVTRawTrackerHits FADCGenericHits HodoReadoutHits HodoCalHits EcalReadoutHits EcalUncalHits HodoGenericClusters VTPBank EcalClusters - all EcalClustersCorr diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanParams.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanParams.java index f11006b5c..45a7333c5 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanParams.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanParams.java @@ -31,7 +31,10 @@ public class KalmanParams { int mxShared; int [] minStereo; int minAxial; - double mxTdif; + double[] mxTdif; // Maximum hit time-range [ns], per iteration (index 0 = iter1, 1 = iter2) + double[] seedTimeSpread; // Max (tmax-tmin) [ns] among the 5 seed hits, per iteration (0 = use mxTdif) + double[] hitTimeWindow; // Max |t_hit - running mean| [ns] to pick up a hit in extension, per iteration (0 = legacy envelope) + double[] maxTotalSpread; // Backstop: max (tMax-tMin) [ns] over all hits on a track, per iteration (0 = disabled) double lowPhThresh; double seedCompThr; // Compatibility threshold for seedTracks helix parameters; ArrayList [] lyrList; @@ -39,6 +42,7 @@ public class KalmanParams { double [] vtxSize; double [] minSeedE; double edgeTolerance; + int debugEvent = -1; // event number to trace in KalmanPatRecHPS, -1 = none static final int numLayers = 14; private int[] Swap = {1,0, 3,2, 5,4, 7,6, 9,8, 11,10, 13,12}; @@ -78,7 +82,10 @@ public void print() { System.out.format(" Maximum residual, in units of detector resolution, for a hit to be shared: %8.2f\n", mxResidShare); System.out.format(" Maximum chi^2 increment to keep a shared hit: %8.2f\n", mxChi2double); System.out.format(" Maximum number of shared hits on a track: %d\n", mxShared); - System.out.format(" Maximum time difference among the hits on a track: %8.2f ns\n", mxTdif); + System.out.format(" Maximum time difference among the hits on a track (iter1, iter2): %8.2f, %8.2f ns\n", mxTdif[0], mxTdif[1]); + System.out.format(" Seed max time spread (iter1, iter2; 0=use mxTdif): %8.2f, %8.2f ns\n", seedTimeSpread[0], seedTimeSpread[1]); + System.out.format(" Hit time window vs running mean (iter1, iter2; 0=legacy): %8.2f, %8.2f ns\n", hitTimeWindow[0], hitTimeWindow[1]); + System.out.format(" Max total track time spread (iter1, iter2; 0=off): %8.2f, %8.2f ns\n", maxTotalSpread[0], maxTotalSpread[1]); System.out.format(" Threshold to remove redundant seeds (-1 to disable): %8.2f\n", seedCompThr); System.out.format(" Maximum chi^2 for 5-hit tracks with a vertex constraint: %8.2f\n", mxChi2Vtx); System.out.format(" Default origin to use for vertex constraints:\n"); @@ -117,6 +124,10 @@ public KalmanParams() { minHitsBot = new int[mxTrials]; mxResid = new double[mxTrials]; minStereo = new int[mxTrials]; + mxTdif = new double[mxTrials]; + seedTimeSpread = new double[mxTrials]; + hitTimeWindow = new double[mxTrials]; + maxTotalSpread = new double[mxTrials]; minSeedE = new double[numLayers]; for (int lyr=0; lyr= 0.0) kPar.setEdgeTolerance(edgeTolerance); + kPar.setDebugEvent(debugEvent); // if (minHits != 0) kPar.setMinHits(minHits); if (minHitsTopIter1 != 0) kPar.setMinHitsTopIter1(minHitsTopIter1); if (minHitsTopIter2 != 0) kPar.setMinHitsTopIter2(minHitsTopIter2); @@ -232,6 +248,22 @@ public void detectorChanged(Detector det) { if (maxSharedHits != 0) kPar.setMaxShared(maxSharedHits); if (maxTimeRange != 0.0) kPar.setMaxTimeRange(maxTimeRange); if (maxTanLambda != 0.0) kPar.setMaxTanL(maxTanLambda); + // Iteration-1 seed-cut overrides. These must come AFTER the tier-2 setters above, + // because setMaxK/setMaxdRho/setMaxdZ/setMaxTanL clamp the iteration-1 tier downward. + if (maxPtInverseIter1 != 0.0) kPar.setMaxKIter1(maxPtInverseIter1); + if (maxD0Iter1 != 0.0) kPar.setMaxdRhoIter1(maxD0Iter1); + if (maxZ0Iter1 != 0.0) kPar.setMaxdZIter1(maxZ0Iter1); + if (maxTanLambdaIter1 != 0.0) kPar.setMaxTanLIter1(maxTanLambdaIter1); + // Per-iteration hit time-range overrides. Must come AFTER setMaxTimeRange above, + // which sets BOTH iterations; these then override a single tier. + if (maxTimeRangeIter1 != 0.0) kPar.setMaxTimeRangeIter1(maxTimeRangeIter1); + if (maxTimeRangeIter2 != 0.0) kPar.setMaxTimeRangeIter2(maxTimeRangeIter2); + if (seedTimeSpreadIter1 != 0.0) kPar.setSeedTimeSpreadIter1(seedTimeSpreadIter1); + if (seedTimeSpreadIter2 != 0.0) kPar.setSeedTimeSpreadIter2(seedTimeSpreadIter2); + if (hitTimeWindowIter1 != 0.0) kPar.setHitTimeWindowIter1(hitTimeWindowIter1); + if (hitTimeWindowIter2 != 0.0) kPar.setHitTimeWindowIter2(hitTimeWindowIter2); + if (maxTotalSpreadIter1 != 0.0) kPar.setMaxTotalSpreadIter1(maxTotalSpreadIter1); + if (maxTotalSpreadIter2 != 0.0) kPar.setMaxTotalSpreadIter2(maxTotalSpreadIter2); if (maxResidual != 0.0) kPar.setMxResid(maxResidual); if (maxChi2Inc != 0.0) kPar.setMxChi2Inc(maxChi2Inc); if (minChi2IncBad != 0.0) kPar.setMinChi2IncBad(minChi2IncBad); @@ -633,6 +665,18 @@ public void setMaxZ0(double maxZ0) { public void setMaxChi2(double maxChi2) { this.maxChi2 = maxChi2; } + /** + * Tolerance in mm on the seed-acceptance check that the extrapolated helix lands + * inside a seed layer's active area. Set very large to disable the check and let + * chi2 decide instead. Default 1 mm, as before. + */ + public void setEdgeTolerance(double edgeTolerance) { + this.edgeTolerance = edgeTolerance; + } + /** Event number to trace seed-by-seed in KalmanPatRecHPS; -1 (default) traces none. */ + public void setDebugEvent(int debugEvent) { + this.debugEvent = debugEvent; + } // public void setMinHits(int minHits) { // this.minHits = minHits; // } @@ -659,9 +703,33 @@ public void setMaxSharedHits(int maxSharedHits) { public void setMaxTimeRange(double maxTimeRange) { this.maxTimeRange = maxTimeRange; } + public void setMaxTimeRangeIter1(double maxTimeRangeIter1) { + this.maxTimeRangeIter1 = maxTimeRangeIter1; + } + public void setMaxTimeRangeIter2(double maxTimeRangeIter2) { + this.maxTimeRangeIter2 = maxTimeRangeIter2; + } + public void setSeedTimeSpreadIter1(double v) { this.seedTimeSpreadIter1 = v; } + public void setSeedTimeSpreadIter2(double v) { this.seedTimeSpreadIter2 = v; } + public void setHitTimeWindowIter1(double v) { this.hitTimeWindowIter1 = v; } + public void setHitTimeWindowIter2(double v) { this.hitTimeWindowIter2 = v; } + public void setMaxTotalSpreadIter1(double v) { this.maxTotalSpreadIter1 = v; } + public void setMaxTotalSpreadIter2(double v) { this.maxTotalSpreadIter2 = v; } public void setMaxTanLambda(double maxTanLambda) { this.maxTanLambda = maxTanLambda; } + public void setMaxPtInverseIter1(double maxPtInverseIter1) { + this.maxPtInverseIter1 = maxPtInverseIter1; + } + public void setMaxD0Iter1(double maxD0Iter1) { + this.maxD0Iter1 = maxD0Iter1; + } + public void setMaxZ0Iter1(double maxZ0Iter1) { + this.maxZ0Iter1 = maxZ0Iter1; + } + public void setMaxTanLambdaIter1(double maxTanLambdaIter1) { + this.maxTanLambdaIter1 = maxTanLambdaIter1; + } public void setMaxResidual(double maxResidual) { this.maxResidual = maxResidual; } diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java index 5210cd3a4..7b7f9bf2d 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java @@ -46,7 +46,13 @@ class KalmanPatRecHPS { private ArrayList XLscat; private int eventNumber; - private static final boolean debug = false; + /** + * Seed-by-seed tracing. Was a compile-time constant, so none of it could be reached + * without editing and rebuilding; it is now armed per event from + * KalmanParams.debugEvent (default -1, i.e. off), which restores the intent of the + * commented-out line below without costing anything when it is not in use. + */ + private boolean debug = false; private int nModules; private KalmanParams kPar; private Logger logger; @@ -101,9 +107,9 @@ void patRecSetVtx(double [] vtx, double [][] vtxCov) { ArrayList kalmanPatRec(EventHeader event, Map hitMapHPS, ArrayList data, int topBottom) { // topBottom = 0 for the bottom tracker (z>0); 1 for the top tracker (z<0) - if (event != null) eventNumber = event.getEventNumber(); + if (event != null) eventNumber = event.getEventNumber(); else eventNumber++; - //debug = (eventNumber == 16276); + debug = (kPar.debugEvent != -1 && eventNumber == kPar.debugEvent); if (debug) startTime = System.nanoTime(); int nCandHits = 0; @@ -292,7 +298,8 @@ ArrayList kalmanPatRec(EventHeader event, Map tmin = Math.min(tmin, ht.hit.time); tmax = Math.max(tmax, ht.hit.time); } - if (tmax - tmin > kPar.mxTdif) { + double seedSpreadMax = (kPar.seedTimeSpread[trial] > 0.) ? kPar.seedTimeSpread[trial] : kPar.mxTdif[trial]; + if (tmax - tmin > seedSpreadMax) { if (debug) { System.out.format("KalmanPatRecHPS: skipping seed %d %d %d %d %d with tdif=%8.2f\n Hits: ", idx[0], idx[1], idx[2], idx[3], idx[4], tmax-tmin); @@ -692,7 +699,7 @@ ArrayList kalmanPatRec(EventHeader event, Map // Check if the track can be improved by removing hits if (removeBadHits(candidateTrack, minHits1, trial)) { if (debug) System.out.format("KalmanPatRecHPS: Refit candidate track %d after removing a hit.\n", candidateTrack.ID); - if (candidateTrack.reFit()) { + if (candidateTrack.reFit(trial)) { if (debug) candidateTrack.print("after refitting and smoothing", false); } else { candidateTrack.good = false; @@ -912,7 +919,7 @@ ArrayList kalmanPatRec(EventHeader event, Map if (!MatrixFeatures_DDRM.hasNaN(aS0.helix.C)) { Collections.sort(tkr.sites, MeasurementSite.SiteComparatorUp); // Occasionally necessary if (tkr.sites.get(0).aS == null) { - if (tkr.reFit()) { + if (tkr.reFit(trial)) { tkr.good = true; if (debug) { System.out.format("KalmanPatRecHPS event %d: resurrecting refit candidate %d with chi2=%9.5f\n", @@ -1138,7 +1145,7 @@ ArrayList kalmanPatRec(EventHeader event, Map } if (site.chi2inc > kPar.mxChi2double) { if (!site.smoothed) logger.log(Level.WARNING,String.format("OOPS, why isn't this site smoothed at layer %d?",site.m.Layer)); - if (tkr.removeHit(site, kPar.mxChi2Inc, kPar.mxTdif)) { + if (tkr.removeHit(site, kPar.mxChi2Inc, Math.max(kPar.mxTdif[0], kPar.mxTdif[1]))) { if (debug) { System.out.format("KalmanPatRecHPS: added a hit after removing one for Track %d, Layer %d\n",tkr.ID, module.Layer); } @@ -1167,7 +1174,7 @@ ArrayList kalmanPatRec(EventHeader event, Map } // Try to add hits on layers with missing hits - int nAdded = tkr.addHits(data, kPar.mxResid[1], kPar.mxChi2Inc, kPar.mxTdif, debug); + int nAdded = tkr.addHits(data, kPar.mxResid[1], kPar.mxChi2Inc, Math.max(kPar.mxTdif[0], kPar.mxTdif[1]), debug); // check that there are enough hits in both views int nStereo = 0; @@ -1538,7 +1545,17 @@ private void filterTrack(TrackCandidate tkrCandidate, int lyrBegin, // layer on } newSite = new MeasurementSite(lyr, m, kPar); int rF; - double [] tRange = {tkrCandidate.tMax - kPar.mxTdif, tkrCandidate.tMin + kPar.mxTdif}; + double [] tRange; + if (kPar.hitTimeWindow[trial] > 0.) { // window vs running mean of accepted hits, with total-spread backstop + double tWin = kPar.hitTimeWindow[trial]; + double tSpreadMax = (kPar.maxTotalSpread[trial] > 0.) ? kPar.maxTotalSpread[trial] : 1.e10; + double tSum = 0.; int nT = 0; + for (KalHit htm : tkrCandidate.hits) { tSum += htm.hit.time; nT++; } + double tMean = (nT > 0) ? tSum/nT : 0.5*(tkrCandidate.tMin + tkrCandidate.tMax); + tRange = new double[]{Math.max(tMean - tWin, tkrCandidate.tMax - tSpreadMax), Math.min(tMean + tWin, tkrCandidate.tMin + tSpreadMax)}; + } else { + tRange = new double[]{tkrCandidate.tMax - kPar.mxTdif[trial], tkrCandidate.tMin + kPar.mxTdif[trial]}; + } if (prevSite == null) { // For first layer use the initializer state vector boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made on last module of the layer rF = newSite.makePrediction(sI, null, hitno, tkrCandidate.nTaken <= kPar.mxShared, pickUp, checkBounds, tRange, trial); diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackCandidate.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackCandidate.java index 61d47fe10..b30fc673c 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackCandidate.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackCandidate.java @@ -218,7 +218,7 @@ void removeHit(KalHit hit, boolean deleteFromList) { if (nstr < 3 || nax < 2) good = false; } - boolean reFit() { + boolean reFit(int trial) { // trial only selects the per-iteration time window; residual cuts unchanged final boolean verbose = false; if (verbose) System.out.format("TrackCandidate.reFit: starting filtering for event %d.\n",eventNumber); @@ -262,7 +262,17 @@ boolean reFit() { boolean allowSharing = nTaken < kPar.mxShared; boolean checkBounds = false; - double [] tRange = {tMax - kPar.mxTdif, tMin + kPar.mxTdif}; + double [] tRange; + if (kPar.hitTimeWindow[trial] > 0.) { + double tWin = kPar.hitTimeWindow[trial]; + double tSpreadMax = (kPar.maxTotalSpread[trial] > 0.) ? kPar.maxTotalSpread[trial] : 1.e10; + double tSum = 0.; int nT = 0; + for (KalHit htm : hits) { tSum += htm.hit.time; nT++; } + double tMean = (nT > 0) ? tSum/nT : 0.5*(tMin + tMax); + tRange = new double[]{Math.max(tMean - tWin, tMax - tSpreadMax), Math.min(tMean + tWin, tMin + tSpreadMax)}; + } else { + tRange = new double[]{tMax - kPar.mxTdif[trial], tMin + kPar.mxTdif[trial]}; + } int rF = currentSite.makePrediction(sH, prevMod, currentSite.hitID, allowSharing, pickupHits, checkBounds, tRange, 0); if (rF < 0) { if (verbose) System.out.format("TrackCandidate.reFit: failed to make prediction at layer %d for event %d!\n",currentSite.m.Layer,eventNumber);