DC v2: includes TFLight fix and new dc beta time walk with constants from /calibration/dc/v2

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java

@@ -708,7 +708,9 @@ private boolean filter(int k, boolean forward, double annealingFactor) {
             return false;
         }
     }
-
+
+    // Since no vertex inforamtion, the starting point for path length is the final point at the last layer.
+    // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken.
     private boolean filter(int k, boolean forward) {
         StateVec sVec = sv.transported(forward).get(k);
         org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec mVec = mv.measurements.get(k);
@@ -857,7 +859,9 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) {
     private void calcFinalChisq(int sector) {
         calcFinalChisq(sector, false);
     }
-
+
+    // Since no vertex inforamtion, the starting point for path length is the final point at the last layer.
+    // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken.
     private void calcFinalChisq(int sector, boolean nofilter) {
         int k = svzLength - 1;
         this.chi2 = 0;
@@ -880,9 +884,9 @@ private void calcFinalChisq(int sector, boolean nofilter) {
             sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward);
 
             StateVec svc = sv.transported(forward).get(0);
-            path += svc.deltaPath;
+            path += (forward ? 1 : -1) * svc.deltaPath;
             svc.setPathLength(path);
-
+            
             double V0 = mv.measurements.get(0).surface.unc[0];
 
             Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z());
@@ -922,7 +926,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
 
                 double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
                 svc = sv.transported(forward).get(k1 + 1);
-                path += svc.deltaPath;
+                path += (forward ? 1 : -1) * svc.deltaPath;
                 svc.setPathLength(path);
                 svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x());
                 svc.setProjectorDoca(h);
@@ -975,7 +979,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
             sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward);
 
             StateVec svc = sv.transported(forward).get(0);
-            path += svc.deltaPath;
+            path += (forward ? 1 : -1) * svc.deltaPath;
             svc.setPathLength(path);
 
             Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z());
@@ -1047,7 +1051,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
                 }
 
                 svc = sv.transported(forward).get(k1 + 1);
-                path += svc.deltaPath;
+                path += (forward ? 1 : -1) * svc.deltaPath;
                 svc.setPathLength(path);
 
                 point = new Point3D(sv.transported(forward).get(k1 + 1).x, sv.transported(forward).get(k1 + 1).y, mv.measurements.get(k1 + 1).surface.measPoint.z());
@@ -1116,6 +1120,10 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) {
     public Matrix propagateToVtx(int sector, double Zf) {
         return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer());
     }
+
+    public double getDeltaPathToVtx(int sector, double Zf) {
+        return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer());
+    }
 
     //Todo: apply the common funciton to replace current function above
     @Override

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java

@@ -406,6 +406,8 @@ private void calcFinalChisq(int sector) {
 		calcFinalChisq(sector, false);
 	}
 
+        // Since no vertex inforamtion, the starting point for path length is the final point at the last layer.
+        // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken.
 	private void calcFinalChisq(int sector, boolean nofilter) {
         int k = svzLength - 1;
         this.chi2 = 0;
@@ -426,11 +428,11 @@ private void calcFinalChisq(int sector, boolean nofilter) {
         kfStateVecsAlongTrajectory = new ArrayList<>();
         if (sVec != null && sVec.CM != null) {
 
-        	boolean forward = false;
+           boolean forward = false;
             sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward);            
 
             StateVec svc = sv.transported(forward).get(0);             
-            path += svc.deltaPath;            
+            path += (forward ? 1 : -1) * svc.deltaPath;           
             svc.setPathLength(path);            
 
             double V0 = mv.measurements.get(0).surface.unc[0];
@@ -473,7 +475,7 @@ private void calcFinalChisq(int sector, boolean nofilter) {
 
                 double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]);
                 svc = sv.transported(forward).get(k1+1);
-                path += svc.deltaPath;
+                path += (forward ? 1 : -1) * svc.deltaPath;
                 svc.setPathLength(path);
                 svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x());
                 svc.setProjectorDoca(h);
