lcsim/src/org/lcsim/util/swim
diff -u -r1.27 -r1.28
--- Helix.java 29 Nov 2007 02:25:55 -0000 1.27
+++ Helix.java 10 Feb 2010 21:57:30 -0000 1.28
@@ -25,7 +25,7 @@
* by Paul Avery.
*
* @author tonyj
- * @version $Id: Helix.java,v 1.27 2007/11/29 02:25:55 jstrube Exp $
+ * @version $Id: Helix.java,v 1.28 2010/02/10 21:57:30 partridge Exp $
*/
public class Helix implements Trajectory {
/**
@@ -118,6 +118,69 @@
return result;
}
+ public double getDistanceToPoint(Hep3Vector point) {
+
+ // Set starting position and direction unit vector
+ Hep3Vector pos = new BasicHep3Vector(origin.v());
+ Hep3Vector u = VecOp.unit(new BasicHep3Vector(px, py, pz));
+ Hep3Vector zhat = new BasicHep3Vector(0., 0., 1.);
+
+ // First estimate distance using z coordinate of the point
+ double s = getDistanceToZPlane(point.z());
+ double stot = s;
+
+ // Propagate to the estimated point
+ Hep3Vector unew = propagateDirection(u, s);
+ Hep3Vector posnew = propagatePosition(pos, u, s);
+ u = unew;
+ pos = posnew;
+
+ // Calculate how far we are from the desired point
+ Hep3Vector delta = VecOp.sub(pos, point);
+
+ // Iteratively close in on the point of closest approach
+ while (delta.magnitude() > eps) {
+
+ // Calculate the coefficients of the indicial equations a*s^2 + b*s + c = 0
+ double c = VecOp.dot(u, delta);
+ double b = 1. - rho * VecOp.dot(VecOp.cross(u, zhat), delta);
+ double a = -0.5 * rho*rho * (c - delta.z()*u.z());
+
+ // Find the two solutions
+ double arg = b*b - 4 * a * c;
+ double s1 = (-b + Math.sqrt(arg)) / (2. * a);
+ double s2 = (-b - Math.sqrt(arg)) / (2. * a);
+
+ // Find the position and distance from desired point for both solutions
+ Hep3Vector pos1 = propagatePosition(pos, u, s1);
+ double d1 = VecOp.sub(pos1, point).magnitude();
+ Hep3Vector pos2 = propagatePosition(pos, u, s2);
+ double d2 = VecOp.sub(pos2, point).magnitude();
+
+ // Pick the closest solution and update the position, direction unit vector, and
+ // path length. If the change in position is small, we have converged.
+ if (d1 < d2) {
+ unew = propagateDirection(u, s1);
+ u = unew;
+ pos = pos1;
+ stot += s1;
+ if (Math.abs(s1) < eps/100.) return stot;
+ } else {
+ unew = propagateDirection(u, s2);
+ u = unew;
+ pos = pos2;
+ stot += s2;
+ if (Math.abs(s2) < eps/100.) return stot;
+ }
+
+ // Update how far we are from the specified point
+ delta = VecOp.sub(pos, point);
+ }
+
+ // If we get here, we found a point within the desired precision
+ return stot;
+ }
+
/**
* Swims the helix along its trajectory to the point of closest approach to
* the given SpacePoint. For more info on swimming see <a
@@ -127,7 +190,7 @@
* Point in Space to swim to.
* @return the length Parameter alpha
*/
- public double getDistanceToPoint(Hep3Vector point) {
+ public double getDistanceToXYPosition(Hep3Vector point) {
double tanLambda = sinLambda / cosLambda;
double pMag = sqrt(px * px + py * py + pz * pz);
Hep3Vector p0 = new BasicHep3Vector(px, py, pz);
@@ -205,7 +268,7 @@
public Hep3Vector getTangentAtDistance(double alpha) {
double p0x = px * cos(rho * alpha) - py * sin(rho * alpha);
double p0y = py * cos(rho * alpha) + px * sin(rho * alpha);
- double p0z = pz * cos(rho * alpha) + pz * (1 - cos(rho * alpha));
+ double p0z = pz;
return new BasicHep3Vector(p0x, p0y, p0z);
}
@@ -243,6 +306,48 @@
}
/**
+ * Propagates the direction unit vector a distance s along the helix
+ *
+ * @param u initial direction unit vector
+ * @param s distance to propagate
+ * @return propagated direction unit vector
+ */
+ private Hep3Vector propagateDirection(Hep3Vector u, double s) {
+ double ux = u.x();
+ double uy = u.y();
+ double uz = u.z();
+ double carg = Math.cos(rho * s);
+ double sarg = Math.sin(rho * s);
+ double uxnew = ux * carg -uy * sarg;
+ double uynew = uy * carg + ux * sarg;
+ double uznew = uz;
+ return new BasicHep3Vector (uxnew, uynew, uznew);
+ }
+
+ /**
+ * Propagate a point on the helix by a specified distance
+ *
+ * @param pos starting position
+ * @param u starting direction unit vector
+ * @param s distance to propagate
+ * @return propagated point on helix
+ */
+ private Hep3Vector propagatePosition(Hep3Vector pos, Hep3Vector u, double s) {
+ double x = pos.x();
+ double y = pos.y();
+ double z = pos.z();
+ double ux = u.x();
+ double uy = u.y();
+ double uz = u.z();
+ double carg = Math.cos(rho * s);
+ double sarg = Math.sin(rho * s);
+ double xnew = x + ux * sarg / rho - uy * (1 - carg) / rho;
+ double ynew = y + uy * sarg / rho + ux * (1 - carg) / rho;
+ double znew = z + uz * s;
+ return new BasicHep3Vector(xnew, ynew, znew);
+ }
+
+ /**
* @param alpha
* The distance along the trajectory in the x-y plane
* @return The unit vector of the momentum
@@ -273,4 +378,7 @@
private double pz;
private double abs_r;
private double rho;
+
+ // Set the desired precision in finding the point closest to the track
+ private double eps = 1.e-6; // 1 nm ought to be good enough for government work!
}
lcsim/src/org/lcsim/util/swim
diff -u -r1.19 -r1.20
--- HelixSwimmer.java 29 Nov 2007 02:25:55 -0000 1.19
+++ HelixSwimmer.java 10 Feb 2010 21:57:30 -0000 1.20
@@ -11,14 +11,10 @@
import org.lcsim.spacegeom.SpaceVector;
import static java.lang.Math.atan;
-import static java.lang.Math.sin;
-import static java.lang.Math.cos;
import static java.lang.Math.sqrt;
import static org.lcsim.constants.Constants.fieldConversion;
-import static org.lcsim.event.LCIOParameters.ParameterName.d0;
import static org.lcsim.event.LCIOParameters.ParameterName.phi0;
import static org.lcsim.event.LCIOParameters.ParameterName.omega;
-import static org.lcsim.event.LCIOParameters.ParameterName.z0;
import static org.lcsim.event.LCIOParameters.ParameterName.tanLambda;
/**
@@ -29,7 +25,7 @@
* by Paul Avery.
*
* @author jstrube
- * @version $Id: HelixSwimmer.java,v 1.19 2007/11/29 02:25:55 jstrube Exp $
+ * @version $Id: HelixSwimmer.java,v 1.20 2010/02/10 21:57:30 partridge Exp $
*/
public class HelixSwimmer {
private double field;
