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 final 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/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); 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;