@@ -504,6 +506,10 @@ public Matrix propagateToVtx(int sector, double Zf) {
         return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer());
     }
 
+    public double getDeltaPathToVtx(int sector, double Zf) {
+        return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer());
+    }
+
     @Override
     public void runFitter(AStateVecs sv, AMeasVecs mv) {
         throw new UnsupportedOperationException("Not supported yet.");

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java

@@ -53,7 +53,7 @@ public void initFromHB(StateVec initSV, double beta) {
         this.trackTrajS.clear();
         this.trackTrajT.put(0, new StateVec(initSV));
     }
-
+    
     /**
      *
      * @param sector
@@ -150,6 +150,104 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, AMeasVecs m
 
         return fVec.CM;
 
+    } 
+
+    /**
+     *
+     * @param sector
+     * @param i initial state vector index
+     * @param Zf
+     * @param iVec state vector at the initial index
+     * @param mv measurements
+     */
+    public double getDeltaPath(int sector, int i, double Zf, StateVec iVec, AMeasVecs mv, Swim swimmer) { // s = signed step-size
+
+        double stepSize = 1.0;
+        StateVec fVec = new StateVec(0);
+        fVec.x = iVec.x;
+        fVec.y = iVec.y;
+        fVec.z = iVec.z;
+        fVec.tx = iVec.tx;
+        fVec.ty = iVec.ty;
+        fVec.Q = iVec.Q;
+        fVec.B = iVec.B;
+        Matrix5x5.copy(iVec.CM, fVec.CM);
+
+        double s = 0;
+        double zInit = mv.measurements.get(i).surface.measPoint.z();
+        double BatMeas = iVec.B;
+
+        double z = zInit;
+
+        while (Math.signum(Zf - zInit) * z < Math.signum(Zf - zInit) * Zf) {
+            z = fVec.z;
+            if (z == Zf) {
+                break;
+            }
+
+            double x = fVec.x;
+            double y = fVec.y;
+            double tx = fVec.tx;
+            double ty = fVec.ty;
+            double Q = fVec.Q;
+            double dPath = fVec.deltaPath;
+            Matrix cMat = new Matrix();
+            Matrix5x5.copy(fVec.CM, cMat);
+            s = Math.signum(Zf - zInit) * stepSize;
+
+            // LOGGER.log(Level.FINE, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s);
+            if (Math.signum(Zf - zInit) * (z + s) > Math.signum(Zf - zInit) * Zf) {
+                s = Math.signum(Zf - zInit) * Math.abs(Zf - z);
+            }
+
+            //rk.RK4transport(sector, Q, x, y, z, tx, ty, s, swimmer, cMat, fVec, dPath);
+            rk.RK4transport(sector, s, swimmer, cMat, fVec);
+
+            // Q  process noise matrix estimate            
+            double p = Math.abs(1. / iVec.Q);
+
+            double X0 = this.getX0(mv.measurements.get(i).surface, z, Z);
+            double t_ov_X0 = Math.abs(s) / X0;//path length in radiation length units = t/X0 [true path length/ X0] ; Ar radiation length = 14 cm
+
+            double beta = this.beta;
+            if (beta > 1.0 || beta <= 0) {
+                beta = 1.0;
+            }
+
+            double sctRMS = 0;
+
+            ////// Todo: Modify multi-scattering or remove it; After update, some parameters, like iteration termintion chonditions, may need to be updated.
+            // Speed of light should be 1
+            // From one measurement site to another, F and Q should be calculated separaetely with multiple steps; and then C' = FTCF + Q
+            if (Math.abs(s) > 0) {
+                sctRMS = ((0.0136) / (beta * PhysicsConstants.speedOfLight() * p)) * Math.sqrt(t_ov_X0)
+                        * (1 + 0.038 * Math.log(t_ov_X0));
+            }
+
+            double cov_txtx = (1 + tx * tx) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS;
+            double cov_tyty = (1 + ty * ty) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS;
+            double cov_txty = tx * ty * (1 + tx * tx + ty * ty) * sctRMS * sctRMS;
+
+            fMS.set(
+                    0, 0, 0, 0, 0,
+                    0, 0, 0, 0, 0,
+                    0, 0, cov_txtx, cov_txty, 0,
+                    0, 0, cov_txty, cov_tyty, 0,
+                    0, 0, 0, 0, 0
+            );
+
+            Matrix5x5.copy(fVec.CM, copyMatrix);
+            Matrix5x5.add(copyMatrix, fMS, fVec.CM);
+
+            if (Math.abs(fVec.B - BatMeas) < 0.0001) {
+                stepSize *= 2;
+            }
+
+            BatMeas = fVec.B;
+        }
+
+        return fVec.deltaPath;
+
     }
 
     /**

reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java

@@ -94,7 +94,7 @@ public static Constants getInstance() {
     public  int     SECTORSELECT = 0;
     public  int     NSUPERLAYERTRACKING = 5;
     private boolean USETSTART = true;
-    private boolean USETIMETBETA = false;
+    //private boolean USETIMETBETA = false;
     private boolean CHECKBETA = false;
     private int     T2D = 1;     // 1=polynomial, 0=exponential
     private boolean USEDOUBLETS = false;
@@ -105,10 +105,10 @@ public static Constants getInstance() {
     public static final String TT             = "/daq/tt/dc";
     public static final String DOCARES        = "/calibration/dc/signal_generation/doca_resolution";
     public static final String TIME2DIST      = "/calibration/dc/time_to_distance/time2dist";
-    public static final String T2DPRESSURE    = "/calibration/dc/time_to_distance/t2d_pressure";
+    public static final String T2DPRESSURE    = "/calibration/dc/v2/t2d_pressure";
     public static final String PRESSURE       = "/hall/weather/pressure";
-    public static final String T2DPRESSUREREF = "/calibration/dc/time_to_distance/ref_pressure";
-    public static final String T0CORRECTION   = "/calibration/dc/time_corrections/T0Corrections";
+    public static final String T2DPRESSUREREF = "/calibration/dc/v2/ref_pressure";
+    public static final String T0CORRECTION   = "/calibration/dc/v2/t0";
     public static final String TDCTCUTS       = "/calibration/dc/time_corrections/tdctimingcuts";
     public static final String WIRESTAT       = "/calibration/dc/tracking/wire_status";
     public static final String TIMEJITTER     = "/calibration/dc/time_jitter";
@@ -296,13 +296,13 @@ public void setUSETSTART(boolean usetstart) {
         USETSTART = usetstart;
     }
 
-    public boolean useUSETIMETBETA() {
-        return USETIMETBETA;
-    }
-
-    public void setUSETIMETBETA(boolean usetimebeta) {
-        USETIMETBETA = usetimebeta;
-    }
+//    public boolean useUSETIMETBETA() {
+//        return USETIMETBETA;
+//    }
+//
+//    public void setUSETIMETBETA(boolean usetimebeta) {
+//        USETIMETBETA = usetimebeta;
+//    }
 
     public double getWIREDIST() {
         return WIREDIST;
@@ -429,7 +429,7 @@ public void printConfig(String engine) {
         LOGGER.log(Level.INFO, "["+engine+"] run with wire ministagger = "  + MINISTAGGERSTATUS.getName());
         LOGGER.log(Level.INFO, "["+engine+"] run with wire feedthroughs = " + FEEDTHROUGHSSTATUS.getName());
         LOGGER.log(Level.INFO, "["+engine+"] run with wire distortions = " + ENDPLATESBOWING);
-        LOGGER.log(Level.INFO, "["+engine+"] run with with time Beta correction (is false for doca Beta correction) = " + USETIMETBETA);
+        //LOGGER.log(Level.INFO, "["+engine+"] run with with time Beta correction (is false for doca Beta correction) = " + USETIMETBETA);
         LOGGER.log(Level.INFO, "["+engine+"] run with with Beta cut = " + CHECKBETA);
         LOGGER.log(Level.INFO, "["+engine+"] run with time to distance function set to exponential/polynomial (0/1) = " + T2D);
         LOGGER.log(Level.INFO, "["+engine+"] run with with hit doublets recovery = " + USEDOUBLETS);

reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java

@@ -791,8 +791,8 @@ private double[] getT0(int sector, int superlayer,
         int cable = this.getCableID1to6(layer, wire);
         int slot = this.getSlotID1to7(wire);
 
-        double t0  = t0Table.getDoubleValue("T0Correction", sector, superlayer, slot, cable);
-        double t0E = t0Table.getDoubleValue("T0Error", sector, superlayer, slot, cable);
+        double t0  = t0Table.getDoubleValue("t0correction", sector, superlayer, slot, cable);
+        double t0E = t0Table.getDoubleValue("t0error", sector, superlayer, slot, cable);
 
         T0Corr[0] = t0;
         T0Corr[1] = t0E;

reconstruction/dc/src/main/java/org/jlab/rec/dc/hit/FittedHit.java

@@ -385,27 +385,28 @@ public void set_TimeToDistance(DataEvent event, double trkAngle, double B, Index
             double deltatime_beta = 0;
             double deltadoca_beta = 0;
 
-            if(Constants.getInstance().useUSETIMETBETA()==true) {
-                deltatime_beta = calcDeltaTimeBetaTFCN(this.get_Time(), tab, beta);
-            } 
-
-            if(event.hasBank("MC::Particle")==false) {
-                distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), 
-                        this.getCorrectedTime(this.get_Time(), deltatime_beta), 
-                        secIdx, slIdx);
-            } else {
-                distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), 
-                        this.getCorrectedTime(this.get_Time(), 0), 
-                        secIdx, slIdx) ;
-            }
-            //get deltadoca
-            if(Constants.getInstance().useUSETIMETBETA()==false) {
-                deltadoca_beta = calcDeltaDocaBeta(distance, tab, beta);
-            }
-
-            distance -=deltadoca_beta;
-            this.set_DeltaTimeBeta(deltatime_beta);
-            this.set_DeltaDocaBeta(deltadoca_beta);
+//            if(Constants.getInstance().useUSETIMETBETA()==true) {
+//                deltatime_beta = calcDeltaTimeBetaTFCN(this.get_Time(), tab, beta);
+//            } 
+//            
+//            if(event.hasBank("MC::Particle")==false) {
+//                distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), 
+//                        this.getCorrectedTime(this.get_Time(), deltatime_beta), 
+//                        secIdx, slIdx);
+//            } else {
+//                distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), 
+//                        this.getCorrectedTime(this.get_Time(), 0), 
+//                        secIdx, slIdx) ;
+//            }
+//            //get deltadoca
+//            if(Constants.getInstance().useUSETIMETBETA()==false) {
+//                deltadoca_beta = calcDeltaDocaBeta(distance, tab, beta);
+//            }
+//            
+//            distance -=deltadoca_beta;
+            //this.set_DeltaTimeBeta(deltatime_beta);
+            //this.set_DeltaDocaBeta(deltadoca_beta);
+            distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), beta, this.get_Time(), secIdx, slIdx);
 
         }
 

reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/T2DFunctions.java

@@ -321,7 +321,7 @@ public static double polyFcnMac(double x, double alpha, double bfield, double v_
         double time = 0;
         // alpha correction 
         double cos30minusalpha=Math.cos(Math.toRadians(30.-alpha));
-        double dmaxalpha = dmax*cos30minusalpha;
+        double dmaxalpha = dmax*cos30minusalpha; 
         double xhatalpha = x/dmaxalpha;
         //   rcapital is an intermediate parameter
         double rcapital = R*dmax;