Fixtrackstates

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;

recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java

@@ -509,12 +509,14 @@ 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(); 
+        BilliorVertexer vtxFitter = new BilliorVertexer(bLocal);
         // TODO: The beam size should come from the conditions database.
         vtxFitter.setBeamSize(beamSize);
         vtxFitter.setBeamPosition(beamPositionToUse);
@@ -880,11 +882,12 @@ 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) 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);
         }

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

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;