lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.9 -r1.10
--- Helix.java 12 Jul 2007 22:02:03 -0000 1.9
+++ Helix.java 12 Jul 2007 23:30:29 -0000 1.10
@@ -23,45 +23,46 @@
* <p>
* For more info on swimming see <a href="doc-files/transport.pdf">this paper</a>
* by Paul Avery.
- *
+ *
* @author tonyj
- * @version $Id: Helix.java,v 1.9 2007/07/12 22:02:03 tknelson Exp $
+ * @version $Id: Helix.java,v 1.10 2007/07/12 23:30:29 tknelson Exp $
*/
-public class Helix implements Trajectory {
+public class Helix implements Trajectory
+{
Hep3Vector origin;
-
+
double xCenter;
-
+
double yCenter;
-
+
private double radius;
-
+
private double sinLambda;
-
+
private double cosLambda;
-
+
private double sinPhi;
-
+
private double cosPhi;
-
+
private double phi;
-
+
// parameterization in terms of 'momentum'
// A helix is a mathematical object and doesn't have "momentum",
// but some of the used algorithms are expressed in terms of it.
// That's OK, it's a private variable.
SpaceVector momentum;
-
+
private double abs_r;
-
+
private double rho;
-
+
/**
* Creates a new instance of Helix. Parameters according to <a
* href="doc-files/L3_helix.pdf">the L3 conventions</a><br />
* Please also have a look at <img src="doc-files/Helix-1.png"
* alt="Helix-1.png"> <img src="doc-files/Helix-2.png" alt="Helix-2.png">
- *
+ *
* @param origin
* A point on the helix
* @param radius
@@ -74,14 +75,15 @@
* @param lambda
* The dip angle w/rt positive part of the x-y plane
*/
- public Helix(SpacePoint org, double r, double Phi, double lambda) {
+ public Helix(SpacePoint org, double r, double Phi, double lambda)
+ {
// if (abs(lambda) > PI/2.)
// throw new IllegalArgumentException("lambda = " + lambda + " is
// outside of -PI/2<lambda<PI/2");
origin = org;
radius = r;
phi = Phi;
-
+
// Calculate some useful quantities
cosPhi = cos(phi);
sinPhi = sin(phi);
@@ -91,34 +93,37 @@
yCenter = origin.y() - radius * cosPhi;
setSpatialParameters();
}
-
+
/**
* returns a point on the Helix at a distance alpha from the origin along
* the trajectory. alpha == 2*PI*radius/cosLambda corresponds to one
* rotation in the x-y plane
*/
- public SpacePoint getPointAtLength(double alpha) {
+ public SpacePoint getPointAtLength(double alpha)
+ {
double angle = phi + alpha * rho;
double x = xCenter - radius * sin(angle);
double y = yCenter + radius * cos(angle);
double z = origin.z() + alpha * sinLambda;
return new CartesianPoint(x, y, z);
}
-
- public SpacePoint getCenterXY() {
+
+ public SpacePoint getCenterXY()
+ {
return new CartesianPoint(xCenter, yCenter, 0);
}
-
+
/**
* Calculates the distance along the Helix until it hits a cylinder centered
* along z
- *
+ *
* @param r
* the radius of the cylinder
* @return the distance along the trajectory or Double.NAN if the cylinder
* can never be hit
*/
- public double getDistanceToInfiniteCylinder(double r) {
+ public double getDistanceToInfiniteCylinder(double r)
+ {
double phiToCenter = atan2(yCenter, xCenter);
double radiusOfCenter = sqrt(xCenter * xCenter + yCenter * yCenter);
// Negative radius doesn't make sense
@@ -137,30 +142,33 @@
result += Math.abs(radius / cosLambda) * 2 * Math.PI;
return result;
}
-
- public double getDistanceToZPlane(double z) {
+
+ public double getDistanceToZPlane(double z)
+ {
return (z - origin.z()) / sinLambda;
}
-
+
/**
* @param alpha
* The distance along the trajectory in the x-y plane
* @return The unit vector of the momentum
*/
- public SpaceVector getUnitTangentAtLength(double alpha) {
+ public SpaceVector getUnitTangentAtLength(double alpha)
+ {
double angle = phi + alpha * rho;
return new CartesianVector(cos(angle), sin(angle), sinLambda
/ cosLambda);
}
-
- public double getRadius() {
+
+ public double getRadius()
+ {
return radius;
}
-
+
/**
* Calculates the <em>signed</em> distance in mm between the Helix and an
* arbitrary point in Space
- *
+ *
* @param point
* the point in space to calculate the distance to
* @return the distance in mm between the point and the helix at the point
@@ -168,24 +176,24 @@
*/
// FIXME The translation in Z should be replaced, because it can lead to an
// infinite loop and the formula is only valid in 2D
- public double getSignedClosestDifferenceToPoint(SpacePoint point) {
+ public double getSignedClosestDifferenceToPoint(SpacePoint point)
+ {
double tanLambda = sinLambda / cosLambda;
Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
Hep3Vector xDiff = VecOp.sub(origin, point);
double pMag = momentum.magnitude();
-
+
// translate along Z because algorithm can handle only numbers in the
// first quadrant
-
+
// modified to get rid of while loops by Tim Nelson and Cosmin Deaconu
double zPos = xDiff.z();
- if (zPos > 0)
- {
- zPos = Math.IEEEremainder(zPos,abs(radius * tanLambda * Math.PI / 4));
- }
- else
+ zPos = Math.IEEEremainder(zPos, abs(radius * tanLambda * Math.PI / 4));
+ if (zPos < 0) zPos += abs(radius * tanLambda * Math.PI / 4);
+
+ if (zPos/abs(radius * tanLambda * Math.PI / 4) < 0 || zPos/abs(radius * tanLambda * Math.PI / 4) > 1)
{
- zPos = Math.abs(radius * tanLambda * Math.PI / 4) - Math.IEEEremainder(zPos,abs(radius * tanLambda * Math.PI / 4));
+ System.out.println("Valid range of zPos/abs(radius*tanLambda*Math.PI/4) [0 and 1], value is: "+zPos/abs(radius * tanLambda * Math.PI / 4));
}
// these are two mutually exclusive cases and two while loops
@@ -197,7 +205,7 @@
// zPos += abs(radius * tanLambda * Math.PI / 4);
// }
xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
-
+
double numerator = (-2 * VecOp.dot(xDiff, pCrossB) + pMag * rho
* (xDiff.magnitudeSquared() - xDiff.z() * xDiff.z()))
/ radius;
@@ -208,88 +216,95 @@
/ radius);
return numerator / denominator;
}
-
- // added by Nick Sinev
- public double getSecondDistanceToInfiniteCylinder(double r) {
- double result = getDistanceToInfiniteCylinder(r);
- SpacePoint first = getPointAtLength(result);
- double angto0 = Math.atan2(-yCenter, -xCenter);
- double angtofirst = Math
- .atan2(first.y() - yCenter, first.x() - xCenter);
- double dang = angtofirst - angto0;
- if (dang < -Math.PI)
- dang += 2. * Math.PI;
- if (dang > Math.PI)
- dang -= 2. * Math.PI;
- double angofarc = 2. * (Math.PI - Math.abs(dang));
- double arc = Math.abs(radius) * angofarc / cosLambda;
- return result + arc;
- }
-
-
- public double getZPeriod() {
- return Math.PI * radius * 2. * sinLambda / cosLambda;
- }
-
- /**
- * Swims the helix along its trajectory to the point of closest approach to
- * the given SpacePoint. For more info on swimming see <a
- * href=doc-files/fitting/transport.pdf> Paul Avery's excellent text</a>
- *
- * @param point
- * Point in Space to swim to.
- * @return the length along the trajectory
- */
- public double getTrackLengthToPoint(SpacePoint point) {
- double tanLambda = sinLambda / cosLambda;
- double pMag = momentum.magnitude();
- Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
-
- Hep3Vector xDiff = VecOp.sub(origin, point);
- int addedQuarterPeriods = 0;
- if (abs(tanLambda) > 1e-10) {
- double zPos = xDiff.z();
- while (abs(zPos) > abs(radius * tanLambda * Math.PI)) {
- zPos -= signum(zPos) * abs(radius * tanLambda * Math.PI);
- ++addedQuarterPeriods;
- }
- // Make sure the helix is in the right quadrant for the atan
- if (zPos > 0 && addedQuarterPeriods > 0)
- addedQuarterPeriods *= -1;
- if (addedQuarterPeriods % 2 != 0)
- addedQuarterPeriods += signum(addedQuarterPeriods);
- xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
- }// int addedHalfPeriods = 0;
- // while (abs(xDiff.z()) > abs(tanLambda * radius * 0.5*Math.PI)) {
- // addedHalfPeriods++;
- // xDiff = VecOp.sub(getPointAtLength(addedHalfPeriods*0.5*Math.PI/rho),
- // point);
- // }
- double pz = momentum.z();
- double factorA1 = pMag - pz * pz / pMag - (VecOp.dot(xDiff, pCrossB))
- * rho;
- double factorA2 = (VecOp.dot(xDiff, momentum) - xDiff.z() * pz) * rho;
- return addedQuarterPeriods * abs(radius / cosLambda * Math.PI)
- - Math.atan2(factorA2, factorA1) / rho;
- }
-
- /**
- * Sets the parametrization in terms of "momentum" and charge
- *
- */
- private void setSpatialParameters() {
- abs_r = abs(radius);
- momentum = new CartesianVector(abs_r * cosPhi, abs_r * sinPhi, abs_r
- * sinLambda / cosLambda);
- rho = -cosLambda / radius;
- }
-
- public String toString() {
- String helix = String.format("Helix:\n");
- helix += String.format("radius: %f\n", radius);
- helix += String.format("phi: %f\n", phi);
- helix += String.format("lambda: %f\n", acos(cosLambda));
- helix += String.format("rho: %f\n", rho);
- return helix;
+
+ // added by Nick Sinev
+ public double getSecondDistanceToInfiniteCylinder(double r)
+ {
+ double result = getDistanceToInfiniteCylinder(r);
+ SpacePoint first = getPointAtLength(result);
+ double angto0 = Math.atan2(-yCenter, -xCenter);
+ double angtofirst = Math
+ .atan2(first.y() - yCenter, first.x() - xCenter);
+ double dang = angtofirst - angto0;
+ if (dang < -Math.PI)
+ dang += 2. * Math.PI;
+ if (dang > Math.PI)
+ dang -= 2. * Math.PI;
+ double angofarc = 2. * (Math.PI - Math.abs(dang));
+ double arc = Math.abs(radius) * angofarc / cosLambda;
+ return result + arc;
+ }
+
+
+ public double getZPeriod()
+ {
+ return Math.PI * radius * 2. * sinLambda / cosLambda;
+ }
+
+ /**
+ * Swims the helix along its trajectory to the point of closest approach to
+ * the given SpacePoint. For more info on swimming see <a
+ * href=doc-files/fitting/transport.pdf> Paul Avery's excellent text</a>
+ *
+ * @param point
+ * Point in Space to swim to.
+ * @return the length along the trajectory
+ */
+ public double getTrackLengthToPoint(SpacePoint point)
+ {
+ double tanLambda = sinLambda / cosLambda;
+ double pMag = momentum.magnitude();
+ Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
+
+ Hep3Vector xDiff = VecOp.sub(origin, point);
+ int addedQuarterPeriods = 0;
+ if (abs(tanLambda) > 1e-10)
+ {
+ double zPos = xDiff.z();
+ while (abs(zPos) > abs(radius * tanLambda * Math.PI))
+ {
+ zPos -= signum(zPos) * abs(radius * tanLambda * Math.PI);
+ ++addedQuarterPeriods;
+ }
+ // Make sure the helix is in the right quadrant for the atan
+ if (zPos > 0 && addedQuarterPeriods > 0)
+ addedQuarterPeriods *= -1;
+ if (addedQuarterPeriods % 2 != 0)
+ addedQuarterPeriods += signum(addedQuarterPeriods);
+ xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
+ }// int addedHalfPeriods = 0;
+ // while (abs(xDiff.z()) > abs(tanLambda * radius * 0.5*Math.PI)) {
+ // addedHalfPeriods++;
+ // xDiff = VecOp.sub(getPointAtLength(addedHalfPeriods*0.5*Math.PI/rho),
+ // point);
+ // }
+ double pz = momentum.z();
+ double factorA1 = pMag - pz * pz / pMag - (VecOp.dot(xDiff, pCrossB))
+ * rho;
+ double factorA2 = (VecOp.dot(xDiff, momentum) - xDiff.z() * pz) * rho;
+ return addedQuarterPeriods * abs(radius / cosLambda * Math.PI)
+ - Math.atan2(factorA2, factorA1) / rho;
+ }
+
+ /**
+ * Sets the parametrization in terms of "momentum" and charge
+ *
+ */
+ private void setSpatialParameters()
+ {
+ abs_r = abs(radius);
+ momentum = new CartesianVector(abs_r * cosPhi, abs_r * sinPhi, abs_r
+ * sinLambda / cosLambda);
+ rho = -cosLambda / radius;
+ }
+
+ public String toString()
+ {
+ String helix = String.format("Helix:\n");
+ helix += String.format("radius: %f\n", radius);
+ helix += String.format("phi: %f\n", phi);
+ helix += String.format("lambda: %f\n", acos(cosLambda));
+ helix += String.format("rho: %f\n", rho);
+ return helix;
+ }
}
-}