hps-java/src/main/java/org/lcsim/hps/recon/tracking

TrackUtils.java 1.11 -> 1.12

diff -u -r1.11 -r1.12
--- TrackUtils.java	26 Nov 2012 18:12:01 -0000	1.11
+++ TrackUtils.java	27 Nov 2012 02:17:29 -0000	1.12
@@ -19,12 +19,11 @@

 import org.lcsim.event.Track;
 import org.lcsim.fit.helicaltrack.*;
 import org.lcsim.hps.users.phansson.WTrack;

-import org.lcsim.hps.users.phansson.WTrackUtils;

 
 /**
  * 
  * @author Omar Moreno <[log in to unmask]>

- * @version $Id: TrackUtils.java,v 1.11 2012/11/26 18:12:01 phansson Exp $

+ * @version $Id: TrackUtils.java,v 1.12 2012/11/27 02:17:29 phansson Exp $

  * TODO: Switch to JLab coordinates
  */

@@ -33,8 +32,7 @@

         private boolean _debug = false;
 	private boolean isTrackSet = false;
 	private double[] trackParameters;

-        private WTrackUtils _wutils = new WTrackUtils();
-

 	
 	/**
 	 *

@@ -48,8 +46,7 @@

 	 */
 	public void setDebug(boolean debug){		
 		this._debug = debug;

-                _wutils.setDebug(debug);
-	}

+        }

         
         /**
 	 *

@@ -262,79 +259,29 @@

 		return new BasicHep3Vector(x, y, z);
 	}

-        public Hep3Vector calculateIterativeHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield) {
-            return this.calculateIterativeHelixInterceptXPlane(helfit, strip, bfield,false);
-        }

-        public Hep3Vector calculateIterativeHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield, boolean debug) {
-            if(debug) System.out.printf("%s: calculateIterativeHelixInterceptXPlane for heltrackfit\n%s \n",this.getClass().getSimpleName(),helfit.toString());

+        public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, Hep3Vector unit_vec_normal_to_plane,Hep3Vector point_on_plane, double bfield) {
+            /*
+             * Use code in WTrack to find the iterative solution to the interception
+             */
+            boolean debug = true;

             WTrack wtrack = new WTrack(helfit,bfield,true); // B-field sign is flipped so flip!

-            if(debug) {
-                System.out.printf("%s: created Wtrack with params for heltrackfit\n%s \n",this.getClass().getSimpleName(),wtrack.toString());
-                System.out.printf("%s: with associated heltrackfit \n%s \n",this.getClass().getSimpleName(),wtrack._htf.toString());
-            }

+            if(debug) System.out.printf("getHelixPlaneIntercept:find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, and WTrack \n%s \n",point_on_plane.toString(),unit_vec_normal_to_plane.toString(), bfield,wtrack.toString());
+            Hep3Vector intercept_point =  wtrack.getHelixAndPlaneIntercept(wtrack,point_on_plane, unit_vec_normal_to_plane, new BasicHep3Vector(0,0,1));
+            if(debug) System.out.printf("getHelixPlaneIntercept: found intercept point at %s\n",intercept_point.toString());
+            return intercept_point;
+        }
+        
+        public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield) {

             Hep3Vector point_on_plane = strip.origin();
             Hep3Vector unit_vec_normal_to_plane = VecOp.cross(strip.u(),strip.v());//strip.w();

-            if(debug) System.out.printf("%s: find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, and track pars:\n%s \n",this.getClass().getSimpleName(),point_on_plane.toString(),unit_vec_normal_to_plane.toString(), bfield,helfit.toString());
-            boolean d = _wutils.getDebug();
-            if(debug) _wutils.setDebug(debug);        
-            Hep3Vector intercept_point =  _wutils.getHelixAndPlaneIntercept(wtrack, point_on_plane, unit_vec_normal_to_plane, new BasicHep3Vector(0,0,1));
-            _wutils.setDebug(d); //turn off if needed

+            Hep3Vector intercept_point = getHelixPlaneIntercept(helfit, unit_vec_normal_to_plane, point_on_plane, bfield);

             return intercept_point;
         }

