Patches to KF-based Alignment Pipeline

tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java

@@ -30,6 +30,7 @@
 import org.lcsim.event.base.BaseTrackState;
 import org.lcsim.event.TrackerHit;
 import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
 import org.lcsim.geometry.Detector;
 import org.lcsim.util.Driver;
 import org.lcsim.util.aida.AIDA;
@@ -152,11 +153,17 @@ protected void detectorChanged(Detector detector) {
     public void process(EventHeader event) {
         List<Track> tracks = event.get(Track.class, trackCollectionName);
 
+        /*
+        System.out.print("GBLOutputDriver tracks (" + trackCollectionName + ") N hits: ");
+        for (Track trk : tracks) {
+          System.out.print(trk.getTrackerHits().size()+" ");
+        }
+        System.out.println();
+         */
+
         int TrackType = 0;
-
         if (trackCollectionName.contains("Kalman") || trackCollectionName.contains("KF")) {
             TrackType = 1;
-            nHits = 2*nHits;
             //System.out.println("PF:: DEBUG :: Found Kalman Tracks in the event");
         }
 
@@ -187,8 +194,13 @@ public void process(EventHeader event) {
             //System.out.println("Track passed chi2");
 
             //Remove tracks with less than 10 hits
-            if (trk.getTrackerHits().size() < nHits)
+            if ((TrackType == 0 && trk.getTrackerHits().size() < nHits) 
+                    || (TrackType == 1 && trk.getTrackerHits().size() < 2*nHits)) {
+                System.out.println("WARNING:: "+trk.getClass().getSimpleName()
+                        +" got to GBLOutputDriver with "+trk.getTrackerHits().size()+" hits"
+                        +" which is below the cut that should have been already applied.");
                 continue;
+            }
 
 
             //System.out.println("Track passed hits");
@@ -414,11 +426,16 @@ private void doBasicGBLtrack(Track trk, Map<HpsSiSensor, TrackerHit> sensorHits)
 
         FillGBLTrackPlot(trkpFolder+"Chi2",isTop,charge,trk.getChi2());
         FillGBLTrackPlot(trkpFolder+"Chi2_vs_p",isTop,charge,trackp,trk.getChi2());
-
-
 
-        aidaGBL.histogram1D(trkpFolder+"nHits" + isTop).fill(trk.getTrackerHits().size());
-        aidaGBL.histogram1D(trkpFolder+"nHits" + isTop+charge).fill(trk.getTrackerHits().size());
+        // deduce multiplication factor for ST-started GBL tracks
+        int nhits = trk.getTrackerHits().size();
+        if (nhits > 0 && trk.getTrackerHits().get(0) instanceof HelicalTrackCross) {
+            // tracks created from cross hits have 2 measurments per hit instead of only
+            // one so we have to double that count for an equal comparison
+            nhits *= 2;
+        }
+        aidaGBL.histogram1D(trkpFolder+"nHits" + isTop).fill(nhits);
+        aidaGBL.histogram1D(trkpFolder+"nHits" + isTop+charge).fill(nhits);
 
         Hep3Vector beamspot = CoordinateTransformations.transformVectorToDetector(TrackUtils.extrapolateHelixToXPlane(trackState, 0));
         if (debug)

tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectoryJna.java

@@ -30,6 +30,7 @@ Pointer GblTrajectoryCtorPtrComposed(Pointer [] points_1, int npoints_1, double
         void GblTrajectory_fit(Pointer self, DoubleByReference Chi2, IntByReference Ndf, DoubleByReference lostWeight, char [] optionList, int aLabel);
         void GblTrajectory_addPoint(Pointer self, Pointer point);
         int  GblTrajectory_isValid(Pointer self);
+        void GblTrajectory_delete(Pointer self);
         int  GblTrajectory_getNumPoints(Pointer self);
         void GblTrajectory_printTrajectory(Pointer self, int level);
         void GblTrajectory_printData(Pointer self);
@@ -62,7 +63,9 @@ public GblTrajectoryJna(List<GblPointJna> points, int flagCurv, int flagU1dir, i
         }
 
         self = GblTrajectoryInterface.INSTANCE.GblTrajectoryCtorPtrArray(ppoints, points.size(), flagCurv, flagU1dir, flagU2dir);
-
+        if (self == Pointer.NULL)
+            System.out.println("Failed generating trajectory");
+
     }
 
     //Simple trajectory constructor with seed 
@@ -157,7 +160,12 @@ public void printPoints(int level) {
         GblTrajectoryInterface.INSTANCE.GblTrajectory_printPoints(self,level);
     }
 
-
+    //Call delete on the underlying objects
+    public void delete() {
+        GblTrajectoryInterface.INSTANCE.GblTrajectory_delete(self);
+    }
+
+
     public void getMeasResults(int aLabel, int numData[], List<Double> aResiduals,List<Double> aMeasErrors, List<Double> aResErrors, List<Double> aDownWeights) {
 
         double[] d_aResiduals  = new double[2];

tracking/src/main/java/org/hps/recon/tracking/gbl/KalmanToGBLDriver.java

@@ -141,11 +141,17 @@ protected void process(EventHeader event) {
             //I'm using the OLD way to refit for the moment.
 
             FittedGblTrajectory fitGbl_traj = HpsGblRefitter.fit(list_kfSCDs, bfac, false);
-            GblTrajectory gbl_fit_trajectory =  fitGbl_traj.get_traj();
+
+            /*
+            System.out.println("DEBUG::Tom::KalmanToGBLDriver - converted KF track to GBL track with "
+                + gbl_fit_trajectory.getNumPoints() + " hits");
+             */
 
 
             // Compute the residuals
-            if (computeGBLResiduals) { 
+            // sometimes the GBL trajectory isn't fitted properly and null is returned
+            // so we should make sure that the trajectory was fit before continuing
+            if (computeGBLResiduals && fitGbl_traj != null) { 
 
                 TrackResidualsData resData  = GblUtils.computeGblResiduals(trk, fitGbl_traj);
                 trackResidualsCollection.add(resData);
@@ -165,6 +171,7 @@ protected void process(EventHeader event) {
             // Get the derivatives
 
             /*
+            GblTrajectory gbl_fit_trajectory =  fitGbl_traj.get_traj();
             for (GblData gbldata : gbl_fit_trajectory.getTrajData()) {
             float vals[] = new float[2];
                 List<Integer> indLocal = new ArrayList<Integer>();

tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java

@@ -81,7 +81,28 @@ public static Pair<Track, GBLKinkData> makeCorrectedTrack(FittedGblTrajectory fi
         // Create a new SeedTrack
         BaseTrack trk = new BaseTrack();
 
-        // Add the hits to the track
+        /*
+         * we could look into converting the points listed with the GBL
+         * trajectory into the actual TrackerHits in the GBL track.
+         * This was confusing so I didn't do it.
+        // Add the hits to the track, converting them to 1D hits
+        // that inherit from TrackerHit
+        for (GblPoint pt : fittedGblTrajectory.get_traj().getSinglePoints()) {
+           * GblPoint has
+           *  - P2P Jacobian
+           *  - Derivative matrices
+           *  - Position projected onto something?
+          SiTrackerHitStrip1D hit = new SiTrackerHitStrip1D(
+              ,// Hep3Vector position
+              ,// SymmetricMatrix covariance matrix
+              ,// double energy
+              ,// double time
+              ,// List<RawTrackerHit> raw hits
+              // TrackerHitType decoded type
+              );
+          trk.addHit(hit);
+        }
+        */
         for (TrackerHit hit : hitsOnTrack) {
             trk.addHit(hit);
         }

tracking/src/main/java/org/hps/recon/tracking/gbl/SimpleGBLTrajAliDriver.java

@@ -111,7 +111,7 @@ public class SimpleGBLTrajAliDriver extends Driver {
     private String trackResidualsRelColName = "TrackResidualsGBLRelations";
     private String rawHitCollectionName = "SVTRawTrackerHits";
     private String gblStripClusterDataRelations = "KFGBLStripClusterDataRelations";
-        
+
     private double bfield;
     private FieldMap bFieldMap;
     private final MultipleScattering _scattering = new MultipleScattering(new MaterialSupervisor());
@@ -345,6 +345,10 @@ public void setEnableStandardCuts(boolean val) {
         enableStandardCuts = val;
     }
 
+    public void setNHitsCut(int val) {
+        nHitsCut = val;
+    }
+
     public void setMinMom(double val) {
         minMom = val;
     }
@@ -497,6 +501,8 @@ protected void process(EventHeader event) {
             //System.out.println("PF:: DEBUG :: Found Kalman Tracks in the event");
         }
 
+        //System.out.println("DEBUG::Tom::Deduced a track type of "+TrackType);
+
         //If using Seed Tracker, get the hits from the event
         if (TrackType == 0) {
             hitToStrips = TrackUtils.getHitToStripsTable(event,helicalTrackHitRelationsCollectionName);
@@ -539,29 +545,40 @@ protected void process(EventHeader event) {
 
 
             if (enableAlignmentCuts) {
-
-
+                //System.out.println("DEBUG::Tom::alignment cuts enabled...");
                 //Get the track parameters
                 double[] trk_prms = track.getTrackParameters();
                 double tanLambda = trk_prms[BaseTrack.TANLAMBDA];
 
                 //Momentum cut: 3.8 - 5.2
                 Hep3Vector momentum = new BasicHep3Vector(track.getTrackStates().get(0).getMomentum());
+                //System.out.print("Enabled alignment cuts with hits cut = ");
                 //Kalman
-                if (TrackType == 1)
-                    nHitsCut = 2*nHitsCut;
+                int actualHitCut = nHitsCut;
+                if (TrackType == 1) {
+                    actualHitCut = 2*nHitsCut;
+                }
+                //System.out.println(actualHitCut);
 
                 if (momentum.magnitude() < minMom || momentum.magnitude() > maxMom) {
                     continue;
                 }
 
                 if (Math.abs(tanLambda) < maxtanL) {
-
                     continue;
                 }
 
                 //Align with tracks with at least 6 hits
-                if ((tanLambda > 0 && track.getTrackerHits().size() < nHitsCut) || (tanLambda < 0 && track.getTrackerHits().size() < nHitsCut))  {
+                if ((tanLambda > 0 && track.getTrackerHits().size() < actualHitCut) 
+                        || (tanLambda < 0 && track.getTrackerHits().size() < actualHitCut))  {
+                    continue;
+                }
+
+                if (TrackType == 1 && track.getTrackerHits().size() % 2 == 1) {
+                    // this is a KF track with an odd number of hits which /cannot/
+                    // be equivalent to a GBL track so we are going to skip it
+                    // in a future where KF-based alignment is the standard,
+                    // we probably want to remove this
                     continue;
                 }
 
@@ -578,6 +595,7 @@ protected void process(EventHeader event) {
                     else if (trackSide == 1 && !TrackUtils.isHoleTrack(track)) 
                         continue;
                 }
+                //System.out.println("DEBUG::Tom::Pass with " + track.getTrackerHits().size() + " hits");
             }
 
 
@@ -890,16 +908,11 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track))
 
                 if (!constrainedFit && !constrainedTanLFit && !constrainedPhi0Fit && !constrainedD0Fit && !constrainedZ0Fit) {
                     trajForMPII =  new GblTrajectoryJna(points_on_traj,1,1,1);
-                }
-
-                else {
-
+                } else {
                     trajForMPII = new GblTrajectoryJna(points_on_traj,1,seedPrecision,1,1,1);
                 }
 
-                if (debugAlignmentDs)
-                    trajForMPII.printData();
-
+                if (debugAlignmentDs) trajForMPII.printData();
 
                 //Fit the trajectory to get the Chi2
                 trajForMPII_unconstrained.fit(Chi2,Ndf, lostWeight,"");
@@ -908,8 +921,13 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track))
                 if (Chi2.getValue() / Ndf.getValue() > writeMilleChi2Cut)
                     continue;
 
-                if (writeMilleBinary)
+                if (writeMilleBinary) {
+                    /*
+                    System.out.println("DEBUG::Tom::Writing track with "
+                        + points_on_traj.size() + " hits to mille binary.");
+                     */
                     trajForMPII.milleOut(mille);
+                }
 
                 if (correctTrack) {                
 
@@ -923,6 +941,8 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track))
                     List<TrackerHit> allHthList = TrackUtils.sortHits(hth);
                     Pair<Track, GBLKinkData>  newTrack = MakeGblTracks.makeCorrectedTrack(fitTraj, TrackUtils.getHTF(track), allHthList, 0, bfield);
                     Track gblTrk = newTrack.getFirst();
+
+                    //System.out.println("DEBUG::Tom::Correct GBL track has "+gblTrk.getTrackerHits().size()+" hits");
 
                     if(computeGBLResiduals) {
 
@@ -986,6 +1006,7 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track))
                             }
 
 
+                            gbl_fit_traj_u.delete();
                         } // loop on sensor map
 
 
@@ -1000,20 +1021,29 @@ else if (trackSide == 1 && !TrackUtils.isHoleTrack(track))
 
 
 
-
-
-
+                    //System.out.println("Refitted track chi2 " + gblTrk.getChi2());
                     refittedTracks.add(gblTrk);
                     kinkDataCollection.add(newTrack.getSecond());
                     kinkDataRelations.add(new BaseLCRelation(newTrack.getSecond(), gblTrk));
                 }
 
+                trajForMPII.delete();
+                trajForMPII_unconstrained.delete();
+
             }// composite Alignment
 
         }//loop on tracks
 
 
         if (correctTrack) {
+            /*
+            System.out.print("Refitted tracks (" + outputCollectionName + ") N hits: ");
+            for (Track trk : refittedTracks) {
+              System.out.print(trk.getTrackerHits().size()+" ");
+            }
+            System.out.println();
+             */
+
             // Put the tracks back into the event and exit
             int flag = 1 << LCIOConstants.TRBIT_HITS;
             event.put(outputCollectionName, refittedTracks, Track.class, flag);