Add track state at target

steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullReconMC_KF_TrackClusterMatcher.lcsim

@@ -12,7 +12,7 @@
         <driver name="RfFitter"/>
 
         <!-- Ecal reconstruction drivers -->        
-        <!--<driver name="EcalRunningPedestal"/>-->
+        <driver name="EcalRunningPedestal"/>
         <driver name="EcalRawConverter" />
         <driver name="HitTimeSmear"/>
         <driver name="ReconClusterer" />
@@ -21,7 +21,7 @@
         <driver name="RawTrackerHitSensorSetup"/>
         <driver name="RawTrackerHitFitterDriver" />
         <driver name="TrackerHitDriver"/>
-        <driver name="HelicalTrackHitDriver"/>
+        <!--<driver name="HelicalTrackHitDriver"/>-->
 
         <!-- <driver name="MergeTrackCollections"/>
         <driver name="GBLRefitterDriver" />
@@ -44,9 +44,9 @@
         <driver name="RfFitter" type="org.hps.evio.RfFitterDriver"/>       
 
         <!-- Ecal reconstruction drivers -->
-        <!--<driver name="EcalRunningPedestal" type="org.hps.recon.ecal.EcalRunningPedestalDriver">
+        <driver name="EcalRunningPedestal" type="org.hps.recon.ecal.EcalRunningPedestalDriver">
             <logLevel>CONFIG</logLevel>
-        </driver>-->
+        </driver>
 
         <!--<driver name="EcalRawConverter" type="org.hps.recon.ecal.EcalRawConverter2Driver">
         </driver>-->
@@ -94,7 +94,7 @@
             <debug>false</debug>
         </driver>
         <driver name="TrackerHitDriver" type="org.hps.recon.tracking.DataTrackerHitDriver">
-            <neighborDeltaT>8.0</neighborDeltaT>
+            <neighborDeltaT>24.0</neighborDeltaT>
         </driver>
         <driver name="HelicalTrackHitDriver" type="org.hps.recon.tracking.HelicalTrackHitDriver">
             <debug>false</debug>
@@ -111,6 +111,9 @@
             <!-- <siHitsLimit>50</siHitsLimit> -->
             <!-- <seedCompThr>0.05</seedCompThr> -->
             <!-- <addResiduals>true</addResiduals> -->
+            <numEvtPlots>40</numEvtPlots>
+            <targetPosition>-4.3</targetPosition>
+            <addTrackStateAtTarget>true</addTrackStateAtTarge>
             <verbose>false</verbose>
           </driver>
 

steering-files/src/main/resources/org/hps/steering/recon/PhysicsRun2016FullRecon_KF_TrackClusterMatcher.lcsim

@@ -83,7 +83,7 @@
             <debug>false</debug>
         </driver>
         <driver name="TrackerHitDriver" type="org.hps.recon.tracking.DataTrackerHitDriver">
-            <neighborDeltaT>8.0</neighborDeltaT>
+            <neighborDeltaT>24.0</neighborDeltaT>
         </driver>
         <driver name="HelicalTrackHitDriver" type="org.hps.recon.tracking.HelicalTrackHitDriver">
             <debug>false</debug>
@@ -126,6 +126,8 @@
             <!-- <addResiduals>true</addResiduals> -->
             <numEvtPlots>40</numEvtPlots>
             <verbose>false</verbose>
+            <targetPosition>-4.3</targetPosition> 
+            <addTrackStateAtTarget>true</addTrackStateAtTarget>
           </driver>
 
         <driver name="LCIOWriter" type="org.lcsim.util.loop.LCIODriver">

tracking/src/main/java/org/hps/recon/tracking/TrackData.java