-        public double calculateHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip) {
-            /*
-            * Calculate the interception path length on the helix 
-            * where the track fit and the plane (where the strip is located) intercepts 
-            */
-            boolean _DEBUG = false;
-
-            HelicalTrackFit cfit = helfit;
-            if(cfit==null) {
-               System.out.println(this.getClass().getSimpleName() + ": no helical track fit was provided for interception calculation");
-               System.exit(1);
-            }
-            if(_DEBUG) {
-                System.out.println(this.getClass().getSimpleName() + ": calculateHelixIntercept ---");
-                System.out.println(this.getClass().getSimpleName() + ": Track curvatute " + cfit.curvature() + " phi0 " + cfit.phi0());
-            }
-            //Get the origin of the plane
-            Hep3Vector org = strip.origin();
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": origin of plane " + org.toString());
-
-            //Get the point on the helix at the origin
-            double xint = org.x();
-            double s = HelixUtils.PathToXPlane(cfit, xint,0,0).get(0);
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": xint " + xint + "  s " + s );
-
-            //Calculate that point in space
-            Hep3Vector pos_xint = HelixUtils.PointOnHelix(cfit, s);
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": xint position " + pos_xint.toString());
-
-            //Create a vector that points from the origin to that point
-            Hep3Vector vec_xint = VecOp.sub(pos_xint, org);
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": vec_xint " + vec_xint.toString());
-
-            //Get the v in global coordinates
-            // v tells how the unmeasured coordinate is oriented in the tracking frame
-            Hep3Vector v = strip.v();
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": v " + v.toString());
-
-            // We are basically looking for the point where v and track intercepts in the x-y plane
-            // Assuming that the track travels parallel to the x axis we can calculate the 
-            // the correction to xint by scaling
-            // NEED TO TAKE INTO ACCOUNT THE REAL TRACK TRAJECTORY!!
-            double x_corr = vec_xint.y() * v.x()/v.y();
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": vy " + vec_xint.y() + " vx/vy " + v.x()/v.y() + " -> " +  "x_corr " + x_corr);
-            double xint_new = xint + x_corr;
-            if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": final xint_new " + xint_new);
-
-            return xint_new;

-        }
-

/**

@@ -369,7 +316,7 @@

         /**
         *
         */

-        public double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){

+        public static double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){

             return .5*(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1); 
         }

@@ -377,7 +324,7 @@

         /**
         *
         */

