lcsim/src/org/lcsim/util/swim
diff -u -r1.9 -r1.10
--- Helix.java 24 Aug 2005 18:15:01 -0000 1.9
+++ Helix.java 26 Aug 2005 08:23:26 -0000 1.10
@@ -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.9 2005/08/24 18:15:01 jstrube Exp $
+ * @version $Id: Helix.java,v 1.10 2005/08/26 08:23:26 jstrube Exp $
*/
public class Helix implements Trajectory
{
@@ -60,13 +60,19 @@
return (z - origin.z())/sinLambda;
}
+ /**
+ * 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)
{
+ // Negative radius doesn't make sense
if (r<0) throw new IllegalArgumentException("radius "+r+"<0");
double darg = r*r/(2.*radius*radiusOfCenter) - radiusOfCenter/(2.*radius) - radius/(2.*radiusOfCenter);
double diff = asin(darg) + phi - phiToCenter;
double result = (radius/cosLambda)*diff;
- if (result < 0) result += Math.abs(radius/cosLambda)*2*Math.PI;
+ while (result < 0) result += Math.abs(radius/cosLambda)*2*Math.PI;
return result;
}
@@ -77,62 +83,59 @@
* @param point Point in Space to swim to.
* @return the length Parameter s
*/
+ // FIXME works only in the first period apparently. Fix: translate along z into first period.
public double getDistanceFromHelixToPoint(Hep3Vector point) {
- // FIXME Tony, could you please implement equals in VecOp or in Hep3Vector and fix this ?
- if (point == origin)
- return 0;
- // 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;
- double rho = -cosLambda/radius;
- // initial momentum
- Hep3Vector p0 = new BasicHep3Vector(px, py, pz);
+ double cosPhi = cos(phi);
+ double sinPhi = sin(phi);
+ double tanLambda = sinLambda/cosLambda;
+ double px = radius*cosPhi;
+ double py = radius*sinPhi;
+ double pz = radius*tanLambda;
+ double rho = -cosLambda/radius;
+ double pMag = sqrt(px*px + py*py + pz*pz);
+ 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 = VecOp.dot(xDiff, p0);
-
- 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 = factorP/2 + sqrt(factorP/2*factorP/2 + factorQ);
- // "minus case"
- 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);
-
- // see if the sufficient requirement of a minimum is met
-
- // Momentum at the new position x B Field;
- Hep3Vector newPCrossB = new BasicHep3Vector(newMomentum.y(), -newMomentum.x(), 0);
- // Spacial difference Vector at the new location
- Hep3Vector newXDiff = VecOp.sub(newPosition, point);
- double secondDerivative = pMag+rho*VecOp.dot(newXDiff, newPCrossB);
- if (secondDerivative > 0 )
- return newXDiff.magnitude();
- return VecOp.sub(getPointAtDistance(length2), point).magnitude();
+ // first, the point needs to be translated into the first period
+ Hep3Vector xDiff = VecOp.sub(origin, point);
+ double zPos = xDiff.z();
+ int addedPeriods = 0;
+ // FIXME translate to first period
+// while (zPos > tanLambda*2*Math.PI) {
+// zPos -= tanLambda*2*Math.PI;
+// ++addedPeriods;
+// }
+ double factorA1 = pMag - pz*pz/pMag - (VecOp.dot(xDiff, pCrossB))*rho;
+ double factorA2 = (VecOp.dot(xDiff, p0) - xDiff.z()*pz)*rho;
+
+ return Math.atan(factorA2/factorA1) / -rho;
+ }
+
+
+ /**
+ * 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 of closest approach
+ */
+ public double getSignedClosestDifferenceToPoint(Hep3Vector point) {
+ double cosPhi = cos(phi);
+ double sinPhi = sin(phi);
+ double tanLambda = sinLambda/cosLambda;
+ double px = radius*cosPhi;
+ double py = radius*sinPhi;
+ double pz = radius*tanLambda;
+ double rho = -cosLambda/radius;
+ Hep3Vector pCrossB = new BasicHep3Vector(py, -px, 0);
+ Hep3Vector xDiff = VecOp.sub(origin, point);
+ double pMag = sqrt(px*px + py*py + pz*pz);
+
+ double numerator = (-2*VecOp.dot(xDiff, pCrossB) + pMag*rho*(xDiff.magnitudeSquared()-xDiff.z()*xDiff.z()))/radius;
+ double denominator = 1 + sqrt(1-2*rho*pMag*VecOp.dot(xDiff, pCrossB)/radius/radius + pMag*pMag*rho*rho*(xDiff.magnitudeSquared()-xDiff.z()*xDiff.z())/radius/radius);
+ return numerator/denominator;
}
// 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
Hep3Vector getMomentumAtDistance(double s) {
double p0x = radius*cos(phi);
double p0y = radius*sin(phi);
@@ -144,7 +147,7 @@
double pz = p0z*cos(rho*s) + p0z*(1-cos(rho*s));
return new BasicHep3Vector(px, py, pz);
}
-
+
private Hep3Vector origin;
private double xCenter;
private double yCenter;
lcsim/test/org/lcsim/util/swim
diff -u -r1.4 -r1.5
--- HelixTest.java 25 Aug 2005 04:52:43 -0000 1.4
+++ HelixTest.java 26 Aug 2005 08:23:26 -0000 1.5
@@ -10,7 +10,7 @@
/**
*
* @author tonyj
- * @version $Id: HelixTest.java,v 1.4 2005/08/25 04:52:43 jstrube Exp $
+ * @version $Id: HelixTest.java,v 1.5 2005/08/26 08:23:26 jstrube Exp $
*/
public class HelixTest extends TestCase
{
@@ -42,13 +42,13 @@
assertTrue(Double.isNaN(circle.getDistanceToInfiniteCylinder(0)));
assertEquals(0,circle.getDistanceToInfiniteCylinder(1),1e-14);
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(1, 1, 1)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(-.5, 0,0)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(.5, 0,0)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,.5,0)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,-.5,0)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,.5,1)));
- System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,-.5,1)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(1, 1, 1)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(-.5, 0,0)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(.5, 0,0)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,.5,0)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,-.5,0)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,.5,1)));
+// System.err.println(circle.getDistanceFromHelixToPoint(new BasicHep3Vector(0,-.5,1)));
}
public void testCircle2()
{
@@ -135,6 +135,18 @@
aida.saveAs("helix.aida");
}
+ public void testGetDistanceFromHelixToPoint() {
+ Hep3Vector origin = new BasicHep3Vector(0,0,0);
+ double radius = 1;
+ double phi = Math.PI/2;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin,radius,phi,lambda);
+
+ Hep3Vector pointOnHelix = helix.getPointAtDistance(2.1);
+ System.out.println(helix.getDistanceFromHelixToPoint(pointOnHelix));
+ System.out.println(helix.getDistanceFromHelixToPoint(helix.getPointAtDistance(21.1)));
+ }
+
public void testGetDistanceToPoint() {
Hep3Vector origin = new BasicHep3Vector(0,0,0);
double radius = 1;
@@ -144,10 +156,10 @@
// assertEquals(helix.getDistanceFromHelixToPoint(origin), 0, 1e-14);
for (int i=0; i<100; ++i) {
- System.err.println("Trying to find" + i*.27);
+// System.err.println("Trying to find" + i*.27);
Hep3Vector pointOnHelix = helix.getPointAtDistance(i*.27);
Hep3Vector pointToFind = new BasicHep3Vector(pointOnHelix.x()+1, pointOnHelix.y(), pointOnHelix.z());
- System.err.println(helix.getDistanceFromHelixToPoint(pointToFind));
+// System.err.println(helix.getDistanceFromHelixToPoint(pointToFind));
}
// System.err.println("Trying to find 2.66");
// Hep3Vector pointOnHelix = helix.getPointAtDistance(2.66);
@@ -155,18 +167,38 @@
// System.err.println("Trying to find 20.66");
// pointOnHelix = helix.getPointAtDistance(20.66);
// helix.getDistanceFromHelixToPoint(pointOnHelix);
- System.err.println("Trying to find 3*Pi/2");
+// System.err.println("Trying to find 3*Pi/2");
Hep3Vector pointOnHelix = helix.getPointAtDistance(radius*3*Math.PI/2);
helix.getDistanceFromHelixToPoint(pointOnHelix);
- System.err.println("Trying to find 1.3");
+// System.err.println("Trying to find 1.3");
pointOnHelix = helix.getPointAtDistance(1.3);
- System.err.println(pointOnHelix);
+// System.err.println(pointOnHelix);
// assertEquals(helix.getDistanceFromHelixToPoint(pointOnHelix), 0, 1e-14);
Hep3Vector newPoint = new BasicHep3Vector(0, 1.0, 0);
// assertEquals(helix.getDistanceFromHelixToPoint(newPoint), 1.0, 1e-14);
Hep3Vector centerOfCircle = new BasicHep3Vector(1, 0, 0);
// assertEquals(helix.getDistanceFromHelixToPoint(centerOfCircle), radius, 1e-14);
}
+
+ public void testGetSignedClosestDifferenceToPoint() throws IOException {
+ Hep3Vector origin = new BasicHep3Vector(0,0,0);
+ double radius = 1;
+ double phi = Math.PI/2;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+ // test "random" points on helix
+ for (int i=0; i<500; ++i) {
+ // System.err.println("Trying to find" + i*.37);
+ Hep3Vector pointOnHelix = helix.getPointAtDistance(i*.27);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(pointOnHelix), 0, 1e-10);
+ }
+ Hep3Vector center = new BasicHep3Vector(radius, 0, 0);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(center), radius, 1e-10);
+ Hep3Vector halfway = new BasicHep3Vector(radius/2, 0, 0);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(halfway), radius/2, 1e-10);
+ Hep3Vector aboveCenter = new BasicHep3Vector(1, 2, 0);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(aboveCenter), -1, 1e-10);
+ }
private void assertEquals(Hep3Vector v1, Hep3Vector v2)
{
