Merge pull request #1115 from JeffersonLab/lcsim-4.5.0-fixtrackstates

bloodyyugo · web-flow · commit 314952083096 · 2025-07-15T11:46:25.000-06:00
lcsim 4.5.0 plus fixtrackstates
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
@@ -0,0 +1,25 @@
+name: build
+
+on:
+  workflow_dispatch:
+  pull_request:
+    branches: [master]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout source code
+        uses: actions/checkout@v3
+
+      - name: Set up JDK 21
+        uses: actions/setup-java@v1
+        with:
+          java-version: 21
+          server-id: github
+
+      - name: Build
+        run: |
+          mvn -B install -DskipTests=true
+        env:
+          GITHUB_TOKEN: ${{secrets.GITHUB_TOKEN }}
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
@@ -12,11 +12,14 @@ jobs:
       - name: Checkout source code
         uses: actions/checkout@v3
 
-      - name: Set up JDK 1.8
+      - name: Set up JDK 21
         uses: actions/setup-java@v1
         with:
-          java-version: 1.8
+          java-version: 21
+          server-id: github
 
       - name: Build and run tests
         run: |
           mvn -B install -T4
+        env:
+          GITHUB_TOKEN: ${{secrets.GITHUB_TOKEN }}
diff --git a/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java b/analysis/src/main/java/org/hps/analysis/MC/TrackToMCParticleRelationsDriver.java
@@ -137,7 +137,7 @@ protected void process(EventHeader event) {
                     truth_trk.getTrackStates().clear();
                     double[] ref = new double[] { 0., 0., 0. };
                     SymmetricMatrix cov = new SymmetricMatrix(5);
-                    TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield);
+                    TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield);
                     truth_trk.getTrackStates().add(stateIP);
                     truth_trk.setChisq(-1);
                     truth_trk.setNDF(-1);
diff --git a/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java b/analysis/src/main/java/org/hps/analysis/examples/StarterAnalysisDriver.java
@@ -147,7 +147,7 @@ private List<LCIOTrackAnalysis> processTracks(EventHeader event) {
 
         for (Track track : tracklist) {   //remember, these tracks are in the lcsim tracking frame!     
             BaseTrackState ts = (BaseTrackState) track.getTrackStates().get(0);
-            ts.computeMomentum(bfield);
+            ts.computeMomentum();
             tkanalList.add(new LCIOTrackAnalysis(track, hittomc, hittostrip, hittorotated));
         }
         return tkanalList;
diff --git a/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java b/analysis/src/main/java/org/hps/analysis/vertex/VertexDebugger.java
@@ -803,8 +803,8 @@ else if (daug.getPDGID() == -11)
             double[] posShiftPars = TrackUtils.getParametersAtNewRefPoint(newRef, posOldTS);
             SymmetricMatrix eleShiftCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, eleOldTS.getReferencePoint(), eleOldTS.getParameters(), new SymmetricMatrix(5, eleOldTS.getCovMatrix(), true));
             SymmetricMatrix posShiftCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, posOldTS.getReferencePoint(), posOldTS.getParameters(), new SymmetricMatrix(5, posOldTS.getCovMatrix(), true));
-            BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleShiftCov.asPackedArray(true), TrackState.AtIP, B_FIELD);
-            BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posShiftCov.asPackedArray(true), TrackState.AtIP, B_FIELD);
+            BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleShiftCov.asPackedArray(true), TrackState.AtPerigee, B_FIELD);
+            BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posShiftCov.asPackedArray(true), TrackState.AtPerigee, B_FIELD);
 //            BaseTrackState eleShiftTS = new BaseTrackState(eleShiftPars, newRef, eleOldTS.getCovMatrix(), TrackState.AtIP, B_FIELD);
 //            BaseTrackState posShiftTS = new BaseTrackState(posShiftPars, newRef, posOldTS.getCovMatrix(), TrackState.AtIP, B_FIELD);
             BilliorTrack electronBTrackShift = toBilliorTrack(eleShiftTS);