-        public boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor){

+        public static boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor){

             boolean debug = false;
             ITransform3D localToGlobal = sensor.getGeometry().getLocalToGlobal();

@@ -385,8 +332,8 @@

             Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid();
             Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
             if(debug){

-                    System.out.println(this.getClass().getSimpleName() + ": Sensor: " + SvtUtils.getInstance().getDescription(sensor));
-                    System.out.println(this.getClass().getSimpleName() + ": Track Position: " + trackPosition.toString());

+                    System.out.println("sensorContainsTrack:  Sensor: " + SvtUtils.getInstance().getDescription(sensor));
+                    System.out.println("sensorContainsTrack:  Track Position: " + trackPosition.toString());

             }
 
             List<Point3D> vertices = new ArrayList<Point3D>();

@@ -399,8 +346,8 @@

                     //vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
                     vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z()));
                     if(debug){

-                            System.out.println(this.getClass().getSimpleName() + ": Vertex 1 Position: " + vertices.get(0).toString());
-                            //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());

+                            System.out.println("sensorContainsTrack:  Vertex 1 Position: " + vertices.get(0).toString());
+                            //System.out.println("sensorContainsTrack:  Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());

                     }
                 } 
                 else if(vertex.y() > 0 && vertex.x() > 0){

@@ -408,8 +355,8 @@

                     //vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
                     vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z()));
                     if(debug){

-                    System.out.println(this.getClass().getSimpleName() + ": Vertex 2 Position: " + vertices.get(1).toString());
-                    //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());

+                    System.out.println("sensorContainsTrack:  Vertex 2 Position: " + vertices.get(1).toString());
+                    //System.out.println("sensorContainsTrack:  Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());

                     }
                 } 
                 else if(vertex.y() > 0 && vertex.x() < 0){

@@ -417,8 +364,8 @@

                     //vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
                     vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z()));
                     if(debug){

-                    System.out.println(this.getClass().getSimpleName() + ": Vertex 3 Position: " + vertices.get(2).toString());
-                    //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());

+                    System.out.println("sensorContainsTrack:  Vertex 3 Position: " + vertices.get(2).toString());
+                    //System.out.println("sensorContainsTrack:  Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());

                     }
                 }             
                 else if(vertex.y() < 0 && vertex.x() < 0){

@@ -426,16 +373,16 @@

                     //vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
                     vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z()));
                     if(debug){

-                    System.out.println(this.getClass().getSimpleName() + ": Vertex 4 Position: " + vertices.get(3).toString());
-                    //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());

+                    System.out.println("sensorContainsTrack:  Vertex 4 Position: " + vertices.get(3).toString());
+                    //System.out.println("sensorContainsTrack:  Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());

                     }
                 } 
             }

-            double area1 = this.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z()); 
-            double area2 = this.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z()); 
-            double area3 = this.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z()); 
-            double area4 = this.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());

+            double area1 = TrackUtils.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z()); 
+            double area2 = TrackUtils.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z()); 
+            double area3 = TrackUtils.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z()); 
+            double area4 = TrackUtils.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());

 
             if((area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0)) return true;

