Updates of forward tracking since November of 2022

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

@@ -14,6 +14,7 @@
 /**
  *
  * @author ziegler
+ * @author Tongtong Cao
  */
 
 public abstract class AMeasVecs {
@@ -22,7 +23,8 @@ public abstract class AMeasVecs {
 
     public void setMeasVecs(List<Surface> measSurfaces) {
         measurements = new ArrayList<>();
-        Collections.sort(measSurfaces);
+        if(measSurfaces.get(0).type != Type.LINEDOCA && measSurfaces.get(0).type != Type.PLANEURWELL) // Measurements from URWell and DC has been sorted 
+        	Collections.sort(measSurfaces);
         for(int i = 0; i < measSurfaces.size(); i++) {
             MeasVec mvec = new MeasVec();
             mvec.k = i ;
@@ -33,9 +35,62 @@ public void setMeasVecs(List<Surface> measSurfaces) {
             mvec.skip = mvec.surface.passive;
             mvec.hemisphere = measSurfaces.get(i).hemisphere;
             measurements.add(mvec);
+            if(measSurfaces.get(i).type == Type.LINEDOCA) {
+            	mvec.region = measSurfaces.get(i).region;
+            	mvec.sector = measSurfaces.get(i).getSector();
+            	mvec.superlayer = measSurfaces.get(i).getSuperLayer();
+            	mvec.layer = measSurfaces.get(i).getLayer();
+            }
+            else if(measSurfaces.get(i).type == Type.PLANEURWELL){
+                mvec.region = measSurfaces.get(i).region;
+                mvec.layer = measSurfaces.get(i).getLayer();
+            	mvec.sector = measSurfaces.get(i).getSector();
+            }
         }
     }
 
+    // For DC, there could be two hits at a measurement
+    public double[] dhDoca(int k, StateVec stateVec) {  
+    	double value[] = {Double.NaN, Double.NaN};
+
+    	Surface surf = this.measurements.get(stateVec.k).surface;
+    	Point3D point = new Point3D(stateVec.x, stateVec.y, stateVec.z);
+    	double h = hDoca(point, surf.wireLine[0]);
+
+    	double signMeas = 1;
+    	double sign = 1;
+    	if(surf.doca[1]!=-99 || !(Math.abs(surf.doca[0])<0.5 && surf.doca[1]==-99 ) ) { //use LR only for double hits or large enough doca for one hit
+    		signMeas = Math.signum(surf.doca[0]);
+    		sign = Math.signum(h);
+    	} else {
+    		signMeas = Math.signum(h);
+    		sign = Math.signum(h);
+    	}
+
+    	value[0] = signMeas*Math.abs(surf.doca[0]) - sign*Math.abs(h);
+
+        //USE THE DOUBLE HIT
+        if(surf.doca[1]!=-99) { 
+            h = hDoca(point, surf.wireLine[1]);
+
+            signMeas = Math.signum(surf.doca[1]);
+            sign = Math.signum(h);
+
+            value[1]= signMeas*Math.abs(surf.doca[1]) - sign*Math.abs(h);
+        } 
+
+    	return value;
+    }
+
+    // Return a signed doca for DC
+    public double hDoca(Point3D point, Line3D wireLine) {              
+
+        Line3D WL = new Line3D();
+        WL.copy(wireLine);
+        WL.copy(WL.distance(point));
+
+        return WL.length()*Math.signum(WL.direction().x());
+    }
 
     public double dh(int k, StateVec stateVec) {
 
@@ -209,6 +264,10 @@ public class MeasVec implements Comparable<MeasVec> {
         public boolean skip  = false;
         public double hemisphere = 1;
 
+        //For DC
+        public int region = -1;
+        public int sector = -1;
+        public int superlayer = -1;
 
         @Override
         public int compareTo(MeasVec arg) {

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

@@ -10,6 +10,14 @@
 import org.jlab.clas.tracking.kalmanfilter.helical.KFitter;
 import org.jlab.clas.tracking.trackrep.Helix;
 import org.jlab.geom.prim.Point3D;
+import org.jlab.jnp.matrix.Matrix;
+import org.jlab.jnp.matrix.Matrix5x5;
+
+/**
+ *
+ * @author ziegler
+ * @author Tongtong Cao
+ */
 
 public abstract class AStateVecs {
 
@@ -245,6 +253,79 @@ public class StateVec {
         public double energyLoss;
         public double dx;
         public double path;
+
+        /////////////////////// extra variables for forward tracking ///////////////////////
+        // tx = px/pz
+        public double ty; //=py/pz
+        public double Q; //q/p
+        public double B;
+        public double deltaPath;   
+        public Matrix CM = new Matrix();;
+        private double _PathLength;
+
+        /////////////////////// For DAF ///////////////////////
+        private double weightDAF_single = 1;
+        private double[] weightDAF_double = {0.5, 0.5};
+
+        public double getWeightDAF_singleHit(){
+            return weightDAF_single;
+        }
+        public void setWeightDAF_singleHit(double weight){
+            this.weightDAF_single = weight;
+        }
+
+        public double[] getWeightDAF_doubleHits(){
+            return weightDAF_double;
+        }
+        public void setWeightDAF_doubleHits(double[] weight){
+            this.weightDAF_double = weight;
+        }
+
+        private double finalDAFWeight = -999;
+        public double getFinalDAFWeight(){
+            return finalDAFWeight;
+        }
+        public void setFinalDAFWeight(double weight){
+            this.finalDAFWeight = weight;
+        }
+
+        private int SorDHit = -1; // 1 : single hit; 0: double hit for DC
+        public int getSorDHit(){
+            return SorDHit;
+        }
+        public void setSorDHit(int SorDHit){
+            this.SorDHit = SorDHit;
+        }
+
+        public double getPathLength() {
+            return _PathLength;
+        }
+
+        public void setPathLength(double _PathLength) {
+            this._PathLength = _PathLength;
+        }
+
+        // KF projector --> get Wire midPoint match
+        private double hw;
+
+        public double getProjector() {
+            return hw;
+        }
+
+        public void setProjector(double h) {
+            this.hw = h;
+        }
+        // KF projector --> get fit doca
+        private double h;
+
+        public double getProjectorDoca() {
+            return h;
+        }
+
+        public void setProjectorDoca(double h) {
+            this.h = h;
+        }
+        /////////////////////// extra variables for forward tracking ///////////////////////
 
         public StateVec(int k) {
             this.k = k;
@@ -294,8 +375,18 @@ public final void copy(StateVec s) {
             this.energyLoss = s.energyLoss;
             this.dx = s.dx;
             this.path = s.path;
-            this.covMat = this.copyMatrix(s.covMat);
-            this.F = this.copyMatrix(s.F);
+            if(s.covMat != null)
+            	this.covMat = this.copyMatrix(s.covMat);
+            if(s.F != null)
+            	this.F = this.copyMatrix(s.F);            
+            this.ty = s.ty;
+            this.Q = s.Q;
+            this.B = s.B;
+            this.deltaPath = s.deltaPath;
+            if(s.CM != null)
+            	Matrix5x5.copy(s.CM, this.CM);
+            this.weightDAF_single = s.weightDAF_single;
+            this.weightDAF_double = s.weightDAF_double;
         }
 
         public void copyCovMat(double[][] c) {
@@ -414,8 +505,7 @@ public Helix getHelix(double xref, double yref) {
             if (Math.abs(phi0) > Math.PI) phi0 -= Math.signum(phi0) * 2 * Math.PI;
             double tanDip   = vec.tanL;
             double z0       = vec.z0 + vec.dz;
-            //double omega    = turningSign * KFitter.polarity / R;  
-            double omega    = Math.signum(vec.kappa) / R;  
+            double omega    = -turningSign / R;
             double d0       = -vec.d_rho;
 
             Helix helix = new Helix(d0, phi0, omega, z0, tanDip, turningSign, bfield, xref, yref, units);

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

@@ -14,6 +14,7 @@
 /**
  *
  * @author ziegler
+ * @author Tongtong Cao
  */
 public class Surface implements Comparable<Surface> {
 
@@ -40,6 +41,57 @@ public class Surface implements Comparable<Surface> {
     public boolean passive = false;
     public double hemisphere = 1;
 
+    // For DC
+    public double z;
+    public double x;
+    public double[] unc = new double[2];
+    public double tilt;
+    public double[] doca = new double[2];
+    public double wireMaxSag;
+    public Line3D[] wireLine = new Line3D[2];
+    public boolean reject = false;
+    public int region;
+    public int superlayer;
+    public int nMeas = 1;
+
+    // For URWell
+    public double y;
+    public double x_err;
+    public double y_err;
+
+    public void setNMeas(int n) {
+    	nMeas = n;
+    }
+
+    public int getNMeas() {
+    	return nMeas;
+    }
+
+    // For URWell
+    public Surface(int sector, double x, double y, double z, double x_err, double y_err) {
+        type = Type.PLANEURWELL;
+        this.sector = sector;
+        this.x = x;
+        this.y = y;
+        this.z = z;
+        this.x_err = x_err;
+        this.y_err = y_err;
+    }
+
+    // For DC
+    public Surface(int region, double z, double x, double tilt, double wireMaxSag, double[] doca, double[] unc, double error, Line3D[] wireLine) {
+    	type = Type.LINEDOCA;
+    	this.region = region;
+    	this.z = z;
+        this.x = x;
+        this.tilt = tilt;
+        this.wireMaxSag = wireMaxSag;
+        this.doca = doca;
+        this.unc = unc;
+        this.error = error;
+        this.wireLine = wireLine;
+    }
+
     public Surface(Plane3D plane3d, Point3D refrPoint, Point3D c1, Point3D c2, double accuracy) {
         type = Type.PLANEWITHPOINT;
         plane = plane3d;
@@ -184,6 +236,20 @@ public int getLayer() {
     public void setLayer(int layer) {
         this.layer = layer;
     }
+
+    /**
+     * @return the superlayer
+     */
+    public int getSuperLayer() {
+        return superlayer;
+    }
+
+    /**
+     * @param superlayer the superlayer to set
+     */
+    public void setSuperLayer(int superlayer) {
+        this.superlayer = superlayer;
+    }
 
     /**
      * @return the sector

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

@@ -3,11 +3,12 @@
 /**
  *
  * @author ziegler
+ * @author Tongtong Cao
  */
 public enum Type {
     UDF(-1), PLANEWITHPOINT(0), PLANEWITHLINE(1), PLANEWITHSTRIP(2),
     CYLINDERWITHPOINT(3), CYLINDERWITHLINE(4), CYLINDERWITHARC(5), 
-    CYLINDERWITHSTRIP(6), LINE(7);
+    CYLINDERWITHSTRIP(6), LINE(7), LINEDOCA(8), PLANEURWELL(9);
     private final int value;
 
     Type(int value) {