From cdeb48301b54182f4c28ad6ba1e540225364b8b0 Mon Sep 17 00:00:00 2001 From: robertprestonjohnson Date: Thu, 6 May 2021 13:25:40 -0700 Subject: [PATCH 1/3] Fix some problems found by Norman when looking for tracks flagged by ECAL clusters. All some hit sharing by candidate tracks that get raised from the dead. Fixed bugs in early cuts made on seed formation, which were meant to speed up the code, and set those cuts better, based on distributions obtained from good tracks. More debug printout and histograms, and small adjustments to some default cuts, including addition of two seed strategies. --- .../org/hps/recon/tracking/kalman/KalHit.java | 4 +- .../recon/tracking/kalman/KalmanParams.java | 23 ++- .../tracking/kalman/KalmanPatRecHPS.java | 152 ++++++++++++------ .../tracking/kalman/KalmanPatRecPlots.java | 32 ++++ .../recon/tracking/kalman/TrackCandidate.java | 4 +- 5 files changed, 148 insertions(+), 67 deletions(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalHit.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalHit.java index b74db082d9..f81b5e2416 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalHit.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalHit.java @@ -33,9 +33,9 @@ String toString(String s) { if (s=="short") { int idx = module.hits.indexOf(hit); if (module.isStereo) { - str = String.format(" {%d %d %d %d %5.1f} ", module.Layer, module.detector, idx, ntks, hit.time); + str = String.format(" {%d %d %d %d %5.1f %6.3f} ", module.Layer, module.detector, idx, ntks, hit.time, hit.energy); } else { - str = String.format(" (%d %d %d %d %5.1f) ", module.Layer, module.detector, idx, ntks, hit.time); + str = String.format(" (%d %d %d %d %5.1f %6.3f) ", module.Layer, module.detector, idx, ntks, hit.time, hit.energy); } } else { str = String.format("Hit %s in layer %d, detector %d, hit %d, value=%10.5f, #tkrs=%d, candidate chi2=", s, module.Layer, module.detector, module.hits.indexOf(hit), hit.v, ntks); 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 e637d9dd41..505d25490e 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 @@ -3,14 +3,9 @@ import java.util.ArrayList; import java.util.logging.Level; import java.util.logging.Logger; -//import org.hps.conditions.beam.BeamPosition; -//import org.hps.conditions.beam.BeamPosition.BeamPositionCollection; -//import org.hps.conditions.database.DatabaseConditionsManager; /** * Parameters used by the Kalman-Filter pattern recognition and fitting - * @author Robert Johnson - * */ public class KalmanParams { static final int mxTrials = 2; // Max number of iterations through the entire pattern recognition; not configurable @@ -119,7 +114,7 @@ public void print() { minSeedE = new double[numLayers]; for (int lyr=0; lyr XLscat; private int eventNumber; - private final static boolean debug = false; + private static boolean debug = false; private int nModules; private KalmanParams kPar; private Logger logger; @@ -99,7 +97,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev // topBottom = 0 for the bottom tracker (z>0); 1 for the top tracker (z<0) this.eventNumber = eventNumber; - //debug = (eventNumber == 176782359); + debug = (eventNumber == 8882492); if (debug) startTime = System.nanoTime(); @@ -152,17 +150,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev } } } - if (debug) { - for (int lyr = 0; lyr < lyrHits.size(); ++lyr) { - ArrayList LL = lyrHits.get(lyr); - System.out.format("KalmanPatRecHPS: layer %d hits:", lyr); - for (KalHit ht : LL) { - ht.print("short"); - } - System.out.format("\n"); - } - } - + if (debug) { System.out.format(" KalmanPatRecHPS: list of the seed strategies to be applied:\n"); for (int[] list : kPar.lyrList[topBottom]) { @@ -187,6 +175,14 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev if (debug) { double runTime = (double)((System.nanoTime() - startTime)/1000000.); System.out.format("\nKalmanPatRecHPS: start of pass %d through the algorithm. Time %10.6f ms\n", trial, runTime); + for (int lyr = 0; lyr < lyrHits.size(); ++lyr) { + ArrayList LL = lyrHits.get(lyr); + System.out.format("KalmanPatRecHPS: layer %d hits:", lyr); + for (KalHit ht : LL) { + ht.print("short"); + } + System.out.format("\n"); + } } // Remove references from the hits to the track candidates of the previous iteration @@ -285,11 +281,12 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev tmax = Math.max(tmax, ht.hit.time); } if (tmax - tmin > kPar.mxTdif) { - //if (debug) { - // System.out.format("KalmanPatRecHPS: skipping seed with tdif=%8.2f\n Hits: ", tmax-tmin); - // for (KalHit ht : hitList) ht.print("short"); - // System.out.format("\n"); - //} + 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); + for (KalHit ht : hitList) ht.print("short"); + System.out.format("\n"); + } continue; } // To avoid wasting time fitting seeds, skip seeds that are entirely contained in already found candidates @@ -306,7 +303,12 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev // Fit the seed to extract helix parameters SeedTrack seed = new SeedTrack(hitList, yOrigin, kPar.beamSpot[1]); - if (!seed.success) continue; + if (!seed.success) { + if (debug) { + System.out.format("Seed %d %d %d %d %d failed fit\n",idx[0], idx[1], idx[2], idx[3], idx[4]); + } + continue; + } // Cuts on the seed quality Vec hp = seed.helixParams(); @@ -430,7 +432,8 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev // Go back to the seed hits on the filtered track and check the detector bounds for (MeasurementSite site : candidateTrack.sites) { - if (debug) System.out.format("KalmamPatRecHPS: check seed hit layer boundaries for candidate %d\n", candidateTrack.ID); + if (debug) System.out.format("KalmamPatRecHPS: check seed hit layer %d boundaries for candidate %d\n", + site.m.Layer, candidateTrack.ID); int nseedht = 0; for (int i=0; i kalmanPatRec(ArrayList data, int topBottom, int ev } if (!seedlyr) continue; HelixState hx = site.aP.helix; - double phi = 0.; //hx.planeIntersect(site.m.p); The helix at this point starts from the intersection, so phi is always 0 here. - if (debug) System.out.format(" phi at layer %d is %10.6f; measurement=%10.5f\n", site.m.Layer, phi, site.m.hits.get(site.hitID).v); + double phi = 0.; // The helix at this point starts from the intersection, so phi is always 0 here. + if (debug) { + System.out.format(" phi at layer %d is %10.6f; measurement=%10.5f\n", site.m.Layer, phi, site.m.hits.get(site.hitID).v); + double phiChk = hx.planeIntersect(site.m.p); + System.out.format(" phi calculated = %10.6f\n", phiChk); + } //if (Double.isNaN(phi)) continue seedLoop; Vec rGlob = hx.toGlobal(hx.atPhi(phi)); Vec rDet = site.m.toLocal(rGlob); if (debug) { System.out.format(" global = %s, local = %s\n", rGlob.toString(), rDet.toString()); - System.out.format(" extents: x=%8.3f->%8.3f y=%8.3f->%8.3f\n", - site.m.xExtent[0], site.m.xExtent[1], site.m.yExtent[0], site.m.yExtent[1]); + System.out.format(" extents: x=%8.3f->%8.3f y=%8.3f->%8.3f, tolerance=%8.3f\n", + site.m.xExtent[0], site.m.xExtent[1], site.m.yExtent[0], site.m.yExtent[1],kPar.edgeTolerance); + } + if (rDet.v[0] > site.m.xExtent[1] + kPar.edgeTolerance) { + if (debug) { + System.out.format(" The x coordinate %8.3f is too large; reject!\n", rDet.v[0]); + hx.print(" the offender"); + seed.print(" the originator"); + } + continue seedLoop; + } + if (rDet.v[0] < site.m.xExtent[0] - kPar.edgeTolerance) { + if (debug) System.out.format(" The x coordinate %8.3f is too small; reject!\n", rDet.v[0]); + continue seedLoop; + } + if (rDet.v[1] > site.m.yExtent[1] + kPar.edgeTolerance) { + if (debug) System.out.format(" The y coordinate %8.3f is too large; reject!\n", rDet.v[1]); + continue seedLoop; + } + if (rDet.v[1] < site.m.yExtent[0] - kPar.edgeTolerance) { + if (debug) System.out.format(" The y coordinate %8.3f is too small; reject!\n", rDet.v[0]); + continue seedLoop; } - if (rDet.v[0] > site.m.xExtent[1] + kPar.edgeTolerance) continue seedLoop; - if (rDet.v[0] < site.m.xExtent[0] - kPar.edgeTolerance) continue seedLoop; - if (rDet.v[1] > site.m.yExtent[1] + kPar.edgeTolerance) continue seedLoop; - if (rDet.v[1] < site.m.yExtent[0] - kPar.edgeTolerance) continue seedLoop; if (site.m.split) { Measurement hit = site.m.hits.get(site.hitID); if (debug) hit.print("lyr 0,1"); @@ -564,7 +587,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev KalHit hitR = hitMap.get(siteR.m.hits.get(siteR.hitID)); if (hitR != null) { //if (!candidateTrack.hits.contains(hitR)) System.out.format("Oops, missing hit!\n"); - candidateTrack.removeHit(hitR); + candidateTrack.removeHit(hitR, true); if (debug) System.out.format("KalmanPatRecHPS event %d candidate %d, removing hit from layer %d detector %d\n", eventNumber, candidateTrack.ID, siteR.m.Layer, siteR.m.detector); removedHit = true; @@ -750,7 +773,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev } } for (TrackCandidate tkr : tksToRemove) { - tkr.removeHit(ht); // This will also remove the reference from the hit to this candidate + tkr.removeHit(ht, true); // This will also remove the reference from the hit to this candidate if (debug) tkr.print("after hit removal", true); } } @@ -803,22 +826,39 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev TrackCandidate tkr = iter.next(); if (!tkr.good) { // Resurrect the candidate if it has enough hits and none of them is shared with a good candidate or finished track - boolean resurrect = true; - for (KalHit ht : tkr.hits) { - if (ht.hit.tracks.size() > 0) { - resurrect = false; - break; - } - int nGood = 0; - for (TrackCandidate tkr2 : ht.tkrCandidates) { - if (tkr2.good) nGood++; + boolean resurrect = (trial == KalmanParams.mxTrials-1); + if (resurrect) { + if (debug) { + System.out.format("KalmanPatRecHPS: considering whether to resurrect candidate %d\n", tkr.ID); } - if (nGood > 0) { - resurrect = false; - break; + Iterator itrht = tkr.hits.iterator(); + while(itrht.hasNext()) { + int nShared = 0; + KalHit ht = itrht.next(); + if (ht.hit.tracks.size() > 0) { + if (debug) System.out.format(" a hit on layer %d was shared with track %d\n", ht.module.Layer, ht.hit.tracks.get(0).ID); + nShared++; + } + for (TrackCandidate tkr2 : ht.tkrCandidates) { + if (tkr2.good) { + if (debug) System.out.format(" a hit on layer %d was shared with track candidate %d\n", ht.module.Layer, tkr2.ID); + nShared++; + } + } + if (nShared > kPar.mxShared) { + if (debug) System.out.format("KalmanPatRecHPS: not resurrecting track candidate %d in event %d for excessive hit sharing.\n", + tkr.ID, eventNumber); + resurrect = false; + break; + } } } if (resurrect) { + if (debug) { + System.out.format(" # hits=%d vs %d, # stereo=%d vs %d, # axial=%d vs %d, helix=%s\n", + tkr.numHits(), kPar.minHits1[trial], tkr.numStereo(), kPar.minStereo[trial], tkr.numHits()-tkr.numStereo(), + kPar.minAxial, tkr.sites.get(0).aS.helix.a.toString()); + } if (tkr.numHits() >= kPar.minHits1[trial] && tkr.numStereo() >= kPar.minStereo[trial]) { int nAxial = tkr.numHits() - tkr.numStereo(); if (nAxial >= kPar.minAxial) { @@ -924,7 +964,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev System.out.format("KalmanPatRecHPS: hit %d removed from track candidate %d on layer %d detector %d\n", hit.module.hits.indexOf(hit.hit), tkr.ID, hit.module.Layer, hit.module.detector); } - tkr.removeHit(hit); + tkr.removeHit(hit, true); } } } @@ -1312,7 +1352,7 @@ private boolean removeBadHits(TrackCandidate tkr, int trial) { KalHit badHit = hitMap.get(badOne); if (badHit != null) { if (debug) System.out.format("KalmanPatRecHPS.removeBadHits: event %d, removing bad hit in layer %d with chi2=%9.3f.\n",eventNumber, tkr.sites.get(idxBad).m.Layer,mxChi2); - tkr.removeHit(badHit); + tkr.removeHit(badHit, true); return true; } } @@ -1677,10 +1717,22 @@ private boolean seedNoGood(int j, int iter) { for (int i=0; i 1.1*kPar.dRhoMax[iter]) return true; - if (Math.abs(slope) > 1.1*kPar.tanlMax[iter]) return true; + if (Math.abs(zIntercept) > kPar.dzMax[iter]) { + if (debug) { + System.out.format("KalmanPatRecHPS.seedNoGood: reject z j=%d, i=%d, zInt=%10.5f vs %10.5f, slope=%10.5f vs %10.5f\n", + j, i, zIntercept, kPar.dzMax[iter], slope, kPar.tanlMax[iter]); + } + return true; + } + if (Math.abs(slope) > kPar.tanlMax[iter]) { + if (debug) { + System.out.format("KalmanPatRecHPS.seedNoGood: reject slope j=%d, i=%d, zInt=%10.5f vs %10.5f, slope=%10.5f vs %10.5f\n", + j, i, zIntercept, kPar.dzMax[iter], slope, kPar.tanlMax[iter]); + } + return true; + } return false; } } diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java index a7e4e68a56..9817fda547 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecPlots.java @@ -189,6 +189,9 @@ class KalmanPatRecPlots { aida.histogram1D("Bad/momentum of bad tracks", 60, 0., 6.); aida.histogram1D("Bad/Number of MC particles associated", 10, 0., 10.); aida.histogram1D("Bad/Number of wrong hits on track", 20, 0., 20.); + aida.histogram1D("seed slope",100,0.,0.3); + aida.histogram1D("seed z intercept",100,0.,10.); + aida.histogram1D("seed y intercept",100,-100.,100.); } void process(EventHeader event, List sensors, ArrayList[] kPatList, @@ -422,6 +425,33 @@ void process(EventHeader event, List sensors, ArrayList[] aida.histogram1D("Kalman track drho").fill(kTk.originHelixParms()[0]); aida.histogram1D("Kalman track dz").fill(kTk.originHelixParms()[3]); aida.histogram1D("Kalman track time range (ns)").fill(kTk.tMax - kTk.tMin); + + // Use good tracks to analyze seed cuts + if (kTk.nHits >=10) { + if (kTk.chi2 < 15.) { + for (MeasurementSite site : kTk.SiteList) { + if (site.m.isStereo) continue; + if (site.hitID < 0) continue; + Measurement hit1 = site.m.hits.get(site.hitID); + double z1 = site.m.toGlobal(new Vec(0.,hit1.v,0.)).v[2]; + double y1 = site.m.p.X().v[1]; + for (MeasurementSite site2 : kTk.SiteList) { + if (site2.m.isStereo) continue; + if (site.m.Layer == site2.m.Layer) continue; + if (site2.hitID < 0) continue; + Measurement hit2 = site2.m.hits.get(site2.hitID); + double z2 = site2.m.toGlobal(new Vec(0.,hit2.v,0.)).v[2]; + double y2 = site2.m.p.X().v[1]; + double slope = (z2 - z1) / (y2 - y1); + double zIntercept = z1 - slope * y1; + double yIntercept = -(zIntercept/slope); + aida.histogram1D("seed slope").fill(Math.abs(slope)); + aida.histogram1D("seed z intercept").fill(Math.abs(zIntercept)); + aida.histogram1D("seed y intercept").fill(yIntercept); + } + } + } + } // Check the covariance matrix boolean badCov = false; @@ -846,6 +876,8 @@ void process(EventHeader event, List sensors, ArrayList[] //int [] badEvents = {51753, 52531, 56183, 57958, 58050, 60199, 80324, 83798, 84933, 86351, 88796, 96749, 97230, 102986, 105578, 106654, // 107191, 108542, 108886, 110453, 120457, 121129, 121311, 121525, 124355, 124910, 127335, 129360, 133951}; + //int [] badEvents = {8788317, 8730045, 8724483, 8716465, 8697779, 8696553, 8565879, 8566151, 8563109, 196192, 155656, 144460, 135200, + // 110480, 69464, 20665, 9000031, 8882086, 8878414, 8750999}; int [] badEvents = {}; if (nPlotted < numEvtPlots) { // && sharedHitTrack) { boolean plotIt = false; 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 fe2566af8d..d809a3094f 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 @@ -83,7 +83,7 @@ Vec originHelix() { return site0.aS.helix.pivotTransform(); } - void removeHit(KalHit hit) { + void removeHit(KalHit hit, boolean deleteFromList) { MeasurementSite siteR = null; SiModule mod = hit.module; for (MeasurementSite site : sites) { @@ -112,7 +112,7 @@ void removeHit(KalHit hit) { tMax = Math.min(tMax, ht.hit.time); } } - hits.remove(hit); + if (deleteFromList) hits.remove(hit); hit.tkrCandidates.remove(this); int nstr = numStereo(); int nax = numHits() - nstr; From f2ce3ea27a733f781fcf98fad8c264d67f47a87f Mon Sep 17 00:00:00 2001 From: robertprestonjohnson Date: Thu, 6 May 2021 13:30:58 -0700 Subject: [PATCH 2/3] Forgot in the previous commit to remove a debug printout. That is fixed here. --- .../java/org/hps/recon/tracking/kalman/KalmanPatRecHPS.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 8fe01367a3..9b52a0ae9d 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 @@ -44,7 +44,7 @@ class KalmanPatRecHPS { private ArrayList XLscat; private int eventNumber; - private static boolean debug = false; + private static final boolean debug = false; private int nModules; private KalmanParams kPar; private Logger logger; @@ -97,7 +97,7 @@ ArrayList kalmanPatRec(ArrayList data, int topBottom, int ev // topBottom = 0 for the bottom tracker (z>0); 1 for the top tracker (z<0) this.eventNumber = eventNumber; - debug = (eventNumber == 8882492); + //debug = (eventNumber == 8882492); if (debug) startTime = System.nanoTime(); From 14d891f887eb52eb92d58a7cdb4fc724f7d34a34 Mon Sep 17 00:00:00 2001 From: robertprestonjohnson Date: Sat, 8 May 2021 16:16:38 -0700 Subject: [PATCH 3/3] Turned off a debug printout --- .../main/java/org/hps/recon/tracking/kalman/StateVector.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/StateVector.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/StateVector.java index 609446d267..0991916f55 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/StateVector.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/StateVector.java @@ -384,7 +384,7 @@ StateVector smooth(StateVector snS, StateVector snP) { vecToM(snS.helix.a.dif(snP.helix.a), tempV); CommonOps_DDRM.mult(tempA, tempV, tempV2); sS.helix.a = helix.a.sum(mToVec(tempV2)); - if (debug || Math.abs(tempV2.unsafe_get(1, 0))>1.5) { + if (debug) { System.out.println("StateVector:smooth, inverse of the covariance:"); Cinv.print("%11.6e"); CommonOps_DDRM.mult(snP.helix.C, Cinv, tempM);