lcsim/src/org/lcsim/util/swim
diff -u -r1.2 -r1.3
--- Helix.java 20 Aug 2005 19:16:20 -0000 1.2
+++ Helix.java 20 Aug 2005 21:11:04 -0000 1.3
@@ -2,13 +2,24 @@
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.abs;
+import static java.lang.Math.atan2;
+import static java.lang.Math.asin;
+import static java.lang.Math.PI;
+
+
/**
* This class represents a helix with its axis aligned along Z.
* 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.2 2005/08/20 19:16:20 tonyj Exp $
+ * @version $Id: Helix.java,v 1.3 2005/08/20 21:11:04 jstrube Exp $
*/
public class Helix implements Trajectory
{
@@ -20,26 +31,26 @@
*/
public Helix(Hep3Vector origin, double radius, double phi, double lambda)
{
- if (Math.abs(lambda) >= Math.PI/2) throw new IllegalArgumentException("lambda="+lambda);
+ if (abs(lambda) >= PI/2) throw new IllegalArgumentException("lambda="+lambda);
this.origin = origin;
this.radius = radius;
this.phi = phi;
// Calculate some useful quantities
- cosLambda = Math.cos(lambda);
- sinLambda = Math.sin(lambda);
- xCenter = origin.x() + radius*Math.sin(phi);
- yCenter = origin.y() - radius*Math.cos(phi);
- phiToCenter = Math.atan2(yCenter,xCenter);
- radiusOfCenter = Math.sqrt(xCenter*xCenter + yCenter*yCenter);
+ cosLambda = cos(lambda);
+ sinLambda = sin(lambda);
+ xCenter = origin.x() + radius*sin(phi);
+ yCenter = origin.y() - radius*cos(phi);
+ phiToCenter = atan2(yCenter,xCenter);
+ radiusOfCenter = sqrt(xCenter*xCenter + yCenter*yCenter);
}
public Hep3Vector getPointAtDistance(double alpha)
{
double darg = alpha*cosLambda/radius - phi;
- double x = xCenter + radius*Math.sin( darg );
- double y = yCenter + radius*Math.cos( darg );
+ double x = xCenter + radius*sin( darg );
+ double y = yCenter + radius*cos( darg );
double z = origin.z() + alpha*sinLambda;
return new BasicHep3Vector(x,y,z);
}
@@ -62,10 +73,104 @@
public double getDistanceToInfiniteCylinder(double r)
{
double darg = r*r/(2.*radius*radiusOfCenter) - radiusOfCenter/(2.*radius) - radius/(2.*radiusOfCenter);
- double diff = Math.asin(darg) + phi - phiToCenter;
+ double diff = asin(darg) + phi - phiToCenter;
return (radius/cosLambda)*diff;
}
+ /**
+ * Swims the helix along its trajectory to the
+ * point of closest approach to the given SpacePoint.
+ * The equation to solve for s is an O(2) Taylor-expansion.
+ * The parameterization is as follows:<br />
+ * p<sub>x</sub> = p<sub>0x</sub> cos(&rho s) - p<sub>0y</sub> sin(&rho s)<br />
+ * p<sub>y</sub> = p<sub>0y</sub> cos(&rho s) + p<sub>0x</sub> sin(&rho s)<br />
+ * p<sub>z</sub> = p<sub>0z</sub><br />
+ * x = x<sub>0</sub> + p<sub>0x</sub>/a sin(&rho s) - p<sub>0y</sub>/a (1-cos(&rho s))<br />
+ * y = y<sub>0</sub> + p<sub>0y</sub>/a sin(&rho s) + p<sub>0x</sub>/a (1-cos(&rho s))<br />
+ * z = z<sub>0</sub> + p<sub>z</sub>/p s<br />
+ * a = -0.299792458 B q, where [B] = T and [q] = e;<br />
+ * &rh0 = a / p;<br />
+ * @param point Point in Space to swim to.
+ * @return the length Parameter s
+ */
+ public double getDOCA(Hep3Vector point) {
+ // Calculate both solutions of the quadratic equation.
+ // Then apply an additional constraint.
+ double cosPhi = cos(phi);
+ double sinPhi = sin(phi);
+ double tanLambda = sinLambda/cosLambda;
+ double px = radius*cosPhi;
+ double py = radius*sinPhi;
+ double pz = radius*tanLambda;
+ 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 factorC = factorCNum/factorCDenom;
+
+ // FIXME need to check for consistency between s and alpha. For now this function also swims.
+
+ // "plus case"
+ double length1 = factorB + sqrt(factorB*factorB + factorC);
+ // "minus case"
+ double length2 = factorB - sqrt(factorB*factorB + factorC);
+
+ // Swim the helix to the two possible positions and check which one is closer
+ Hep3Vector newMomentum = getMomentumAtLength(length1);
+ Hep3Vector newPosition = getPositionAtLength(length1);
+
+ // 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);
+ // Spacial difference Vector at the new location
+ Hep3Vector newXDiff = VecOp.sub(newMomentum, point);
+ double rho = cosLambda/radius;
+ if (pMag-rho*VecOp.dot(newXDiff, pCrossB) > 0 )
+ return newXDiff.magnitude();
+ return VecOp.sub(getPositionAtLength(length2), point).magnitude();
+ }
+
+ // Returns the "momentum" at the length s from the starting point.
+ // This uses the definition in http://www.phys.ufl.edu/~avery/fitting/transport.pdf
+ // FIXME if alpa and s can be transformed into one another, this function can be obsoleted
+ // NOTE to FIXME: x and p are orthogonal in the x-y plane (assuming x0 = (0,0,0))
+ // That should be enough of a requirement to obsolete this function
+ private Hep3Vector getMomentumAtLength(double s) {
+ double p0x = radius*cos(phi);
+ double p0y = radius*sin(phi);
+ double p0z = radius*sinLambda/cosLambda;
+ double rho = cosLambda/radius;
+
+ double px = p0x*cos(rho*s) - p0y*sin(rho*s);
+ double py = p0y*cos(rho*s) + p0x*sin(rho*s);
+ double pz = p0z*cos(rho*s) + p0z*(1-cos(rho*s));
+ return new BasicHep3Vector(px, py, pz);
+ }
+
+ // Returns the point at the length s from the starting point.
+ // FIXME if s and alpha are the same or can be transformed into one another, this is obsolete
+ private Hep3Vector getPositionAtLength(double s) {
+ double p0x = radius*cos(phi);
+ double p0y = radius*sin(phi);
+ double p0z = radius*sinLambda/cosLambda;
+ double rho = cosLambda/radius;
+ double p0Mag = sqrt(p0x*p0x + p0y*p0y + p0z*p0z);
+
+ double x = origin.x() + p0x/rho/p0Mag*sin(rho*s) - p0y/rho/p0Mag*(1-cos(rho*s));
+ double y = origin.y() + p0y/rho/p0Mag*sin(rho*s) + p0x/rho/p0Mag*(1-cos(rho*s));
+ double z = origin.z() + p0z/rho/p0Mag*sin(rho*s) + p0z/p0Mag*(s-sin(rho*s)/rho);
+ return new BasicHep3Vector(x, y, z);
+ }
+
private Hep3Vector origin;
private double xCenter;
private double yCenter;
lcsim/src/org/lcsim/util/swim
diff -u -r1.9 -r1.10
--- HelixSwim.java 18 Aug 2005 06:22:29 -0000 1.9
+++ HelixSwim.java 20 Aug 2005 21:11:04 -0000 1.10
@@ -16,9 +16,10 @@
* This implementation works for charged and neutral particles alike.
*
* @author W.Walkowiak
- * @version $Id: HelixSwim.java,v 1.9 2005/08/18 06:22:29 jstrube Exp $
+ * @version $Id: HelixSwim.java,v 1.10 2005/08/20 21:11:04 jstrube Exp $
+ * @deprecated Please use @link{#HelixSwimmer} and @link{#Helix} instead
*/
-public class HelixSwim {
+@Deprecated public class HelixSwim {
/**
* Create a helix swimmmer.
*
