From 44ee66611145b33a029c6e599381262471d6f5a6 Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:37:20 -0700 Subject: [PATCH 1/7] Kalman: per-iteration hit timing cuts and iteration-1 seed overrides Make mxTdif a per-iteration array and add three new per-iteration timing knobs alongside it: seedTimeSpread gate on tmax-tmin across the 5 seed hits hitTimeWindow pick up hits within a window of the running mean of accepted hit times, replacing the tMin/tMax envelope maxTotalSpread backstop on total tMax-tMin over the track All three default to 0, which selects the existing behaviour, and setMaxTimeRange() still sets both iterations, so current steering files are unaffected. removeHit()/addHits() take max(mxTdif[0], mxTdif[1]) since they have no iteration context. TrackCandidate.reFit() gains a trial argument to select the time window. Also expose iteration-1 seed cuts (setMaxKIter1, setMaxdRhoIter1, setMaxdZIter1, setMaxTanLIter1), which must be applied after the tier-2 setters because those clamp the iteration-1 tier downward; and edgeTolerance, the seed detector-bounds check that acts as an acceptance cut at low tan(lambda), previously hard-coded at 1 mm. KalmanPatRecHPS.debug is now armed per event from KalmanParams.debugEvent (default -1, off) instead of being a compile-time constant. --- .../recon/tracking/kalman/KalmanParams.java | 125 +++++++++++++++++- .../tracking/kalman/KalmanPatRecDriver.java | 70 +++++++++- .../tracking/kalman/KalmanPatRecHPS.java | 35 +++-- .../recon/tracking/kalman/TrackCandidate.java | 14 +- 4 files changed, 225 insertions(+), 19 deletions(-) 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 f11006b5cb..45a7333c55 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 5210cd3a4c..7b7f9bf2d1 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 61d47fe106..b30fc673c5 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); From 4969d1cbf7cc8479f31b7b8caab1c8ee968bde6a Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:41:35 -0700 Subject: [PATCH 2/7] SVT digi: truth-relation diagnostics and optional hit provenance Under pulser overlay many raw hits reach readout with no truth relation. Two independent gates can drop one: G1 (totalContrib > 4*meanNoise), which attaches simHits to the channel, and G2 (readoutCuts), evaluated on the combined waveform whose baseline need not match the DB pedestal. Add counters for both, indexed by whether the channel carried a pulser hit, plus a pulser-baseline-vs-pedestal measurement, reported at endOfData(). These are always accumulated and change no behaviour; setDebug() additionally prints per-hit loss records. Add opt-in writeHitOriginCollections (default false), writing SVTHitOriginPulser and SVTHitOriginMCContrib. Combined with SVTTrueHitRelations these resolve "no truth relation" into NOISE, MC_PURE, MC_PURE_SUBTHRESH, PULSER_PURE, MERGED or MERGED_SUBTHRESH, separating a lost relation from an absent one. Enabling it widens the output MCParticle collection to cover sub-threshold contributors, so output is only bit-identical with it disabled. --- ...ionWithPulserDataMergingReadoutDriver.java | 370 +++++++++++++++++- 1 file changed, 363 insertions(+), 7 deletions(-) diff --git a/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java b/digi/src/main/java/org/hps/digi/SvtDigitizationWithPulserDataMergingReadoutDriver.java index 41a69282e3..cef66c006c 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. From 409c0815ae89846b414807e96b33254ae585796e Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:43:42 -0700 Subject: [PATCH 3/7] Add SVT hit provenance and event forensics drivers Two analysis drivers consuming the provenance relations written by SvtDigitizationWithPulserDataMergingReadoutDriver; both require writeHitOriginCollections=true. SvtHitProvenanceDriver answers the aggregate question: are the untruthed hits on a track adjacent to channels that did receive MC charge? Since high occupancy makes adjacency alone meaningless, it is measured against the baseline rate for any pulser hit in the event. SvtEventForensicsDriver dumps named events hit by hit, separating a track whose particle never produced a truthed raw hit from one where truthed raw hits existed but pattern recognition used different ones. --- .../analysis/MC/SvtEventForensicsDriver.java | 397 ++++++++++++++++++ .../analysis/MC/SvtHitProvenanceDriver.java | 297 +++++++++++++ 2 files changed, 694 insertions(+) create mode 100644 analysis/src/main/java/org/hps/analysis/MC/SvtEventForensicsDriver.java create mode 100644 analysis/src/main/java/org/hps/analysis/MC/SvtHitProvenanceDriver.java 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 0000000000..aa03472ccf --- /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 0000000000..feccedad1c --- /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(); + } +} From 3568f7780891f92ce4354c2072661c1d5f984bab Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:44:11 -0700 Subject: [PATCH 4/7] UnbiasedTriggerFilterDriver: add counters and accept/reject summary Count events by the reason they were kept or dropped (EPICS, scaler, pulser, FaradayCup, missing TSBank, other) and print a summary at endOfData(). Optional per-event logging via debug=true in the steering file. Pulser and FaradayCup were tested in a single condition and are now separate branches so they can be counted apart; the accept/reject outcome is unchanged. --- .../UnbiasedTriggerFilterDriver.java | 77 +++++++++++++++++-- 1 file changed, 72 insertions(+), 5 deletions(-) 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 2511ec6817..1b95fcf9d8 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("========================================================="); + } } From 6ed19f7678feca6d17c01bacfa7f68639f699513 Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:44:48 -0700 Subject: [PATCH 5/7] Remove stray closing tag in pass5 physics recon steering The extra tag followed the v0 skim driver block and did not match any open element. --- .../PhysicsRun2021_pass5_recon_skimmed_dataqual_physics.lcsim | 1 - 1 file changed, 1 deletion(-) 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 464458884c..e344df17c7 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 From 18827f60f2a57af108be8740580698ea6c75467d Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:44:50 -0700 Subject: [PATCH 6/7] Regenerate detector-model dependency-reduced-pom (junit 4.13.1 -> 4.13.2) Shade plugin output picking up the junit version already declared in the parent pom; no source change. --- detector-model/dependency-reduced-pom.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/detector-model/dependency-reduced-pom.xml b/detector-model/dependency-reduced-pom.xml index a57ace5df4..6be2e55332 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 From b9ee85fe526ccb2366c4035115a3a7609228a0a6 Mon Sep 17 00:00:00 2001 From: Matthew Gignac Date: Thu, 6 Aug 2026 10:52:43 -0700 Subject: [PATCH 7/7] adding lowp MC steering file --- ..._lowp_physics_HitSmear_HitKill_14272.lcsim | 285 ++++++++++++++++++ 1 file changed, 285 insertions(+) create mode 100644 steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2021MCRecon_KF_WithSpacing_pass6dev_lowp_physics_HitSmear_HitKill_14272.lcsim 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 0000000000..e96d2cb736 --- /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 + + + + + +