From 5686a0588bc40689ae8d126766e44c8fee3c7bf5 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 22 Sep 2025 12:33:36 -0700 Subject: [PATCH 1/3] clean up a bunch of commented code --- .../tracking/kalman/KalmanInterface.java | 58 +------------------ 1 file changed, 2 insertions(+), 56 deletions(-) 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 2449118e8..a9939a4db 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 @@ -101,7 +101,6 @@ public class KalmanInterface { private int nBigEvents; private int eventNumber; private static double target_pos = -999.9; - private static boolean addTrackStateAtTarget = false; private double[] beamPosition = null; private static final boolean debug = false; @@ -653,27 +652,10 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { ts = null; int loc = TrackState.AtOther; - //HpsSiSensor hssd = (HpsSiSensor) moduleMap.get(site.m).getSensor(); - //int lay = hssd.getMillepedeId(); - // System.out.printf("ssp id %d \n", hssd.getMillepedeId()); - if (i == firstHit_idx) { loc = TrackState.AtFirstHit; } else if (i == lastHit_idx) - loc = TrackState.AtLastHit; - - /* - //DO Not att the missing layer track states yet. - if (storeTrackStates) { - for (int k = 1; k < lay - prevID; k++) { - // uses new lcsim constructor - BaseTrackState dummy = new BaseTrackState(dummyCounter); - newTrack.getTrackStates().add(dummy); - dummyCounter--; - } - prevID = lay; - } - */ + loc = TrackState.AtLastHit; if (loc == TrackState.AtFirstHit || loc == TrackState.AtLastHit || storeTrackStates) { ts = createTrackState(site, loc, true, saveTrackStateAtIntercept); @@ -694,24 +676,6 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { double BAtEcal = BfieldAtEcal.mag(); double alphaAtEcal = conFac/ BAtEcal; - /* - - DMatrixRMaj ecalCov = new DMatrixRMaj(getCovarianceFromHelix(helixAtEcal)); - double[] ecalParams = helixAtEcal.a.v.clone(); - double[] ecalLCSimParams = getLCSimParams(ecalParams, alphaAtEcal); - double[] ecalLCSimCov = getLCSimCov(ecalCov, alphaAtEcal).asPackedArray(true); - double[] refAtEcal = localKalToHps(helixAtEcal.origin); - if(debug)System.out.println(this.getClass().getName()+":: reference to TrackState @ ecal = "+ - +refAtEcal[0]+", "+refAtEcal[1]+","+ refAtEcal[2]); - TrackState ts_ecal_helix = new BaseTrackState(ecalLCSimParams,refAtEcal, ecalLCSimCov, TrackState.AtCalorimeter, BAtEcal); - if(debug)System.out.println(this.getClass().getName()+":: Uncorrected track state: curvature = "+ts_ecal_helix.getOmega() - +" bField = "+ts_ecal_helix.getBLocal()+" momentum z = "+ts_ecal_helix.getMomentum()[0]); - System.out.println("Helix at ECal from kT.getHelixAtPlane and by hand conversions "); - System.out.println(ts_ecal_helix.toString()); - */ - //newTrack.getTrackStates().add(ts_ecal_helix); - - // this is how createTrackState (above) goes from measurement to TrackState // from the measurements. It calls toHPShelix. // the helix state here must already be propagated @@ -719,15 +683,7 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { TrackState ts_ecal=helixAtEcal.toTrackState(alphaAtEcal, ecalPlane, TrackState.AtCalorimeter, saveTrackStateAtIntercept); newTrack.getTrackStates().add(ts_ecal); - //if(debug)System.out.println("Helix at ECal from helix.toTrackState"); - // if(debug)System.out.println(ts_toTrackState.toString()); - // System.out.println("Helix at ECal from helix.toTrackState"); - // System.out.println(ts_toTrackState.toString()); - - // Extrapolate to the ECAL and make a new trackState there. - //BaseTrackState ts_ecal = new BaseTrackState(); - //ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber); - + Vec targetFace=origin; Plane targetPlane = new Plane(targetFace, new Vec(0., 1., 0.)); HelixState helixAtTarget=kT.getHelixAtPlane(targetPlane,saveTrackStateAtIntercept); // this propagates (via RK) the helix to the plane @@ -739,16 +695,6 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) { if(debug)System.out.println("Helix at Target from helix.toTrackState"); if(debug)System.out.println(ts_target.toString()); newTrack.getTrackStates().add(ts_target); - // Extrapolate to Target and make a new trackState there. - /* - BaseTrackState ts_target = new BaseTrackState(); - if (target_pos != -999.9 && addTrackStateAtTarget){ - ts_target = TrackUtils.getTrackExtrapAtTargetRK(newTrack, target_pos, beamPosition, fM, 0); - if (ts_target != null){ - newTrack.getTrackStates().add(ts_target); - } - } - */ // other track properties newTrack.setChisq(kT.chi2); From c733f6dc675450d212d704f82756dd5f102e7f5b Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Tue, 23 Sep 2025 15:23:12 -0700 Subject: [PATCH 2/3] add hole accounting to subdetectorHitNumbers as -1 --- .../hps/recon/tracking/kalman/KalTrack.java | 53 +++++++++------- .../tracking/kalman/KalmanInterface.java | 62 ++++++++++++++++++- .../tracking/kalman/MeasurementSite.java | 8 +++ 3 files changed, 97 insertions(+), 26 deletions(-) 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 72baaf0b8..4c1fef118 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 @@ -295,27 +295,49 @@ public Pair unbiasedIntersect(MeasurementSite site, boolean lo Vec intPnt = HelixState.atPhi(site.aS.helix.X0, aStar, phiInt, site.aS.helix.alpha); // Transform the intersection back into the global coordinates Vec globalInt = site.aS.helix.toGlobal(intPnt); + // System.out.println("KalTrack::unbiasedIntersect global x = "+globalInt.v[0]+ + // " y = "+globalInt.v[1]+ + // " z = "+globalInt.v[2]); + if (!local) { - new Pair<>(new Double[]{globalInt.v[0], globalInt.v[1], globalInt.v[2]}, 999); //I didn't include global variance since I'm lazy + new Pair<>(new Double[]{globalInt.v[0], globalInt.v[1], globalInt.v[2]}, 999); //I didn't include global variance since I'm lazy } // Transform the intersection point to the local sensor system Vec localInt = site.m.toLocal(globalInt); CommonOps_DDRM.mult(Cstar, site.H, tempV); Double varUmeas = CommonOps_DDRM.dot(site.H, tempV); return new Pair<>(new Double[]{localInt.v[0], localInt.v[1], localInt.v[2]}, varUmeas); } else { - if(debug)System.out.println("no phi-interect with this layer found!"); + if(debug)System.out.println("KalTrack::unbiasedIntersect no phi-interect with this layer found!"); } return new Pair<>(new Double[]{666., 666., 666.}, 666.); } + public boolean isTrackHole(int layer){ + if (lyrMap == null) { + makeLyrMap(); + } + if (lyrMap.containsKey(layer)) { + MeasurementSite ms=lyrMap.get(layer); + return isTrackHole(ms); + }else{ + return false; + } + } + public boolean isTrackHole(MeasurementSite ms){ + if(ms.hitID >=0 && ms.m.hits.size()>0) + return false; //track has a hit in this layer + Pair lInt=unbiasedIntersect(ms,true); + boolean trkInSensor=ms.isInSensor(lInt.getFirstElement()); + return trkInSensor; + } private void makeLyrMap() { lyrMap = new HashMap(nHits); for (MeasurementSite site : SiteList) { lyrMap.put(site.m.Layer, site); } } - + private void makeMillipedeMap() { millipedeMap = new HashMap(nHits); for (MeasurementSite site : SiteList) { @@ -512,7 +534,9 @@ public Pair unbiasedResidual(MeasurementSite site) { CommonOps_DDRM.mult(Cstar, site.H, tempV); varResid = sigma * sigma + CommonOps_DDRM.dot(site.H, tempV); - } + }else{ + System.out.println("KalTrack::unbiasedResidual phiInt is NaN"); + } } return new Pair(resid, varResid); } @@ -908,26 +932,7 @@ public HelixState originConstraint(double[] vtx, double[][] vtxCov) { Vec v = helixAtOrigin.toLocal(new Vec(3, vtx)); SquareMatrix Cov = helixAtOrigin.Rot.rotate(new SquareMatrix(3, vtxCov)); Vec X0 = helixAtOrigin.X0; - double phi = phiDOCA(helixAtOrigin.a, v, X0, alpha); - /* if (debug) { // Test the DOCA algorithm - Vec rDoca = HelixState.atPhi(X0, helixAtOrigin.a, phi, alpha); - System.out.format("originConstraint: phi of DOCA=%10.5e\n", phi); - rDoca.print(" DOCA point"); - double doca = rDoca.dif(v).mag(); - System.out.format(" Minimum doca=%10.7f\n", doca); - for (double p=phi-0.0001; p m.xExtent[1] + tol) return false; + if (rLocal[1] < m.yExtent[0] - tol || rLocal[1] > m.yExtent[1] + tol) return false; + return true; + } + // Comparator functions for sorting measurement sites by layer number static Comparator SiteComparatorUp = new Comparator() { public int compare(MeasurementSite s1, MeasurementSite s2) { From e4fba71e0e7de2e2e677bf5763d7d7c7e2ac59c7 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Thu, 25 Sep 2025 08:32:01 -0700 Subject: [PATCH 3/3] see if converting tab to spaces makes format better --- .../org/hps/recon/tracking/kalman/MeasurementSite.java | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/MeasurementSite.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/MeasurementSite.java index a306e0945..eb2e78237 100644 --- a/tracking/src/main/java/org/hps/recon/tracking/kalman/MeasurementSite.java +++ b/tracking/src/main/java/org/hps/recon/tracking/kalman/MeasurementSite.java @@ -76,7 +76,6 @@ String toString(String s) { } return str; } - MeasurementSite(int thisSite, SiModule data, KalmanParams kPar) { this.thisSite = thisSite; this.kPar = kPar; @@ -718,11 +717,10 @@ private void buildH(StateVector S, DMatrixRMaj H) { } public boolean isInSensor(Double[] rLocal){ - double tol = kPar.edgeTolerance; // Tolerance on the check, in mm - // double tol=0.0; - if (rLocal[0] < m.xExtent[0] - tol || rLocal[0] > m.xExtent[1] + tol) return false; - if (rLocal[1] < m.yExtent[0] - tol || rLocal[1] > m.yExtent[1] + tol) return false; - return true; + double tol = kPar.edgeTolerance; // Tolerance on the check, in mm + if (rLocal[0] < m.xExtent[0] - tol || rLocal[0] > m.xExtent[1] + tol) return false; + if (rLocal[1] < m.yExtent[0] - tol || rLocal[1] > m.yExtent[1] + tol) return false; + return true; } // Comparator functions for sorting measurement sites by layer number