fix bugs for sign of d0, calculation of residual, and calculation of position and momentum at fixed radius during CVT track finding

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/HelicalTrackFitter.java

@@ -152,7 +152,7 @@ public FitStatus fit(List<Double> X, List<Double> Y, List<Double> Z, List<Double
         _circlefitpars = _circlefit.getFit();
 
         // get the parameters of the helix representation of the track
-        double fit_dca = _circlefitpars.doca(); //check sign convention
+        double fit_dca = -_circlefitpars.doca(); //sign of d0 in the Helix class is oppsite to here
         double fit_phi_at_dca = _circlefitpars.phi();
         double fit_curvature = _circlefitpars.rho();
         double fit_tandip = _linefitpars.slope();
@@ -184,8 +184,8 @@ public FitStatus fit(List<Double> X, List<Double> Y, List<Double> Z, List<Double
         //covr[5] =  delta_dca.delta_dca;
 
         fit_covmatrix[0][0] = _circlefitpars.cov()[5];
-        fit_covmatrix[1][0] = _circlefitpars.cov()[4];
-        fit_covmatrix[2][0] = _circlefitpars.cov()[2];
+        fit_covmatrix[1][0] = -_circlefitpars.cov()[4];
+        fit_covmatrix[2][0] = -_circlefitpars.cov()[2];
         fit_covmatrix[0][1] = _circlefitpars.cov()[4];
         fit_covmatrix[1][1] = _circlefitpars.cov()[3];
         fit_covmatrix[2][1] = _circlefitpars.cov()[1];

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java

@@ -4,12 +4,14 @@
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
+import org.jlab.geom.prim.Point3D;
 import org.jlab.detector.base.DetectorType;
 import org.jlab.rec.cvt.Constants;
 import org.jlab.rec.cvt.bmt.BMTType;
 import org.jlab.rec.cvt.cross.Cross;
 import org.jlab.rec.cvt.fit.CircleFitter;
 import org.jlab.rec.cvt.fit.CircleFitPars;
+import org.jlab.rec.cvt.trajectory.Helix;
 
 public class TrackSeederXY {
 
@@ -235,19 +237,18 @@ private int countType(List<Cross> cand, DetectorType dt) {
 
         return countsvt;
     }
+
+    // A method to calculate residual for a cross relative to a circle from fit 
     private double calcResi(double rho, double d0, double phi0, double xc, double yc) {
-        double r = Math.sqrt(xc*xc+yc*yc);
-        double par = 1. - ((r * r - d0 * d0) * rho * rho) / (2. * (1. + d0 * Math.abs(rho)));
-        double newPathLength = Math.abs(Math.acos(par) / rho);
-        int charge = (int) Math.signum(rho);
+        Helix helix = new Helix(-d0, phi0, rho, 0, 0, 0, 0); // d0 from circular fitting is opposite to d0 in Helix
+
+        double r = Math.sqrt(xc*xc+yc*yc);        
+        double newPathLength = helix.getPathLengthAtRadius(r);
         double alpha = -newPathLength * rho;
-
-        double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(rho)) 
-                * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0));
-        double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(rho)) 
-                * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0));
+
+        double[] xy = helix.getXYAtRadius(alpha); 
 
-        double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc));
+        double res2 = ((xy[0]-xc)*(xy[0]-xc)+(xy[1]-yc)*(xy[1]-yc));
         return Math.sqrt(res2);
     }
 

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java

@@ -255,61 +255,56 @@ public int getCharge() {
 
         return charge;
     }