diff --git a/pom.xml b/pom.xml
@@ -19,7 +19,7 @@
     <skipSite>false</skipSite>
     <skipPlugin>true</skipPlugin>
     <skipCheckstyle>true</skipCheckstyle>
-    <lcsimVersion>4.4.0</lcsimVersion>
+    <lcsimVersion>4.5.0</lcsimVersion>
     <etVersion>14.1</etVersion>
     <evioVersion>4.4.6</evioVersion>
     <junitVersion>4.13.1</junitVersion>
diff --git a/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java b/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java
@@ -16,6 +16,7 @@
 import org.hps.recon.vertexing.BilliorVertex;
 import org.hps.recon.vertexing.BilliorVertexer;
 import org.hps.record.StandardCuts;
+import org.hps.recon.tracking.TrackStateUtils;
 import org.lcsim.detector.tracker.silicon.HpsSiSensor;
 import org.lcsim.event.EventHeader;
 import org.lcsim.event.RawTrackerHit;
@@ -509,12 +510,22 @@ private BilliorVertex fitVertex(Constraint constraint, ReconstructedParticle ele
         BilliorTrack electronBTrack = toBilliorTrack(electron.getTracks().get(0));
         BilliorTrack positronBTrack = toBilliorTrack(positron.getTracks().get(0));
 
-        // Create a vertex fitter from the magnetic field.
+	// Create a vertex fitter from the magnetic field.
         // Note that the vertexing code uses the tracking frame coordinates
         // HPS X => TRACK Y
         // HPS Y => TRACK Z
         // HPS Z => TRACK X
-        BilliorVertexer vtxFitter = new BilliorVertexer(bField);
+	//first get the field @ perigee reference from the first tracks state (doesn't matter which one)
+	//	double bLocal=(electron.getTracks().get(0)).getTrackStates().get(TrackState.AtPerigee).getBLocal();
+	double bLocal=TrackStateUtils.getTrackStatesAtLocation(electron.getTracks().get(0),TrackState.AtPerigee).get(0).getBLocal();
+	//	System.out.println("For Vertexer:  "+TrackStateUtils.getTrackStatesAtLocation(electron.getTracks().get(0),TrackState.AtPerigee).get(0).toString());
+	// if we are using GBL tracks (trackType=0), set bLocal to bField (i.e. at SVT center)
+	if(trackType==0)
+	    bLocal=bField;
+
+	//set up vertexer with field
+	//	System.out.println("track type = "+trackType+"; setting vertex field = "+bLocal); 
+        BilliorVertexer vtxFitter = new BilliorVertexer(bLocal);
         // TODO: The beam size should come from the conditions database.
         vtxFitter.setBeamSize(beamSize);
         vtxFitter.setBeamPosition(beamPositionToUse);
@@ -880,11 +891,13 @@ private List<BilliorTrack> shiftTracksToVertex(List<ReconstructedParticle> parti
         double[] newRef = {vtxPos.z(), vtxPos.x(), 0.0};//the  TrackUtils.getParametersAtNewRefPoint method only shifts in xy tracking frame
         List<BilliorTrack> newTrks = new ArrayList<BilliorTrack>();
         for (ReconstructedParticle part : particles) {
-            BaseTrackState oldTS = (BaseTrackState) part.getTracks().get(0).getTrackStates().get(0);
+	    BaseTrackState oldTS=(BaseTrackState)TrackStateUtils.getTrackStatesAtLocation(part.getTracks().get(0),TrackState.AtPerigee).get(0);
+	    //	    BaseTrackState oldTS = (BaseTrackState) part.getTracks().get(0).getTrackStates().get(TrackState.AtPerigee);
             double[] newParams = TrackUtils.getParametersAtNewRefPoint(newRef, oldTS);
             SymmetricMatrix newCov = TrackUtils.getCovarianceAtNewRefPoint(newRef, oldTS.getReferencePoint(), oldTS.getParameters(), new SymmetricMatrix(5, oldTS.getCovMatrix(), true));
             //mg...I don't like this re-casting, but toBilliorTrack only takes Track as input
-            BaseTrackState newTS = new BaseTrackState(newParams, newRef, newCov.asPackedArray(true), TrackState.AtIP, bField);
+	    //make the state with bfield = field at original z-position because TrackUtils.getParametersAtNewRefPoint does not change curvature!
+            BaseTrackState newTS = new BaseTrackState(newParams, newRef, newCov.asPackedArray(true), TrackState.AtPerigee, oldTS.getBLocal());
             BilliorTrack electronBTrackShift = this.toBilliorTrack(newTS);
             newTrks.add(electronBTrackShift);
         }
diff --git a/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java b/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java
@@ -153,6 +153,10 @@ public void setApplyClusterCorrections(boolean val) {
      * LCIO collection name for tracks.
      */
     protected String matcherTrackCollectionName = "GBLTracks";
+    /**
+     * track type:  Kalman = 1; GBL = 0  obtained from matcherTrackCollectionName
+     */
+    protected int trackType = 0; 
     /**
      * Track Cluster Algorithm set to Kalman or GBL Tracks
      */
@@ -446,6 +450,10 @@ protected void detectorChanged(Detector detector) {
         matcher.setTrackCollectionName(matcherTrackCollectionName);
         matcher.enablePlots(enableTrackClusterMatchPlots);
 
+	//set the track type; default is 0 (GBL)
+	if(matcherTrackCollectionName.contains("Kalman") || matcherTrackCollectionName.contains("KF"))
+	    trackType=1;
+	
         // Set the magnetic field parameters to the appropriate values.
         Hep3Vector ip = new BasicHep3Vector(0., 0., 500.0);
         bField = detector.getFieldMap().getField(ip).y();
diff --git a/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java b/recon/src/main/java/org/hps/recon/utils/TrackTruthRelationsDriver.java
@@ -1479,7 +1479,7 @@ protected void process(EventHeader event) {
                     truth_trk.getTrackStates().clear();
                     double[] ref = new double[] { 0., 0., 0. };
                     SymmetricMatrix cov = new SymmetricMatrix(5);
-                    TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield);
+                    TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield);
                     truth_trk.getTrackStates().add(stateIP);
                     truth_trk.setChisq(-1);
                     truth_trk.setNDF(-1);
@@ -1538,7 +1538,7 @@ public double[] getMCPTrackParameters(MCParticle mcp, double bfield){
         truth_trk.getTrackStates().clear();
         double[] ref = new double[] { 0., 0., 0. };
         SymmetricMatrix cov = new SymmetricMatrix(5);
-        TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtIP,bfield);
+        TrackState stateIP = new BaseTrackState(mcp_htf.parameters(),ref,cov.asPackedArray(true),TrackState.AtPerigee,bfield);
         truth_trk.getTrackStates().add(stateIP);
         truth_trk.setChisq(-1);
         truth_trk.setNDF(-1);
diff --git a/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java b/tracking/src/main/java/org/hps/recon/tracking/TrackTweakDriver.java
@@ -214,9 +214,9 @@ public void process(EventHeader event) {
 
             // Extrapolate the tweaked track to the face of the Ecal and get the
             // track state
-            TrackState stateIP = TrackUtils.getTrackStateAtLocation(track, TrackState.AtIP);
+            TrackState stateIP = TrackUtils.getTrackStateAtLocation(track, TrackState.AtPerigee);
             if (stateIP == null) {
-                throw new RuntimeException("IP track state for GBL track was not found");
+                throw new RuntimeException("Perigee track state for GBL track was not found");
             }
             TrackState stateEcalIP = TrackUtils.extrapolateTrackUsingFieldMap(stateIP, extStartPos, ecalPosition, stepSize, bFieldMap);
 
diff --git a/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java b/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
@@ -1750,18 +1750,18 @@ public static Double[] getIsolations(Track trk, RelationalTable hitToStrips, Rel
      * .
      */
     public static BaseTrackState extrapolateTrackUsingFieldMap(Track track, double startPositionX, double endPositionX, double stepSize, FieldMap fieldMap) {
-        TrackState stateAtIP = null;
+        TrackState stateAtPerigee = null;
         for (TrackState state : track.getTrackStates()) {
-            if (state.getLocation() == TrackState.AtIP) {
-                stateAtIP = state;
+            if (state.getLocation() == TrackState.AtPerigee) {
+                stateAtPerigee = state;
             }
         }
-        if (stateAtIP == null) {
-            throw new RuntimeException("No track state at IP was found so this function shouldn't be used.");
+        if (stateAtPerigee == null) {
+            throw new RuntimeException("No track state at Perigee was found so this function shouldn't be used.");
         }
 
         // Extrapolate this track state
-        return extrapolateTrackUsingFieldMap(stateAtIP, startPositionX, endPositionX, stepSize, fieldMap);
+        return extrapolateTrackUsingFieldMap(stateAtPerigee, startPositionX, endPositionX, stepSize, fieldMap);
     }
 
     /**
diff --git a/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java b/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
@@ -123,9 +123,14 @@ public static Pair<Track, GBLKinkData> makeCorrectedTrack(FittedGblTrajectory fi
 
         trk.setTrackParameters(correctedHelixParams.getFirst(), bfield);// hack to set the track charge
         trk.getTrackStates().clear();
-        TrackState stateIP = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
-        trk.getTrackStates().add(stateIP);
 
+	// Set the track state to AtPerigee to agree with new track state location definitions as of hps-lcsim version 4.5.0
+        TrackState statePerigee = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtPerigee, bfield);
+        trk.getTrackStates().add(statePerigee);
+	//also add state at the IP for backward compatibility, but it's the same as the perigee for gbl tracks 
+	TrackState stateIP = new BaseTrackState(correctedHelixParams.getFirst(), ref, correctedHelixParams.getSecond().asPackedArray(true), TrackState.AtIP, bfield);
+	trk.getTrackStates().add(stateIP);
+	
         if (!storeTrackStates) {
             // just store last state
             Pair<double[], SymmetricMatrix> correctedHelixParamsLast = fittedGblTrajectory.getCorrectedPerigeeParameters(helicalTrackFit, FittedGblTrajectory.GBLPOINT.LAST, bfield);
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
@@ -322,9 +322,11 @@ static Vec pivotTransform(Vec pivot, Vec a, Vec X0, double alpha, double deltaEo
 
         return new Vec(5, aP);
     }
-    
+    HelixState propagateRungeKutta(Plane pln, ArrayList<Double> yScat, ArrayList<Double> XL, org.lcsim.geometry.FieldMap fM, double [] arcLength){
+	return propagateRungeKutta(pln, yScat, XL, fM, arcLength, false); 
+    }
     // Propagate a helix by Runge-Kutta integration to an arbitrary plane
-    HelixState propagateRungeKutta(Plane pln, ArrayList<Double> yScat, ArrayList<Double> XL, org.lcsim.geometry.FieldMap fM, double [] arcLength) {
+    HelixState propagateRungeKutta(Plane pln, ArrayList<Double> yScat, ArrayList<Double> XL, org.lcsim.geometry.FieldMap fM, double [] arcLength, boolean pivotAtIntersect) {
         // pln   = plane to where the extrapolation is taking place in global coordinates.  
         //         The origin of pln will be the new helix pivot point in global coordinates and the origin of the B-field system.
         // yScat = input array of y values where scattering in silicon will take place. Only those between the start and finish points
@@ -393,6 +395,7 @@ HelixState propagateRungeKutta(Plane pln, ArrayList<Double> yScat, ArrayList<Dou
             helixAtTarget.print("helix with pivot at final plane");
             pln.X().print("origin of final field system in global system");
             Rot.print("rotation matrix from global to final field system");
+	    System.out.println("Return parameters with pivot at the intersection:: "+pivotAtIntersect); 
         }
 
         // The covariance matrix is transformed assuming a sequence of pivot transforms (not Runge Kutta)
@@ -407,14 +410,18 @@ HelixState propagateRungeKutta(Plane pln, ArrayList<Double> yScat, ArrayList<Dou
             System.out.println("    Errors: ");
             for (int i = 0; i < 5; ++i) { System.out.format(" %10.7f", FastMath.sqrt(C.get(i,i))); }
             System.out.println("\n");
-            newCovariance.print("transformed covariance");
+            newCovariance.print();
             System.out.println("    Errors: ");
             for (int i = 0; i < 5; ++i) { System.out.format(" %10.7f", FastMath.sqrt(newCovariance.get(i,i))); }
             System.out.println("\n");
         }
-        
-        HelixState newHelixState = new HelixState(helixAtTarget, new Vec(0.,0.,0.), pln.X(), newCovariance, Bmag, tB);
-        newHelixState.xPlaneRK = xPlane;
+        HelixState newHelixState=null; 
+	if(pivotAtIntersect)
+	    newHelixState = new HelixState(helixAtIntersect, new Vec(0.,0.,0.), xPlane, newCovariance, Bmag, tB);
+	else
+	    newHelixState = new HelixState(helixAtTarget, new Vec(0.,0.,0.), pln.X(), newCovariance, Bmag, tB);
+	
+	newHelixState.xPlaneRK = xPlane;
         return newHelixState;
     }
 
@@ -630,9 +637,10 @@ boolean helixStepper(double maxStep, ArrayList<Double> yScat, ArrayList<Double>
     // In the returned TrackState the reference point gets set to the point on the helix closest to the
     // original pivot point (e.g. the helix intersection with the plane of silicon).
     // The pivot of the returned TrackState is always the origin (0,0,0)
-    TrackState toTrackState(double alphaCenter, Plane pln, int loc) {
+    TrackState toTrackState(double alphaCenter, Plane pln, int loc, boolean pivotAtIntercept) {
+	//    TrackState toTrackState(double alphaCenter, Plane pln, int loc) {
         // See TrackState for the different choices for loc (e.g. TrackState.atOther)
-        return KalmanInterface.toTrackState(this, pln, alphaCenter, loc);    
+        return KalmanInterface.toTrackState(this, pln, alphaCenter, loc,  pivotAtIntercept);    
     }
      
     // Transform a helix from one pivot to another through a non-uniform B field in several steps
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
@@ -762,7 +762,7 @@ class HelixTest3 { // Program for testing the Kalman fitting code
                 else if (kF.sites.indexOf(site) == kF.sites.size()-1) {
                     loc = TrackState.AtLastHit;
                 }
-                TrackState ts = KI.createTrackState(site, loc, true);    
+                TrackState ts = KI.createTrackState(site, loc, true,true);    
                 states.add(ts);
             }
             TrackState lastState = null;
diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KFOutputDriver.java
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
diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanDriverHPS.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanDriverHPS.java
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
diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanKinkFit.java
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
diff --git a/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java b/tracking/src/main/java/org/hps/recon/tracking/kalman/PropagatedTrackState.java

Original file line number	Diff line number	Diff line change
`@@ -147,7 +147,7 @@ private List<LCIOTrackAnalysis> processTracks(EventHeader event) {`
`147`	`147`
`148`	`148`	`for (Track track : tracklist) { //remember, these tracks are in the lcsim tracking frame!`
`149`	`149`	`BaseTrackState ts = (BaseTrackState) track.getTrackStates().get(0);`
`150`		`- ts.computeMomentum(bfield);`
	`150`	`+ ts.computeMomentum();`
`151`	`151`	`tkanalList.add(new LCIOTrackAnalysis(track, hittomc, hittostrip, hittorotated));`
`152`	`152`	`}`
`153`	`153`	`return tkanalList;`