diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java index 1168bc8672..14d07ca340 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixState.java @@ -93,8 +93,8 @@ String toString(String s) { String str; str = String.format("HelixState %s: helix parameters=%s, pivot=%s\n", s, a.toString(), X0.toString()); str = str + String.format(" Origin=%s, B=%10.6f in direction %s\n", origin.toString(), B, tB.toString()); - str = str + "Covariance:" + C.toString(); - str = str + Rot.toString("from global coordinates to field coordinates"); + if (C != null) str = str + "Covariance:" + C.toString(); + if (Rot != null) str = str + Rot.toString("from global coordinates to field coordinates"); str = str + "End of HelixState dump\n"; return str; } diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java index 3cc6f21c2b..11fd010444 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixTest3.java @@ -38,7 +38,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code // Control parameters // Units are Tesla, GeV, mm - int nTrials = 10000; // The number of test events to generate for fitting + int nTrials = 100; // The number of test events to generate for fitting int startLayer = 10; // Where to start the Kalman filtering int nIteration = 2; // Number of filter iterations int nAxial = 3; // Number of axial layers needed by the linear fit @@ -49,6 +49,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code boolean verbose = nTrials < 2; double executionTime = 0.; + int nBadCov = 0; double eCalLoc = 1394.; @@ -97,9 +98,9 @@ class HelixTest3 { // Program for testing the Kalman fitting code Vec B = new Vec(3, fM.getField(new Vec(0., y, z))); System.out.format("x=0 y=%6.1f z=%6.1f: %s\n", y, z, B.toString()); } - System.out.format("B field map vs x at ECAL:\n"); - for (double x=-200.; x<200.; x+=5.) { - double y=eCalLoc + 10.; + System.out.format("B field map vs x at center:\n"); + for (double x=-300.; x<300.; x+=5.) { + double y=505.; double z=20.; Vec B = new Vec(3, fM.getField(new Vec(x, y, z))); System.out.format("x=%6.1f y=%6.1f z=%6.1f: %s\n", x, y, z, B.toString()); @@ -418,7 +419,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code hResidS4[i] = new Histogram(100, -0.1, 0.002, String.format("Smoothed true residual for plane %d", i), "mm", "hits"); hResidX[i] = new Histogram(100, -0.8, 0.016, String.format("True residual in global X for plane %d", i), "mm", "hits"); hResidZ[i] = new Histogram(100, -0.1, 0.002, String.format("True residual in global Z for plane %d", i), "mm", "hits"); - hUnbias[i] = new Histogram(100, -0.2, 0.004, String.format("Unbiased residual for plant %d", i), "mm", "hits"); + hUnbias[i] = new Histogram(100, -0.2, 0.004, String.format("Unbiased residual for plane %d", i), "mm", "hits"); hUnbiasSig[i] = new Histogram(100, -10., 0.2, String.format("Unbiased residuals for layer %d", i), "sigmas", "hits"); } Histogram hpropx = new Histogram(100,-5.,0.1,"projected track-state x error","mm","track"); @@ -740,11 +741,14 @@ class HelixTest3 { // Program for testing the Kalman fitting code Matrix C = new Matrix(KalmanTrack.originCovariance()); EigenvalueDecomposition eED= new EigenvalueDecomposition(C); double [] ev = eED.getRealEigenvalues(); + boolean badCov = false; for (int i=0; i<5; ++i) { if (ev[i] < 0.) { System.out.format("Event %d, eigenvalue %d of covariance is negative!", iTrial, i); + badCov = true; } } + if (badCov) nBadCov++; if (iTrial < 10) { Vec evV = new Vec(5,ev); evV.print("Eigenvalues of covariance"); @@ -1107,6 +1111,7 @@ else if (kF.sites.indexOf(site) == kF.sites.size()-1) { } } + System.out.format("Number of track fits with bad covariance matrix = %d\n", nBadCov); long grdef = rnd.nextLong(); System.out.format("New seed = %d\n", grdef); timestamp = Instant.now(); diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java index 50d4107fb5..388d09fb3e 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalTrack.java @@ -100,8 +100,8 @@ public class KalTrack { for (int idx=firstSite; idx<=lastSite; ++idx) { MeasurementSite site = SiteList.get(idx); if (site.aS == null) { // This should never happen - logger.log(Level.WARNING, String.format("Event %d: site is missing smoothed state vector for layer %d detector %d", - eventNumber, site.m.Layer, site.m.detector)); + logger.log(Level.SEVERE, String.format("Event %d: site of track %d is missing smoothed state vector for layer %d detector %d", + eventNumber, ID, site.m.Layer, site.m.detector)); logger.log(Level.WARNING, site.toString("bad site")); bad = true; continue; diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java index 7f6778c375..363e136289 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java @@ -118,7 +118,7 @@ static Vec getField(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) { static double [] getFielD(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) { // Field map for stand-alone running - if (FieldMap.class.isInstance(hpsFm)) { return ((FieldMap) (hpsFm)).getField(kalPos); } + if (FieldMap.class.isInstance(hpsFm)) return ((FieldMap) (hpsFm)).getField(kalPos); // Standard field map for running in hps-java //System.out.format("Accessing HPS field map for position %8.3f %8.3f %8.3f\n", kalPos.v[0], kalPos.v[1], kalPos.v[2]); @@ -130,6 +130,8 @@ static Vec getField(Vec kalPos, org.lcsim.geometry.FieldMap hpsFm) { } else { if (hpsPos[1] > 70.0) hpsPos[1] = 70.0; // To avoid getting a field returned that is identically equal to zero if (hpsPos[1] < -70.0) hpsPos[1] = -70.0; + if (hpsPos[0] < -225.) hpsPos[0] = -225.; + if (hpsPos[0] > 270.) hpsPos[0] = 270.; } double[] hpsField = hpsFm.getField(hpsPos); if (uniformB) { 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 e0f0b68e73..e637d9dd41 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 @@ -157,7 +157,7 @@ public void print() { mxTdif = 30.; // Maximum time difference of hits in a track firstLayer = 0; // First layer in the tracking system (2 for pre-2019 data) lowPhThresh = 0.25; // Residual improvement ratio necessary to use a low-ph hit instead of high-ph - seedCompThr = -1.; // Remove SeedTracks with all Helix params within relative seedCompThr . If -1 do not apply duplicate removal + seedCompThr = 0.05; // Remove SeedTracks with all Helix params within relative seedCompThr . If -1 do not apply duplicate removal // Load the default search strategies // Index 0 is for the bottom tracker (+z), 1 for the top (-z) @@ -169,42 +169,24 @@ public void print() { // A S A S A S A S A S A S A S top // 0,1,2,3,4,5,6,7,8,9,10,11,12,13 // S A S A S A S A S A S A S A bottom - final int[] list0 = {6, 7, 8, 9, 10}; - final int[] list1 = {4, 5, 6, 7, 8}; - final int[] list2 = {5, 6, 8, 9, 10}; - final int[] list3 = {5, 6, 7, 8, 10}; - final int[] list4 = { 3, 6, 8, 9, 10 }; - final int[] list5 = { 4, 5, 8, 9, 10 }; - final int[] list6 = { 4, 6, 7, 8, 9 }; - final int[] list7 = { 4, 6, 7, 9, 10 }; - final int[] list8 = { 2, 5, 8, 9, 12}; - final int[] list9 = { 8, 10, 11, 12, 13}; - final int[] list10 = {6, 9, 10, 11, 12}; - final int[] list11 = {6, 7, 9, 10, 12}; - final int[] list12 = {2, 3, 4, 5, 6}; - final int[] list13 = {2, 4, 5, 6, 7}; - final int[] list14 = {6, 7, 8, 10, 11}; - final int[] list15 = {1, 2, 3, 4, 6}; - final int[] list16 = {0, 2, 3, 4, 5}; - final int[] list17 = {0, 3, 4, 5, 6}; - lyrList[0].add(list0); - lyrList[0].add(list1); - lyrList[0].add(list2); - lyrList[0].add(list3); - lyrList[0].add(list4); - lyrList[0].add(list5); - lyrList[0].add(list6); - lyrList[0].add(list7); - lyrList[0].add(list8); - lyrList[0].add(list9); - lyrList[0].add(list10); - lyrList[0].add(list11); - lyrList[0].add(list12); - lyrList[0].add(list13); - lyrList[0].add(list14); - lyrList[0].add(list15); - lyrList[0].add(list16); - lyrList[0].add(list17); + addStrategy("000BBS0"); + addStrategy("00BBS00"); + addStrategy("00ASBS0"); + addStrategy("00ABSS0"); + addStrategy("0A0SBS0"); + addStrategy("00B0BS0"); + addStrategy("00SBSA0"); + addStrategy("00SBB00"); + addStrategy("00SBAS0"); + addStrategy("0SA0BS0"); + addStrategy("0000SBB"); + addStrategy("000SABS"); + addStrategy("0BBS000"); + addStrategy("0SBB000"); + addStrategy("000BSSA"); + addStrategy("ABSS000"); + addStrategy("SBB0000"); + addStrategy("SABS000"); maxListIter1 = 14; // The maximum index for lyrList for the first iteration beamSpot = new double[3]; @@ -231,22 +213,6 @@ public void print() { // beamSpot[1] = beamPosKal.v[1]; // beamSpot[2] = beamPosKal.v[2]; // } - - // Swap axial/stereo in list entries for the top tracker - for (int[] list: lyrList[0]) { - int [] listTop = new int[5]; - for (int i=0; i<5; ++i) { - listTop[i] = Swap[list[i]]; - } - for (int i=0; i<4; ++i) { - if (listTop[i] > listTop[i+1]) { // Sorting entries. No more than one swap should be necessary. - int tmp = listTop[i]; - listTop[i] = listTop[i+1]; - listTop[i+1] = tmp; - } - } - lyrList[1].add(listTop); - } } public void setLowPhThreshold(double cut) { @@ -489,6 +455,10 @@ public void setNumSeedIter1(int num) { maxListIter1 = n-1; } + private boolean addStrategy(String strategy) { + return addStrategy(strategy,"top") && addStrategy(strategy,"bottom"); + } + // Add a seed search strategy for the bottom or top tracker public boolean addStrategy(String strategy, String topBottom) { if (!(topBottom == "top" || topBottom == "bottom")) { 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 8b2674a2d0..bae61a3cac 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 @@ -1436,7 +1436,7 @@ private void filterTrack(TrackCandidate tkrCandidate, int lyrBegin, // layer on int rF; double [] tRange = {tkrCandidate.tMax - kPar.mxTdif, tkrCandidate.tMin + kPar.mxTdif}; 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 if hitno=-1 + 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); if (rF > 0) { if (m.hits.get(newSite.hitID).tracks.size() > 0) tkrCandidate.nTaken++; @@ -1460,7 +1460,7 @@ private void filterTrack(TrackCandidate tkrCandidate, int lyrBegin, // layer on return; } } else { - boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made if hitno=-1 + boolean checkBounds = imod < moduleList.get(lyr).size() - 1; // Note: boundary check is not made on last module of layer rF = newSite.makePrediction(prevSite.aF, prevSite.m, hitno, tkrCandidate.nTaken <= kPar.mxShared, pickUp, checkBounds, tRange, trial); if (rF > 0) { if (m.hits.get(newSite.hitID).tracks.size() > 0) tkrCandidate.nTaken++; 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 9bebb66dbb..a7e4e68a56 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 @@ -43,10 +43,9 @@ import hep.aida.IHistogram1D; import hep.aida.IHistogramFactory; import hep.physics.vec.Hep3Vector; + /** * Histograms and plots for Kalman Filter pattern recognition development - * @author Robert Johnson - * */ class KalmanPatRecPlots { private KalmanInterface KI; @@ -93,10 +92,9 @@ class KalmanPatRecPlots { // arguments to histogram1D: name, nbins, min, max aida.histogram1D("Kalman number of tracks", 10, 0., 10.); - aida.histogram1D("Kalman Track Chi2", 50, 0., 100.); - aida.histogram1D("Kalman Track Chi2, >=12 hits", 50, 0., 100.); + aida.histogram1D("Kalman Track Chi2", 100, 0., 200.); + aida.histogram1D("Kalman Track Chi2, >=12 hits", 100, 0., 200.); aida.histogram1D("Kalman Track simple Chi2, >=12 hits", 50, 0., 100.); - aida.histogram1D("Kalman Track Chi2 with >=12 good hits", 50, 0., 100.); aida.histogram2D("number tracks Kalman vs GBL", 20, 0., 5., 20, 0., 5.); aida.histogram1D("helix chi-squared at origin", 100, 0., 25.); aida.histogram1D("GBL track chi^2", 50, 0., 100.); @@ -184,6 +182,13 @@ class KalmanPatRecPlots { hnh = aida.histogram1D("MC number hits",15,0.,15.); hnhf = aida.histogram1D("MC number hits, found",15,0.,15.); pEff = new Efficiency(40,0.,0.1,"Track efficency vs momentum","momentum (GeV)","efficiency"); + aida.histogram1D("Bad/Number of hits on bad tracks", 20, 0., 20.); + aida.histogram1D("Bad/Chi-squared of bad tracks", 100, 0., 200.); + aida.histogram1D("Bad/drho of bad tracks", 50, -8., 8.); + aida.histogram1D("Bad/dz of bad tracks", 50, -4., 4.); + 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.); } void process(EventHeader event, List sensors, ArrayList[] kPatList, @@ -419,14 +424,24 @@ void process(EventHeader event, List sensors, ArrayList[] aida.histogram1D("Kalman track time range (ns)").fill(kTk.tMax - kTk.tMin); // Check the covariance matrix + boolean badCov = false; Matrix C = new Matrix(kTk.originCovariance()); EigenvalueDecomposition eED= new EigenvalueDecomposition(C); double [] e = eED.getRealEigenvalues(); for (int i=0; i<5; ++i) { if (e[i] < 0.) { logger.warning(String.format("Event %d, eigenvalue %d of covariance is negative!", event.getEventNumber(), i)); + System.out.format("Event %d, eigenvalue %d of covariance is negative for track %d!", event.getEventNumber(), i, kTk.ID); + badCov = true; } } + if (badCov || kTk.bad) { + aida.histogram1D("Bad/Number of hits on bad tracks").fill(kTk.nHits); + aida.histogram1D("Bad/Chi-squared of bad tracks").fill(kTk.chi2); + aida.histogram1D("Bad/drho of bad tracks").fill(kTk.originHelixParms()[0]); + aida.histogram1D("Bad/dz of bad tracks").fill(kTk.originHelixParms()[3]); + aida.histogram1D("Bad/momentum of bad tracks").fill(pMag); + } if (kTk.nHits >= 10) { DMatrixRMaj thisCov = kTk.SiteList.get(0).aS.helix.C; double ptInv1 = kTk.SiteList.get(0).aS.helix.a.v[2]; @@ -516,6 +531,7 @@ void process(EventHeader event, List sensors, ArrayList[] } aida.histogram1D("Kalman track number of shared hits").fill(nShared); aida.histogram1D("Kalman track number MC particles").fill(mcParts.size()); + // Which MC particle is the best match? int idBest = -1; int nMatch = 0; @@ -548,6 +564,10 @@ void process(EventHeader event, List sensors, ArrayList[] if (!goodHit) nBad++; } aida.histogram1D("Kalman number of wrong hits on track").fill(nBad); + if (badCov || kTk.bad) { + aida.histogram1D("Bad/Number of MC particles associated").fill(mcParts.size()); + aida.histogram1D("Bad/Number of wrong hits on track").fill(nBad); + } if (kTk.nHits >= 10) aida.histogram1D("Kalman number of wrong hits on track, >= 10 hits").fill(nBad); MCParticle mcBest = null; @@ -824,21 +844,29 @@ void process(EventHeader event, List sensors, ArrayList[] } //loop on GBL Tracks } //check if event has GBLTracks + //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 = {}; if (nPlotted < numEvtPlots) { // && sharedHitTrack) { - KI.plotKalmanEvent(outputGnuPlotDir, event, kPatList); - //KI.plotGBLtracks(outputGnuPlotDir, event); - nPlotted++; - } - - double maxTruErr = simHitRes(event); - if (maxTruErr < 0.02) { - for (int topBottom=0; topBottom<2; ++topBottom) { - for (KalTrack kTk : kPatList[topBottom]) { - if (kTk.nHits < 12) continue; - aida.histogram1D("Kalman Track Chi2 with >=12 good hits").fill(kTk.chi2); - } + boolean plotIt = false; + if (badEvents.length > 0) { + for (int i=0; i m.xExtent[1] + tol) { return -2; } - if (rLocal.v[1] < m.yExtent[0] - tol || rLocal.v[1] > m.yExtent[1] + tol) { return -2; } - } - double deltaE = 0.; // dEdx*thickness/ct; Vec origin = m.p.X(); @@ -197,6 +182,7 @@ int makePrediction(StateVector pS, SiModule mPs, int hitNumber, boolean sharingO } // Move pivot point to X0 to generate the predicted helix + // First we need the momentum direction to calculate how much silicon we pass through Vec pMom = pS.helix.Rot.inverseRotate(pS.helix.getMom(0.)); double XL; if (mPs == null) { @@ -216,6 +202,29 @@ int makePrediction(StateVector pS, SiModule mPs, int hitNumber, boolean sharingO } } aP = pS.predict(thisSite, X0, B, tB, origin, XL, deltaE); + + // This calculates the H corresponding to using aP in h( , , ). It is kind of trivial, because + // aP are helix parameters for a pivot right at the predicted intersection point, on the helix. + // Hence the prediction at that point does not depend on the helix parameters at all. + buildH(aP, H); + CommonOps_DDRM.mult(aP.helix.C, H, tempV); + double Rextrap = CommonOps_DDRM.dot(H, tempV); + + // Check whether the intersection is within the bounds of the detector, with some margin + // If not, then the pattern recognition may look in another detector in the layer. + // Don't do the check if the hit number is already specified. Also, the check will not + // be called for in the last sensor of the layer, because we need to run the Kalman prediction + // to this layer before moving to the next, even if there is no hit. + double tol = kPar.edgeTolerance + FastMath.sqrt(Rextrap); // Tolerance on the check, in mm + if (debug) System.out.format("MeasurementSite.predict: edge tolerance = %10.4f mm\n", tol); + Vec rLocal = null; + if (checkBounds && hitNumber < 0) { + Vec rGlobal = pS.helix.toGlobal(X0); // Transform from field coordinates to global coordinates + rLocal = m.toLocal(rGlobal); // Rotate into the detector coordinate system + if (rLocal.v[0] < m.xExtent[0] - tol || rLocal.v[0] > m.xExtent[1] + tol) return -2; + if (rLocal.v[1] < m.yExtent[0] - tol || rLocal.v[1] > m.yExtent[1] + tol) return -2; + } + if (debug) { pS.helix.a.print("original helix in MeasurementSite.makePrediction"); pS.helix.X0.print("original helix pivot point"); @@ -230,9 +239,6 @@ int makePrediction(StateVector pS, SiModule mPs, int hitNumber, boolean sharingO // Vec X02 = aP.atPhi(phi2); // X02.print("intersection in local coordinates from new helix"); // aP.toGlobal(X02).print("intersection in global coordinates from new helix"); - } - - if (debug) { System.out.format("MeasurementSite.makePrediction: old helix intersects plane at phi=%10.7f\n", phi); Vec rGlobalOld = pS.helix.toGlobal(pS.helix.atPhi(phi)); rGlobalOld.print("global intersection with old helix from measurementSite.makePrediction"); @@ -244,11 +250,6 @@ int makePrediction(StateVector pS, SiModule mPs, int hitNumber, boolean sharingO aP.mPred = h(pS, m, phi); - // This calculates the H corresponding to using aP in h( , , ). It is kind of trivial, because - // aP are helix parameters for a pivot right at the predicted intersection point, on the helix. - // Hence the prediction at that point does not depend on the helix parameters at all. - buildH(aP, H); - // Test of the H vector, by comparing with a numerical difference. If this is done with the H // calculated from aP, comparing with h calculated from aP, then the difference is always zero. // By comparing with the less trivial case of h and H calculated from pS, we verify here that @@ -380,8 +381,7 @@ int makePrediction(StateVector pS, SiModule mPs, int hitNumber, boolean sharingO System.out.format("MeasurementSite.makePrediction: intersection with new helix is at phi=%10.7f, z=%10.7f\n", phi2, mPred2); } - CommonOps_DDRM.mult(aP.helix.C, H, tempV); - aP.R = m.hits.get(theHit).sigma * m.hits.get(theHit).sigma + CommonOps_DDRM.dot(H, tempV); + aP.R = m.hits.get(theHit).sigma * m.hits.get(theHit).sigma + Rextrap; chi2inc = aP.r * aP.r / aP.R; double cutVal = Math.abs(aP.r / m.hits.get(theHit).sigma); 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 52f14fe2c6..609446d267 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 @@ -148,6 +148,7 @@ StateVector predict(int newSite, Vec pivot, double B, Vec t, Vec originPrime, do } else { aPrime.helix.a = this.helix.pivotTransform(pivot, deltaEoE); } + if (debug) { // drho and dz are indeed always zero here aPrime.helix.a.print("StateVector predict: pivot transformed helix; should have zero drho and dz"); helix.a.print("old helix");