-
+
+    // a method to return path length at a fixed radius in xy plane
+    public double getPathLengthAtRadius(double r){
+        double par = (2. * (1. - _dca * _curvature) - (r * r - _dca * _dca) * _curvature * _curvature) / (2. * Math.abs(1. - _dca * _curvature));
+        return Math.abs(Math.acos(par) / _curvature);
+    } 
+
+    // a method to return a point (as a vector) at a given radial position.
+    public double[] getXYAtRadius(double alpha){
+        double x = xb - _dca * Math.sin(_phi_at_dca) + 1./_curvature * (Math.sin(_phi_at_dca) - Math.sin(_phi_at_dca + alpha));
+        double y = yb + _dca * Math.cos(_phi_at_dca) - 1./_curvature * (Math.cos(_phi_at_dca) - Math.cos(_phi_at_dca + alpha));
+        return new double[]{x, y};
+    }     
+
+    // a method to return a point (as a vector) at a given radial position.	
     public Point3D getPointAtRadius(double r) {
-        // a method to return a point (as a vector) at a given radial position.	
-        double d0 = _dca;
-        double omega = _curvature;
-        double charge = this.getCharge();
-        double phi0 = _phi_at_dca;
-        double tandip = _tandip;
-        double z0 = _Z0;
-
         if (Math.abs(_curvature)<1E-5) { // R > 100 m, assume it's straight track
-            double x = -d0 * Math.sin(phi0) + r*Math.cos(phi0)+xb;
-            double y =  d0 * Math.cos(phi0) + r*Math.sin(phi0)+yb;
-            double z =  z0 + r*tandip;
+            double x = -_dca * Math.sin(_phi_at_dca) + r*Math.cos(_phi_at_dca)+xb;
+            double y =  _dca * Math.cos(_phi_at_dca) + r*Math.sin(_phi_at_dca)+yb;
+            double z =  _Z0 + r*_tandip;
             return new Point3D(x, y, z);
         }
 
+        double newPathLength = getPathLengthAtRadius(r);
 
-        double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega)));
-        double newPathLength = Math.abs(Math.acos(par) / omega);
-
-        double alpha = -newPathLength * omega;
-
-        double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(omega)) 
-                * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0))+xb;
-        double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(omega)) 
-                * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0))+yb;
-        double z = z0 + newPathLength * tandip;
+        double alpha = -newPathLength * _curvature;
+
+        double[] xy = getXYAtRadius(alpha);
+        double z = _Z0 + newPathLength * _tandip;
 
-        return new Point3D(x, y, z);
+        return new Point3D(xy[0], xy[1], z);
     }
 
-    public Vector3D getTrackDirectionAtRadius(double r) {
-        // a method to return a point (as a vector) at a given radial position.
-        double d0 = _dca;
-        double omega = _curvature;
-        double phi0 = _phi_at_dca;
-        double tandip = _tandip;
-
+    // a method to return a point (as a vector) at a given radial position. 
+    public Vector3D getTrackDirectionAtRadius(double r) {       
         if (Math.abs(_curvature) <1E-5) { // R > 100 m, assume it's straight track
-            double ux =  Math.cos(phi0);
-            double uy =  Math.sin(phi0);
-            double uz =  tandip;
-            return new Vector3D(ux, uy, uz);
+            double ux =  Math.cos(_phi_at_dca);
+            double uy =  Math.sin(_phi_at_dca);
+            double uz =  _tandip;
+            Vector3D trkDir = new Vector3D(ux, uy, uz);
+            return trkDir.asUnit();
         }
 
-        double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega)));
-        double newPathLength = Math.abs(Math.acos(par) / omega);
-
-        double alpha = newPathLength * omega;
+        double newPathLength = getPathLengthAtRadius(r);
+        double alpha = -newPathLength * _curvature;
 
         double sintheta = Math.abs(sinTheta());
 
-        double ux = Math.cos(-alpha + phi0) * sintheta;
-        double uy = Math.sin(-alpha + phi0) * sintheta;
+        double ux = Math.cos(_phi_at_dca + alpha) * sintheta;
+        double uy = Math.sin(_phi_at_dca + alpha) * sintheta;
         double uz = cosTheta();
 
         Vector3D trkDir = new Vector3D(ux, uy, uz);