diff -u -r1.8 -r1.9
--- Helix.java	23 Aug 2005 22:06:02 -0000	1.8
+++ Helix.java	24 Aug 2005 18:15:01 -0000	1.9
@@ -19,7 +19,7 @@

  * All quantities in this class are dimensionless. It has no dependencies
  * except for Hep3Vector (which could easily be removed).
  * @author tonyj

- * @version $Id: Helix.java,v 1.8 2005/08/23 22:06:02 tonyj Exp $

+ * @version $Id: Helix.java,v 1.9 2005/08/24 18:15:01 jstrube Exp $

  */
 public class Helix implements Trajectory
 {

@@ -89,25 +89,29 @@

 	   double px = radius*cosPhi;
 	   double py = radius*sinPhi;
 	   double pz = radius*tanLambda;

+	   double rho = -cosLambda/radius;
+       // initial momentum
+	   Hep3Vector p0 = new BasicHep3Vector(px, py, pz);
+       Hep3Vector pCrossB = new BasicHep3Vector(py, -px, 0);

 	   Hep3Vector xDiff = VecOp.sub(origin, point);
 	   
 	   double pMag = sqrt(px*px + py*py + pz*pz);

-	   double factorA = - xDiff.x()*sinPhi + xDiff.y()*cosPhi;  
-	   
-	   double numerator = (pMag + cosLambda*factorA)*radius;
-	   double denominator = cosLambda*cosLambda* factorA;
-	          
-	   double factorB = numerator/denominator;
-	   double factorCNum = 2*radius*radius*(xDiff.x()*cosPhi + xDiff.y()*sinPhi + xDiff.z()*tanLambda);
-	   double factorCDenom = cosLambda*cosLambda*(xDiff.x()*cosPhi + xDiff.y()*sinPhi);

+	   double factorA = VecOp.dot(xDiff, p0);  

-	   double factorC = factorCNum/factorCDenom;

+       double factorB = pMag - rho* VecOp.dot(xDiff, pCrossB);
+	   double factorC = xDiff.z()*p0.z();
+	   double factorD = .5 * rho * rho;
+       
+       double denominator = (factorA - factorC)*factorD;
+       double factorP = factorB / denominator;
+       double factorQ = factorA / denominator;

 	   
 	   // "plus case"

-	   double length1 = factorB + sqrt(factorB*factorB + factorC);

+	   double length1 = factorP/2 + sqrt(factorP/2*factorP/2 + factorQ);

 	   // "minus case"

-	   double length2 = factorB - sqrt(factorB*factorB + factorC);

+	   double length2 = factorP/2 - sqrt(factorP/2*factorP/2 + factorQ);

+       System.err.printf("Found length1: %.3f\tFound length2: %.3f\n", length1, length2);

 	   // Swim the helix to the two possible positions and check which one is closer
 	   Hep3Vector newMomentum = getMomentumAtDistance(length1);
 	   Hep3Vector newPosition = getPointAtDistance(length1);

@@ -115,11 +119,11 @@

 	   // see if the sufficient requirement of a minimum is met
 	   
 	   // Momentum at the new position x B Field;

-	   Hep3Vector pCrossB = new BasicHep3Vector(newMomentum.y(), -newMomentum.x(), 0);

+	   Hep3Vector newPCrossB = new BasicHep3Vector(newMomentum.y(), -newMomentum.x(), 0);

 	   // Spacial difference Vector at the new location
 	   Hep3Vector newXDiff = VecOp.sub(newPosition, point);

-	   double rho = -cosLambda/radius;
-	   if (pMag-rho*VecOp.dot(newXDiff, pCrossB) > 0 )

+       double secondDerivative = pMag+rho*VecOp.dot(newXDiff, newPCrossB);
+	   if (secondDerivative > 0 )

 		   return newXDiff.magnitude();
 	   return VecOp.sub(getPointAtDistance(length2), point).magnitude();
    }