Iss930 Kalman code for incorporating ECAL energy information into track fits

online-monitoring/src/main/java/org/hps/monitoring/drivers/trackrecon/KFSVTOpeningAlignment.java

@@ -280,7 +280,7 @@ protected void detectorChanged(Detector detector) {
         kPar = new KalmanParams();
         kPar.print();
 
-        KI = new KalmanInterface(false, kPar, fm);
+        KI = new KalmanInterface(false, kPar, detector, fm);
         KI.createSiModules(detPlanes);
 
 //        plotterParsBot = pfac.create("Bot Track Pars");

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

@@ -1,8 +1,11 @@
 package org.hps.recon.tracking.gbl;
 
+import java.io.IOException;
 import java.util.ArrayList;
 import java.util.List;
 import java.util.Set;
+import java.util.logging.Level;
+import java.util.logging.Logger;
 import java.util.Collections;
 import java.util.Comparator;
 
@@ -11,21 +14,28 @@
 //import org.lcsim.event.TrackState;
 
 import org.lcsim.event.Track;
-
+import org.lcsim.event.TrackState;
 import org.lcsim.event.EventHeader;
 import org.lcsim.util.Driver;
-import org.lcsim.geometry.Detector;
+import org.lcsim.util.aida.AIDA;
 
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+
+import org.lcsim.geometry.Detector;
 import org.lcsim.event.LCRelation;
 import org.lcsim.event.RelationalTable;
 import org.lcsim.event.base.BaseRelationalTable;
-
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
 import org.lcsim.constants.Constants;
+import org.apache.commons.math3.util.Pair;
 import org.hps.recon.tracking.TrackResidualsData;
+import org.hps.recon.tracking.TrackStateUtils;
 import org.lcsim.event.base.BaseLCRelation;
 
 //import org.lcsim.fit.helicaltrack.HelicalTrackFit;
 import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.TrackingReconstructionPlots;
 
 /**
  * A Driver which refits Kalman Tracks using GBL
@@ -42,19 +52,46 @@ public class KalmanToGBLDriver extends Driver {
     private String kinkDataCollectionName = GBLKinkData.DATA_COLLECTION;
     private String kinkDataRelationsName  = GBLKinkData.DATA_RELATION_COLLECTION;
     private Boolean _debug = false;
+    private Boolean _analysis = false;
     private Boolean constrainedBSFit = false;
     private double bfield;
     private Boolean computeGBLResiduals = true;
+    public AIDA aida;
+    private IHistogram1D hGBLurt, hGBLurb, hKFurt, hKFurb, hGBLrt, hGBLrb;
+    private IHistogram1D hDelPhi0, hDelTanL, hDelOmega, hDelD0, hDelZ0;
+    private IHistogramFactory ahf;
 
-    void setDebug(boolean val) {
+    public void setDebug(boolean val) {
         _debug = val;
     }
 
+    public void setAnalysis(boolean val) {
+        _analysis = val;
+    }
+
 
+    void definePlots() {
+        if (aida == null) aida = AIDA.defaultInstance();
+        aida.tree().cd("/");
+        ahf = aida.histogramFactory();
+        hGBLurt = aida.histogram1D("GBL unbiased residuals top detector", 100, -0.25, 0.25);
+        hGBLurb = aida.histogram1D("GBL unbiased residuals bottom detector", 100, -0.25, 0.25);
+        hGBLrt = aida.histogram1D("GBL biased residuals top detector", 100, -1., 1.);
+        hGBLrb = aida.histogram1D("GBL biased residuals bottom detector", 100, -1., 1.);
+        hKFurt = aida.histogram1D("KF unbiased residuals top detector", 100, -1., 1.);
+        hKFurb = aida.histogram1D("KF unbiased residuals bottom detector", 100, -1., 1.);
+        hDelPhi0 = aida.histogram1D("KF minus GBL phi0", 100, -0.02, 0.02);
+        hDelD0 = aida.histogram1D("KF minus GBL D0", 100, -5., 5.);
+        hDelOmega = aida.histogram1D("KF minus GBL Omega", 100, -3.E-4, 3.E-4);
+        hDelTanL = aida.histogram1D("KF minus GBL tanLamba", 100, -0.01, 0.01);
+        hDelZ0 = aida.histogram1D("KF minus GBL Z0", 100, -0.1, 0.1);
+    }
+
     @Override 
     protected void detectorChanged(Detector detector) {
         bfield = Math.abs(TrackUtils.getBField(detector).magnitude());
 
+        if (_analysis) definePlots();
     }
 
     @Override
@@ -69,19 +106,66 @@ protected void process(EventHeader event) {
         List<GBLKinkData> kinkDataCollection = new ArrayList<GBLKinkData>();
         List<LCRelation> kinkDataRelations = new ArrayList<LCRelation>();
 
+        List<Track> GBLtracks = null;
+        if (_analysis) {
+            if (event.hasCollection(Track.class, "GBLTracks")) {
+                GBLtracks = event.get(Track.class, "GBLTracks");
+                if (_debug) {
+                    System.out.format("KalmanToGBLDriver: number of input GBL tracks = %d\n", GBLtracks.size());
+                }
+            }
+        }
+
         //Get the track collection from the event
 
-        if (!event.hasCollection(Track.class, inputCollectionName)) 
+        if (!event.hasCollection(Track.class, inputCollectionName)) {
+            if (_debug) System.out.format("KalmanToGBLDriver: collection %s not found\n", inputCollectionName);
             return;
+        }
 
         List<Track> tracks = event.get(Track.class,inputCollectionName);
 
-        if (_debug)
-            System.out.println("Found tracks: " + inputCollectionName+" " + tracks.size());
-
+        String residualsKFcollectionName = "KFUnbiasRes";
+        List<TrackResidualsData> residualsKF = null;
+        if (event.hasCollection(TrackResidualsData.class, residualsKFcollectionName)) {
+            if (_debug) System.out.format("KalmanToGBLDriver: collection %s found\n", residualsKFcollectionName);
+            residualsKF = event.get(TrackResidualsData.class, residualsKFcollectionName);
+        }
+        String residualsRelationsKFcollectionName = "KFUnbiasResRelations";
+        List<LCRelation> residualsRelKF = null;
+        RelationalTable kfResidsRT = null;
+        if (event.hasCollection(LCRelation.class, residualsRelationsKFcollectionName)) {
+            if (_debug) System.out.format("KalmanToGBLDriver: collection %s found\n", residualsRelationsKFcollectionName);
+            residualsRelKF = event.get(LCRelation.class, residualsRelationsKFcollectionName);
+            kfResidsRT = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+            for (LCRelation relation : residualsRelKF) {
+                if (relation != null && relation.getFrom() != null && relation.getTo() != null) { 
+                    kfResidsRT.add(relation.getFrom(), relation.getTo());
+                }
+            }
+        }
+
+        if (_debug) System.out.println("KalmanToGBLDriver, found tracks: " + inputCollectionName+" " + tracks.size());
 
-        //Loop on Kalman Tracks 
+        RelationalTable kfSCDsRT = null;
+        List<LCRelation> kfSCDRelation = new ArrayList<LCRelation>();       
+        if (event.hasCollection(LCRelation.class,"KFGBLStripClusterDataRelations")) {
+            kfSCDsRT = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+            kfSCDRelation = event.get(LCRelation.class,"KFGBLStripClusterDataRelations");
+            for (LCRelation relation : kfSCDRelation) {
+                if (relation != null && relation.getFrom() != null && relation.getTo() != null) { 
+                    kfSCDsRT.add(relation.getFrom(), relation.getTo());
+                }
+            }   
+        } else {
+            System.out.println("null KFGBLStripCluster Data Relations."); 
+            return;
+        }
 
+        List<Track> newGBLtks = null;
+        if (_analysis) newGBLtks = new ArrayList<Track>(tracks.size());
+
+        //Loop on Kalman Tracks        
         for (Track trk : tracks ) {
 
             //Remove tracks with less than 10 hits
@@ -93,34 +177,12 @@ protected void process(EventHeader event) {
             //Remove tracks where there is high mis-tracking rate
             //TrackState trackState = trk.getTrackStates().get(0);
             //if (Math.abs(trackState.getTanLambda()) < 0.02)
-            //    continue;
-
-
-            //Get the GBLStripClusterData
-            RelationalTable kfSCDsRT = null;
-            List<LCRelation> kfSCDRelation = new ArrayList<LCRelation>();
-
-
-            if (event.hasCollection(LCRelation.class,"KFGBLStripClusterDataRelations")) {
-                kfSCDsRT = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
-                kfSCDRelation = event.get(LCRelation.class,"KFGBLStripClusterDataRelations");
-                for (LCRelation relation : kfSCDRelation) {
-                    if (relation != null && relation.getFrom() != null && relation.getTo() != null) { 
-                        kfSCDsRT.add(relation.getFrom(), relation.getTo());
-                    }
-                }   
-            } else {
-                System.out.println("null KFGBLStripCluster Data Relations."); 
-                return;
-            }
-
+            //    continue;                                 
+
             //Get the strip cluster data
             Set<GBLStripClusterData> kfSCDs = kfSCDsRT.allFrom(trk);
 
-
-            if (_debug) 
-                System.out.println("Kalman Strip Cluster data size: " + kfSCDs.size());
-
+            if (_debug) System.out.println("Kalman Strip Cluster data size: " + kfSCDs.size());            
 
             //Convert the set to a list for sorting it
 
@@ -142,6 +204,7 @@ protected void process(EventHeader event) {
 
             FittedGblTrajectory fitGbl_traj = HpsGblRefitter.fit(list_kfSCDs, bfac, false);
 
+
             /*
             System.out.println("DEBUG::Tom::KalmanToGBLDriver - converted KF track to GBL track with "
                 + gbl_fit_trajectory.getNumPoints() + " hits");
@@ -166,8 +229,7 @@ protected void process(EventHeader event) {
 
             }//computeGBLResiduals
 
-
-
+
             // Get the derivatives
 
             /*
@@ -191,6 +253,15 @@ protected void process(EventHeader event) {
             }
             */
 
+            // Create a GBL track from the fitted trajectory
+            if (_analysis) {
+                int trackType = 57;    // randomly chosen track type
+                TrackState atIP = TrackStateUtils.getTrackStateAtIP(trk);
+                HelicalTrackFit htkft = TrackUtils.getHTF(atIP);
+                Pair<Track, GBLKinkData> GBLtkrPair = MakeGblTracks.makeCorrectedTrack(fitGbl_traj, htkft, trk.getTrackerHits(), trackType, bfield, true);
+                newGBLtks.add(GBLtkrPair.getFirst());
+            }
+
         }// track loop
 
         if (computeGBLResiduals) {
@@ -200,6 +271,63 @@ protected void process(EventHeader event) {
             event.put(kinkDataRelationsName, kinkDataRelations, LCRelation.class, 0);
 
         }
+
+        // Analysis added by R. Johnson for debugging etc
+        if (_debug) {
+            if (residualsKF == null) System.out.format("KalmanToGBLDriver: residualsKF is null\n");
+            if (residualsRelKF == null) System.out.format("KalmanToGBLDriver: residualsRelKF is null\n");
+        }
+        if (_analysis && residualsKF != null && residualsRelKF != null) {
+
+            RelationalTable gblResidsRT = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+            for (LCRelation relation : trackResidualsRelations) {
+                if (relation != null && relation.getFrom() != null && relation.getTo() != null) { 
+                    gblResidsRT.add(relation.getFrom(), relation.getTo());
+                }
+            }
+            for (Track trk : tracks ) {
+                Track trkGBL = newGBLtks.get(tracks.indexOf(trk));
+                Set<TrackResidualsData> gblResids = gblResidsRT.allTo(trk);
+                TrackState atIP_GBL = TrackStateUtils.getTrackStateAtIP(trkGBL);
+                TrackState atIP_KF = TrackStateUtils.getTrackStateAtIP(trk);
+                hDelPhi0.fill(atIP_KF.getPhi() - atIP_GBL.getPhi());
+                hDelOmega.fill(atIP_KF.getOmega() - atIP_GBL.getOmega());
+                hDelTanL.fill(atIP_KF.getTanLambda() - atIP_GBL.getTanLambda());
+                hDelD0.fill(atIP_KF.getD0() - atIP_GBL.getD0());
+                hDelZ0.fill(atIP_KF.getZ0() - atIP_GBL.getZ0());
+                if (_debug) {                    
+                    System.out.println("GBL Track state: " + atIP_GBL.toString());
+                    System.out.println("KF Track state:  " + atIP_KF.toString());
+                    System.out.format("Track %d has %d set of GBL residuals:\n", tracks.indexOf(trk), gblResids.size());
+                }
+                for (TrackResidualsData resData : gblResids) {
+                    int vol = resData.getIntVal(resData.getNInt()-1);
+                    if (_debug) System.out.format("  GBL residuals for tracker volume %d\n", vol);
+                    for (int i=0; i<resData.getNDouble(); ++i) {
+                        if (i%2 == 0) {
+                            if (vol == 0) hGBLurt.fill(resData.getDoubleVal(i));
+                            else hGBLurb.fill(resData.getDoubleVal(i));
+                            if (_debug) System.out.format("   Layer %d, biased residual   = %12.6e +- %12.6e\n", resData.getIntVal(i), resData.getDoubleVal(i), resData.getFloatVal(i));
+                        } else {
+                            if (vol == 0) hGBLrt.fill(resData.getDoubleVal(i));
+                            else hGBLrb.fill(resData.getDoubleVal(i));
+                            if (_debug) System.out.format("   Layer %d, unbiased residual = %12.6e +- %12.6e\n", resData.getIntVal(i), resData.getDoubleVal(i), resData.getFloatVal(i));
+                        }
+                    }
+                }
+                Set<TrackResidualsData> kfResids = kfResidsRT.allTo(trk);
+                if (_debug) System.out.format("Kalman track %d has %d set of residuals:\n", tracks.indexOf(trk), kfResids.size());
+                for (TrackResidualsData resData : kfResids) {
+                    int vol = resData.getIntVal(resData.getNInt()-1);
+                    if (_debug) System.out.format("  KF residuals for tracker volume %d\n", vol);
+                    for (int i=0; i<resData.getNDouble(); ++i) {
+                        if (vol == 0) hKFurt.fill(resData.getDoubleVal(i));
+                        else hKFurb.fill(resData.getDoubleVal(i));
+                        if (_debug) System.out.format("   Layer %d, unbiased residual = %12.6e +- %12.6e\n", resData.getIntVal(i), resData.getDoubleVal(i), resData.getFloatVal(i));
+                    }
+                }
+            }
+        }
     }
 
     @Override