@@ -444,16 +391,16 @@

 
         public HelicalTrackFit makeHelicalTrackFit(double[] helix_parameters) {

-        SymmetricMatrix cov = new SymmetricMatrix(5);
-        double[] chisq = new double[2];
-        int[] ndf = new int[2];
-        Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>();
-        Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
-        HelicalTrackFit fit = new HelicalTrackFit(helix_parameters, cov, chisq, ndf,smap,msmap);
-        return fit;

+            SymmetricMatrix cov = new SymmetricMatrix(5);
+            double[] chisq = new double[2];
+            int[] ndf = new int[2];
+            Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>();
+            Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
+            HelicalTrackFit fit = new HelicalTrackFit(helix_parameters, cov, chisq, ndf,smap,msmap);
+            return fit;

-        public Map<String,Double> calculateTrackHitResidual(HelicalTrackHit hth,HelicalTrackFit track, boolean includeMS ) {

+        public static Map<String,Double> calculateTrackHitResidual(HelicalTrackHit hth,HelicalTrackFit track, boolean includeMS ) {

             Map<String,Double> residuals = new HashMap<String,Double>();
 
             Map<HelicalTrackHit, MultipleScatter> msmap = track.ScatterMap();

hps-java/src/main/java/org/lcsim/hps/users/phansson

WTrack.java 1.4 -> 1.5

diff -u -r1.4 -r1.5
--- WTrack.java	7 Nov 2012 20:54:13 -0000	1.4
+++ WTrack.java	27 Nov 2012 02:17:29 -0000	1.5
@@ -7,6 +7,8 @@

 import hep.physics.vec.BasicHep3Vector;
 import hep.physics.vec.Hep3Vector;
 import hep.physics.vec.VecOp;

+import java.util.ArrayList;
+import java.util.List;

 import org.lcsim.constants.Constants;
 import org.lcsim.fit.helicaltrack.HelicalTrackFit;

@@ -21,11 +23,13 @@

     public enum PARAM{TEST;}
     private double[] _parameters = new double[7];
     public HelicalTrackFit _htf = null;

-    double _bfield;
-    double _bfield_constant;
-

+    private double _bfield;
+    private double _bfield_constant;

     int _q = 0;

+    static int max_iterations_intercept = 10;
+    static double epsilon_intercept = 1e-4;
+

     public WTrack(double [] params, double bfield, int q) {
         _bfield = bfield;
         _q = q;

@@ -36,8 +40,10 @@

     public WTrack(WTrack trk) {
         _bfield = trk._bfield;
         _q = trk.getCharge();

-        _bfield_constant = -1*Constants.fieldConversion*_bfield*_q;

+        _bfield_constant = trk._bfield_constant;

         _parameters = trk.getParameters();

+        _htf = trk._htf;
+        _debug = trk._debug;

     }
     
     public WTrack(HelicalTrackFit track, double bfield) {

@@ -113,13 +119,16 @@

         return ( new BasicHep3Vector(_parameters[4],_parameters[5],_parameters[6]));
     }

+    public String paramsToString() {
+        String str = "";
+        for(int i=0;i<7;++i) str += _parameters[i] + ", ";
+        return str;
+    }

     public String toString() {

-        String str = "WTrack params [";
-        for(int i=0;i<7;++i) str += _parameters[i] + ", ";
-        str += "]";

+        String str = "WTrack params [" + paramsToString() + "]";

         if(this._htf!=null) {

-            str += "\n WTrack corresponding HelicalTrackFit params:\n";

+            str += "\n with corresponding HelicalTrackFit:\n";

             str += this._htf.toString(); 
         }
        return str;

@@ -127,6 +136,212 @@

+    
+    
+    
+     private Hep3Vector getMomentumOnHelix(WTrack track, double s) {
+        double a = track.getBfieldConstant();
+        Hep3Vector p0 = track.getP0();
+        double rho = a / p0.magnitude();
+        double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s);
+        double py = p0.y()*Math.cos(rho*s) + p0.x()*Math.sin(rho*s);
+        double pz = p0.z(); 
+        return (new BasicHep3Vector(px,py,pz));
+    }
+    
+    private Hep3Vector getPointOnHelix(WTrack track, double s) {
+        double a = track.getBfieldConstant();
+        Hep3Vector p0 = track.getP0();
+        Hep3Vector x0 = track.getX0();
+        double rho = a / p0.magnitude();
+        double x = x0.x() + p0.x()/a*Math.sin(rho*s) - p0.y()/a*(1-Math.cos(rho*s));
+        double y = x0.y() + p0.y()/a*Math.sin(rho*s) + p0.x()/a*(1-Math.cos(rho*s));
+        double z = x0.z() + p0.z()/p0.magnitude()*s;
+        return (new BasicHep3Vector(x,y,z));
+    }
+    
+    private double getPathLengthToPlaneApprox(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+        /*
+         * Find the approximate path length to the point xp 
+         * in arbitrary oriented, constant magnetic field with unit vector h
+         */
+        boolean debug = true;
+        double a = track.getBfieldConstant();
+        Hep3Vector p0 = track.getP0();
+        Hep3Vector x0 = track.getX0();
+        double p = p0.magnitude();
+        double rho = a / p;
+        double A = VecOp.dot(eta,VecOp.cross(p0, h))/p*0.5*rho;
+        double B = VecOp.dot(p0,eta)/p;
+        double C = VecOp.dot(VecOp.sub(x0,xp),eta);
+        double t = B*B-4*A*C;
+        if(t<0) {
+            System.out.println(" getPathLengthToPlaneApprox ERROR t is negative (" + t + ")!" );
+            System.out.println(" p " + p + " rho " + rho + " A " + A + " B " + B + " C " + C);
+            System.out.println(" track params\n" + track.toString());
+            System.out.println(" xp " + xp.toString());
+            System.out.println(" eta " + eta.toString());
+            System.out.println(" h " + h.toString());
+            System.exit(1);
+        }
+        double root1 = (-B + Math.sqrt(t)) /(2*A);
+        double root2 = (-B - Math.sqrt(t)) /(2*A);
+
+        // choose the smallest positive solution
+        // if both negative choose the smallest negative ???
+        //if(root1==0 || root2==0) root=0;
+        double root = Math.abs(root1) <= Math.abs(root2) ? root1 : root2;
+//        else if(Math.signum(root1)>0 && Math.signum(root2)<0) root = root1;
+//        else if(Math.signum(root2)>0 && Math.signum(root1)<0) root = root2;
+//        else if(Math.signum(root1)>0 && Math.signum(root2)>0) root =  root1 > root2 ? root2 : root1;
+//        else if(Math.signum(root1)<0 && Math.signum(root2)<0) root =  root1 < root2 ? root2 : root1;
+//        else {
+//            System.out.println(" I should never get here! (root1=" + root1 + " root2=" + root2+")");
+//            System.exit(1);
+//        }
+        if((debug)) {
+                System.out.println(" getPathLengthToPlaneApprox ");
+                System.out.println(" " + track.toString());
+                System.out.println(" xp " + xp.toString());
+                System.out.println(" eta " + eta.toString());
+                System.out.println(" h " + h.toString());
+                System.out.println(" p " + p + " rho " + rho + " t " + t + " A " + A + " B " + B + " C " + C);
+                System.out.println(" root1 " + root1 + " root2 " + root2 + " -> root " + root);
+        }
+        return root;
+    
+    }
+    
+    
+    private Hep3Vector getPointOnHelix(WTrack track, double s, Hep3Vector h) {
+        /*
+         * Get point on helix at path lenght s 
+         * in arbitrary oriented, constant magnetic field with unit vector h
+         */
+        
+        double a = track.getBfieldConstant();
+        Hep3Vector p0 = track.getP0();
+        double p = p0.magnitude();
+        Hep3Vector x0 = track.getX0();
+        double rho = a / p0.magnitude();
+        double srho = s*rho;
+        Hep3Vector a1 = VecOp.mult(1/a*Math.sin(srho), p0);
+        Hep3Vector a2 = VecOp.mult(1/a*(1-Math.cos(srho)),VecOp.cross(p0,h));
+        Hep3Vector a3 = VecOp.mult(VecOp.dot(p0, h)/p,h);
+        Hep3Vector a4 = VecOp.mult(s-Math.sin(srho)/rho, a3);
+        Hep3Vector x = VecOp.add(x0,a1);
+        x = VecOp.sub(x,a2);
+        x = VecOp.add(x,a4);
+        return x;
+    }
+    
+    private Hep3Vector getMomentumOnHelix(WTrack track, double s, Hep3Vector h) {
+        /*
+         * Get point on helix at path lenght s 
+         * in arbitrary oriented, constant magnetic field with unit vector h
+         */
+        
+        double a = track.getBfieldConstant();
+        Hep3Vector p0 = track.getP0();
+        double rho = a / p0.magnitude();
+        double srho = s*rho;
+        Hep3Vector a1 = VecOp.mult(Math.cos(srho), p0);
+        Hep3Vector a2 = VecOp.cross(p0, VecOp.mult(Math.sin(srho),h));
+        Hep3Vector a3 = VecOp.mult(VecOp.dot(p0,h),VecOp.mult(1-Math.cos(srho),h));
+        Hep3Vector p  = VecOp.sub(a1,a2);
+        p = VecOp.add(p,a3);
+        return p;
+    }
+    
+    
+    private double[] getHelixParametersAtPathLength(WTrack track, double s, Hep3Vector h) {
+        /*
+         * Calculate the exact position of the new helix parameters at path length s
+         * in an arbitrarily oriented, constant magnetic field
+         * 
+         * point xp is the point 
+         * h is a unit vector in the direction of the magnetic field
+         */
+        boolean debug = true;
+        // Find track parameters at that path length
+        Hep3Vector p = getMomentumOnHelix(track, s, h);
+        Hep3Vector x = getPointOnHelix(track,s, h);
+        
+        Hep3Vector p_tmp = getMomentumOnHelix(track, s);
+        Hep3Vector x_tmp = getPointOnHelix(track, s);
+        
+        if((debug)) {
+            System.out.println(" point on helix at s");
+            System.out.println(" p  " + p.toString() + "   p_tmp " + p_tmp.toString());
+            System.out.println(" x  " + x.toString() + "   x_tmp " + x_tmp.toString());
+        }
+        
+        
+        //Create the new parameter array
+        double [] pars = new double[7];
+        pars[0] = p.x();
+        pars[1] = p.y();
+        pars[2] = p.z();
+        pars[3] = track.getParameters()[3]; //E is unchanged
+        pars[4] = x.x();
+        pars[5] = x.y();
+        pars[6] = x.z();
+        return pars;
+    }
+    
+    
+    public Hep3Vector getHelixAndPlaneIntercept(Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+        return this.getHelixAndPlaneIntercept(this, xp, eta, h);
+    }
+    public Hep3Vector getHelixAndPlaneIntercept(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+    
+        /*
+         * Find the interception point between the helix and plane
+         * xp: point on the plane
+         * eta: unit vector of the plane
+         * h: unit vector of magnetic field
+         */
+        boolean debug = true;
+        
+        
+        int iteration = 1;
+        double delta_s = 9999;
+        double s_prev = 999999;
+        //List<WTrack> tracks = new ArrayList<WTrack>();
+        WTrack trk = new WTrack(track);
+        while(iteration<=max_iterations_intercept && Math.abs(s_prev)>epsilon_intercept) {
+            if((debug)) System.out.printf("%s: Iteration %d:  s_pref =%.3f delta_s=%.3f\ntrk params: [%s] \n",this.getClass().getSimpleName(),iteration,s_prev,delta_s,paramsToString());
+                
+            // Start by approximating the path length to the point
+            double s = getPathLengthToPlaneApprox(trk, xp, eta, h);
+            
+            if((debug)) System.out.printf("%s: Iteration %d: path length step s=%.3f\n",this.getClass().getSimpleName(),iteration,s);
+            
+            // Find the track parameters at this point
+            double[] params = getHelixParametersAtPathLength(trk, s, h);
+            // update the track parameters           
+            trk.setTrackParameters(params);
+            
+            if((debug)) System.out.printf("%s: Iteration %d: updated trak params: [%s]\n",this.getClass().getSimpleName(),iteration,trk.paramsToString());
+            
+            //tracks.add(trk);
+            iteration++;
+            delta_s = Math.abs(s-s_prev);
+            s_prev = s;
+            
+            //Save distance between point and estimate
+            //Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
+            
+            
+        }
+        
+        if((debug)) System.out.printf("%s: final s=%.3f and ds=%.3f after %d iterations with track params: [%s]\n",this.getClass().getSimpleName(),s_prev,delta_s,iteration,trk.paramsToString());
+        
+        Hep3Vector point = trk.getX0();
+        return point;
+    
+    }
+

Commit in `hps-java/src/main/java/org/lcsim/hps` on MAIN
`recon/tracking/HelicalTrackHitResidualsDriver.java`	+21	-18	1.1 -> 1.2
`/DataTrackerFakeHitDriver.java`	+7	-30	1.6 -> 1.7
`/TrackUtils.java`	+41	-94	1.11 -> 1.12
`users/phansson/StripMPAlignmentInput.java`	+15	-25	1.1 -> 1.2
`/MPAlignmentInputCalculator.java`		-2	1.1 -> 1.2
`/WTrack.java`	+223	-8	1.4 -> 1.5
`/ModuleMPAlignmentInput.java`	+2	-1	1.1 -> 1.2
	+309	-178