lcsim/src/org/lcsim/fit/helicaltrack

diff -u -r1.8 -r1.9
--- HelixUtils.java	30 Jun 2009 00:41:12 -0000	1.8
+++ HelixUtils.java	17 Dec 2009 20:32:29 -0000	1.9
@@ -12,6 +12,7 @@

 import hep.physics.matrix.BasicMatrix;
 import hep.physics.matrix.Matrix;
 

+import hep.physics.matrix.SymmetricMatrix;

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

@@ -33,6 +34,7 @@

     private static double _minslope = 1.e-6;
 	private static double _maxpath  = 2400; // todo, add a more "dynamic" maxpath
 	private static double _epsilon  = 0.01; // value below witch a z vector Or an angle is considered as null

+        private static double _tol = 1.e-4; // Maximum aspect ratio of hit covariance terms for hit to measure x-y coord

     
     /**
      * Creates a new instance of HelixUtils.

@@ -48,8 +50,11 @@

      * @return path length from hit1 to hit2
      */
     public static double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {

+        //  Find the position on the circle for the hits
+        Hep3Vector pos1 = getPositionOnCircle(cfit, hit1);
+        Hep3Vector pos2 = getPositionOnCircle(cfit, hit2);

         //  Return the path length between hit1 and hit2

-        return PathCalc(xc(cfit), yc(cfit), RC(cfit), hit1.x(), hit1.y(), hit2.x(), hit2.y());

+        return PathCalc(xc(cfit), yc(cfit), RC(cfit), pos1.x(), pos1.y(), pos2.x(), pos2.y());

     }
     
     /**

@@ -59,8 +64,10 @@

      * @return path length from the DCA to the hit
      */
     public static double PathLength(CircleFit cfit, HelicalTrackHit hit) {

+        //  Find the position on the circle for the hit
+        Hep3Vector pos = getPositionOnCircle(cfit, hit);

         //  Return the path length from the DCA

-        return PathCalc(xc(cfit), yc(cfit), RC(cfit), x0(cfit), y0(cfit), hit.x(), hit.y());

+        return PathCalc(xc(cfit), yc(cfit), RC(cfit), x0(cfit), y0(cfit), pos.x(), pos.y());

     }
     
     /**

@@ -175,7 +182,48 @@

         //  Return the position as a Hep3Vector
         return new BasicHep3Vector(x, y, z);
     }

-    

+
+    private static Hep3Vector getPositionOnCircle(CircleFit cfit, HelicalTrackHit hit) {
+
+        //  Get hit position and errors
+        Hep3Vector pos = hit.getCorrectedPosition();
+        SymmetricMatrix cov = hit.getCorrectedCovMatrix();
+        double dxdx = cov.diagonal(0);
+        double dydy = cov.diagonal(1);
+        if ((dxdx < _tol * dydy) && (dydy < _tol * dxdx)) return pos;
+
+        //  Get circle parameters
+        double x0 = xc(cfit);
+        double y0 = yc(cfit);
+        double R = Math.abs(RC(cfit));
+
+        //  Get hit coordinates
+        double x = pos.x();
+        double y = pos.y();
+
+
+        if (_tol * dxdx > dydy && Math.abs(y-y0) < R) {
+            double xnew1 = Math.sqrt(R*R - (y-y0)*(y-y0)) + x0;
+            double xnew2 = -xnew1 + 2. * x0;
+            if (Math.abs(xnew1 - x) < Math.abs(xnew2 - x)) {
+                x = xnew1;
+            } else {
+                x = xnew2;
+            }
+        }
+        if (_tol * dydy > dxdx && Math.abs(x-x0) < R) {
+            double ynew1 = Math.sqrt(R*R - (x-x0)*(x-x0)) + y0;
+            double ynew2 = -ynew1 + 2. * y0;
+            if (Math.abs(ynew1 - y) < Math.abs(ynew2 - y)) {
+                y = ynew1;
+            } else {
+                y = ynew2;
+            }
+        }
+
+        return new BasicHep3Vector(x, y, pos.z());
+    }
+

     private static double PathCalc(double xc, double yc, double RC, double x1, double y1, double x2, double y2) {
         //  Find the angle between these points measured wrt the circle center
         double phi1 = Math.atan2(y1 - yc, x1 - xc);