@@ -208,7 +336,14 @@ protected void startOfData() {
 
     @Override 
     protected void endOfData() {
-
+        if (_analysis) {
+            try {
+                System.out.println("KalmanToGBLDriver: Outputting the histogram plots now.");
+                aida.saveAs("KalmanToGBLDriverPlots.root");
+            } catch (IOException ex) {
+                Logger.getLogger(TrackingReconstructionPlots.class.getName()).log(Level.SEVERE, null, ex);
+            }
+        }
     }
 
     static Comparator<GBLStripClusterData>  arcLComparator = new Comparator<GBLStripClusterData>() {

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

@@ -133,7 +133,12 @@ public FieldMap(String FileName, String type, boolean uniform, double xOffset, d
         offsets.print("field map offsets");
     }
 
-    double [] getField(Vec r) { // Interpolate the 3D field map
+    /**
+     * Interpolate the 3D field map
+     * @param r       3-vector position
+     * @return        array of 3 field components
+     */
+    double [] getField(Vec r) {
         Vec rHPS;
         if (uniform) {
             rHPS = new Vec(0., 0., 505.57);
@@ -192,6 +197,17 @@ public FieldMap(String FileName, String type, boolean uniform, double xOffset, d
         return Bout; 
     }
 
+    /**
+     * Trilinear interpolation
+     * @param i     x index into array
+     * @param j     y index into array
+     * @param k     z index into array
+     * @param xd    x value
+     * @param yd    y value
+     * @param zd    z value
+     * @param f     3D array to interpolate
+     * @return      interpolated value of the array f at the given coordinates
+     */
     private double triLinear(int i, int j, int k, double xd, double yd, double zd, double[][][] f) {
         // System.out.format(" triLinear: xd=%10.7f, yd=%10.7f, zd=%10.7f\n", xd,yd,zd);
         // System.out.format(" 000=%12.4e, 100=%12.4e\n", f[i][j][k],f[i+1][j][k]);
@@ -210,7 +226,11 @@ private double triLinear(int i, int j, int k, double xd, double yd, double zd, d
         return c;
     }
 
-    void writeBinaryFile(String fName) { // Make a binary field map file that can be read much more quickly
+    /**
+     * Make a binary field map file that can be read much more quickly
+     * @param fName    Name of output file
+     */
+    void writeBinaryFile(String fName) { 
         FileOutputStream ofile;
         try {
             ofile = new FileOutputStream(fName);