@@ -18,11 +18,14 @@ public class TrackData implements GenericObject {
     public static final int L1_ISOLATION_INDEX = 0;
     public static final int L2_ISOLATION_INDEX = 1;
     public static final int N_ISOLATIONS = 14;  //Default
-    public static final int N_TRACK_PARAMS = 4;
+    public static final int N_TRACK_PARAMS = 7;
     public static final int TRACK_TIME_INDEX = 0;
     public static final int PX_INDEX = 1;
     public static final int PY_INDEX = 2;
     public static final int PZ_INDEX = 3;
+    public static final int ORIGIN_BFY_INDEX = 4; //BFieldY at Origin
+    public static final int TARGET_BFY_INDEX = 5; //BFieldY at Target TrackState
+    public static final int ECAL_BFY_INDEX = 6;   //BFieldY at ECal TrackState
     public static final int TRACK_VOLUME_INDEX = 0;
     public static final String TRACK_DATA_COLLECTION = "TrackData";
     public static final String TRACK_DATA_RELATION_COLLECTION = "TrackDataRelations";
@@ -81,6 +84,25 @@ public TrackData(int trackVolume, float trackTime, double[] isolations, float[]
         this.floats = new float[]{trackTime,momentum[0],momentum[1],momentum[2]};
         this.ints = new int[]{trackVolume};
     }
+
+    /**
+     * Constructor 
+     *
+     * @param trackVolume : SVT volume associated with the track
+     * @param trackTime : The track time
+     * @param isolations : an array of doubles containing isolations for every
+     * sensor layer
+     * @param momentum   : an array of floats containing track momentum in the form (px,py,pz)
+     * @param origin_bfieldY : value of BfieldY at origin ref
+     * @param target_bfieldY : value of BfieldY at Target Track State
+     * @param ecal_bfieldY   : value of BfieldY at ECal Track State
+     */
+    public TrackData(int trackVolume, float trackTime, double[] isolations, float[] momentum, float origin_bfieldY, float target_bfieldY, float ecal_bfieldY) {
+
+        this.doubles = isolations;
+        this.floats = new float[]{trackTime,momentum[0],momentum[1],momentum[2],origin_bfieldY,target_bfieldY,ecal_bfieldY};
+        this.ints = new int[]{trackVolume};
+    }
 
     /**
      * Get isolation value for the hit in the given sensor layer.

tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java

@@ -755,6 +755,60 @@ public static BaseTrackState getTrackExtrapAtHodoRK(Track trk, FieldMap fM, int
         return getTrackExtrapAtHodoRK(trk, fM, 0, hodoLayer);
     }
 
+    public static BaseTrackState getTrackExtrapAtTargetRK(Track track, double target_z, double[] beamPosition, FieldMap fM, double stepSize) {
+
+        TrackState ts = track.getTrackStates().get(0);
+
+        //if track passed to extrapolateHelixToXPlane, uses first track state
+        //by default, else if trackstate is passed, uses trackstate params.
+        //Forms HTF, projects to location, then returns point on helix
+        Hep3Vector startPos = extrapolateHelixToXPlane(ts, 0.0);
+        Hep3Vector startPosTrans = CoordinateTransformations.transformVectorToDetector(startPos);
+        double distanceZ = target_z;
+        double charge = -1.0 * Math.signum(getR(ts));
+
+        //extrapolateTrackUsingFieldMapRK gets HTF of input track/trackstate
+        //if track is passed, defaults to first track state
+        //if trackstate is passed, uses trackstate params
+        org.hps.util.Pair<Hep3Vector, Hep3Vector> RKresults = extrapolateTrackUsingFieldMapRK(ts, startPosTrans, distanceZ, stepSize, fM);
+        double bFieldY = fM.getField(RKresults.getFirstElement()).y();
+        Hep3Vector posTrans = CoordinateTransformations.transformVectorToTracking(RKresults.getFirstElement());
+        Hep3Vector momTrans = CoordinateTransformations.transformVectorToTracking(RKresults.getSecondElement());
+
+        Hep3Vector finalPos = posTrans;
+        if (RKresults.getFirstElement().z() != target_z) {
+            Hep3Vector mom = RKresults.getSecondElement();
+            double dz = target_z - RKresults.getFirstElement().z();
+            double dy = dz * mom.y() / mom.z();
+            double dx = dz * mom.x() / mom.z();
+            Hep3Vector dPos = new BasicHep3Vector(dx, dy, dz);
+            finalPos = CoordinateTransformations.transformVectorToTracking(VecOp.add(dPos, RKresults.getFirstElement()));
+        }
+        bFieldY = fM.getField(CoordinateTransformations.transformVectorToDetector(finalPos)).y();
+        //params are calculated with respect to ref = {trackX, trackY, z=0)
+        double[] params = getParametersFromPointAndMomentum(finalPos, momTrans, (int) charge, bFieldY);
+        BaseTrackState bts = new BaseTrackState(params, bFieldY);
+        //reference point is set to track position in X Y Z
+        bts.setReferencePoint(finalPos.v());
+        //Define new reference point, to which track parameters are calc wrt
+        double[] newRef = {target_z, beamPosition[0], beamPosition[1]};
+        params = getParametersAtNewRefPoint(newRef, bts);
+        bts.setParameters(params, bFieldY);
+        //Reference point records final position of track.
+        //This does not hold the reference point to which the track params are
+        //calculated from 
+        bts.setReferencePoint(finalPos.v());
+        bts.setLocation(TrackState.LastLocation);
+
+        //Get covariance matrix at target. This does not use RK extrap, but
+        //simple change of reference point
+        SymmetricMatrix originCovMatrix = new SymmetricMatrix(5, ts.getCovMatrix(),true);
+        SymmetricMatrix covtrans = getCovarianceAtNewRefPoint(newRef, ts.getReferencePoint(), ts.getParameters(), originCovMatrix); 
+        bts.setCovMatrix(covtrans.asPackedArray(true));
+
+        return bts;
+    }
+
     /**
      * Extrapolate track to given position. For backwards compatibility.
      *
@@ -1869,6 +1923,10 @@ public static TrackState getTrackStateAtLocation(Track trk, int location) {
         return null;
     }
 
+    public static TrackState getTrackStateAtTarget(Track trk){
+        return getTrackStateAtLocation(trk, TrackState.LastLocation);
+    }
+
     public static TrackState getTrackStateAtECal(Track trk) {
         return getTrackStateAtLocation(trk, TrackState.AtCalorimeter);
     }

tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanInterface.java

@@ -79,6 +79,9 @@ public class KalmanInterface {
     private int maxHits;
     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;    
     private static final double SVTcenter = 505.57;
@@ -97,6 +100,18 @@ public void setSiHitsLimit(int limit) {
     public int getSiHitsLimit() {
         return _siHitsLimit;
     }
+
+    public void setTargetPosition(double target_pos){
+        this.target_pos = target_pos;   
+    }
+
+    public void setAddTrackStateAtTarget(boolean input){
+        this.addTrackStateAtTarget = input;
+    }
+
+    public void setBeamPosition(double[] beamPosition){
+        this.beamPosition = beamPosition;
+    }
 
     // Get the HPS tracker hit corresponding to a Kalman hit
     public TrackerHit getHpsHit(Measurement km) {
@@ -520,7 +535,7 @@ public List<GBLStripClusterData> createGBLStripClusterData(KalTrack kT) {
     public PropagatedTrackState propagateTrackState(TrackState stateHPS, double [] location, double [] direction) {
         return new PropagatedTrackState(stateHPS, location, direction, detPlanes, fM);
     }
-    
+
     // Create an HPS track from a Kalman track
     public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
         if (kT.SiteList == null) {
@@ -596,9 +611,18 @@ public BaseTrack createTrack(KalTrack kT, boolean storeTrackStates) {
         }
 
         // Extrapolate to the ECAL and make a new trackState there.
-
-        BaseTrackState ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber);
+        BaseTrackState ts_ecal = new BaseTrackState();
+        ts_ecal = TrackUtils.getTrackExtrapAtEcalRK(newTrack, fM, runNumber);
         newTrack.getTrackStates().add(ts_ecal);
+
+        // 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);

tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanPatRecDriver.java

@@ -19,6 +19,7 @@
 import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
 import org.hps.recon.tracking.MaterialSupervisor.SiStripPlane;
 import org.hps.recon.tracking.gbl.GBLStripClusterData;
+import org.hps.recon.tracking.TrackUtils;
 import org.hps.util.Pair;
 import org.lcsim.detector.tracker.silicon.HpsSiSensor;
 import org.lcsim.event.EventHeader;
@@ -56,6 +57,8 @@ public class KalmanPatRecDriver extends Driver {
     private KalmanParams kPar;
     private KalmanPatRecPlots kPlot;
     private static Logger logger;
+    private static double target_pos = -999.9;
+    private static boolean addTrackStateAtTarget = false;
 
     // Parameters for the Kalman pattern recognition that can be set by the user in the steering file:
     private ArrayList<String> strategies;     // List of seed strategies for both top and bottom trackers, from steering
@@ -130,6 +133,14 @@ public void setSiHitsLimit(int input) {
     public void setAddResiduals(boolean input) {
         addResiduals = input;
     }
+
+    public void setTargetPosition(double target_pos){
+        this.target_pos = target_pos;
+    }
+
+    public void setAddTrackStateAtTarget(boolean input){
+        this.addTrackStateAtTarget = input;
+    }
 
     @Override
     public void detectorChanged(Detector det) {
@@ -248,17 +259,28 @@ public void detectorChanged(Detector det) {
 
         // Setup optional usage of beam positions from database.
         final DatabaseConditionsManager mgr = DatabaseConditionsManager.getInstance();
-        if (useBeamPositionConditions && mgr.hasConditionsRecord("beam_positions")) {
-            logger.config("Using Kalman beam position from the conditions database");
+        double[] beamPositionArr = {beamPositionX, beamPositionY, beamPositionZ};
+        if (mgr.hasConditionsRecord("beam_positions")){
             BeamPositionCollection beamPositions = 
                     mgr.getCachedConditions(BeamPositionCollection.class, "beam_positions").getCachedData();
             BeamPosition beamPositionCond = beamPositions.get(0); 
-            if (!useFixedVertexZPosition) kPar.setBeamSpotY(beamPositionCond.getPositionZ());  
-            else logger.config("Using fixed Kalman beam Z position: " + kPar.beamSpot[1]);
-            kPar.setBeamSpotX(beamPositionCond.getPositionX());   // Includes a transformation to Kalman coordinates
-            kPar.setBeamSpotZ(-beamPositionCond.getPositionY());
-        } else {
-            logger.config("Using Kalman beam position from the steering file or default");
+            beamPositionArr[0] = beamPositionCond.getPositionX();
+            beamPositionArr[1] = beamPositionCond.getPositionY();
+            beamPositionArr[2] = beamPositionCond.getPositionZ();
+            System.out.println("beamPosition[0]: "+beamPositionArr[0]);
+            System.out.println("beamPosition[1]: "+beamPositionArr[1]);
+            System.out.println("beamPosition[2]: "+beamPositionArr[2]);
+            if (useBeamPositionConditions) {
+                logger.config("Using Kalman beam position from the conditions database");
+                if (!useFixedVertexZPosition) kPar.setBeamSpotY(beamPositionCond.getPositionZ());  
+                else logger.config("Using fixed Kalman beam Z position: " + kPar.beamSpot[1]);
+                kPar.setBeamSpotX(beamPositionCond.getPositionX());   // Includes a transformation to Kalman coordinates
+                kPar.setBeamSpotZ(-beamPositionCond.getPositionY());
+            } 
+            else {
+                logger.config("Using Kalman beam position from the steering file or default");
+            }
+
         }
         logger.config("Using Kalman beam position [ Z, X, Y ]: " + String.format("[ %f, %f, %f ]",
                        kPar.beamSpot[0], -kPar.beamSpot[2], kPar.beamSpot[1]) + " in HPS coordinates.");      
@@ -268,6 +290,12 @@ public void detectorChanged(Detector det) {
         kPar.print();
 
         KI = new KalmanInterface(uniformB, kPar, fm);
+        //Track State at Target uses beam position as track param reference
+        KI.setBeamPosition(beamPositionArr);
+        if (target_pos != -999.9 && addTrackStateAtTarget) {
+            KI.setAddTrackStateAtTarget(addTrackStateAtTarget);
+            KI.setTargetPosition(target_pos);
+        }
         KI.setSiHitsLimit(siHitsLimit);
         KI.createSiModules(detPlanes);
         decoder = det.getSubdetector("Tracker").getIDDecoder();
@@ -398,6 +426,7 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
                 //double pt = Math.abs(1./ptInv_check);
 
                 outputFullTracks.add(KalmanTrackHPS);
+
                 List<GBLStripClusterData> clstrs = KI.createGBLStripClusterData(kTk);
                 if (verbose) {
                     for (GBLStripClusterData clstr : clstrs) {
@@ -431,9 +460,28 @@ private ArrayList<KalTrack>[] prepareTrackCollections(EventHeader event, List<Tr
                 momentum_f[0] = (float) momentum.x();
                 momentum_f[1] = (float) momentum.y();
                 momentum_f[2] = (float) momentum.z();
-
+
+                //Get Bfield at origin
+                double origin_bFieldY = fm.getField((new BasicHep3Vector(0.0,0.0,0.0))).y();
+
+                //Get Bfield at target
+                double target_bFieldY = -999.9;
+                if (TrackUtils.getTrackStateAtTarget(KalmanTrackHPS) != null)
+                {
+                    Hep3Vector target_pos = new BasicHep3Vector(TrackUtils.getTrackStateAtTarget(KalmanTrackHPS).getReferencePoint());
+                    target_bFieldY = fm.getField(CoordinateTransformations.transformVectorToDetector(target_pos)).y();
+
+                }
+                //Get Bfield at ecal
+                double ecal_bFieldY = -999.9;
+                if (TrackUtils.getTrackStateAtECal(KalmanTrackHPS) != null)
+                {
+                    Hep3Vector ecal_pos = new BasicHep3Vector(TrackUtils.getTrackStateAtECal(KalmanTrackHPS).getReferencePoint());
+                    ecal_bFieldY = fm.getField(CoordinateTransformations.transformVectorToDetector(ecal_pos)).y();
+                }
+
                 //Add the Track Data 
-                TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f);
+                TrackData KFtrackData = new TrackData(trackerVolume, (float) kTk.getTime(), qualityArray, momentum_f, (float) origin_bFieldY, (float) target_bFieldY, (float) ecal_bFieldY);
                 trackDataCollection.add(KFtrackData);
                 trackDataRelations.add(new BaseLCRelation(KFtrackData, KalmanTrackHPS));