lcsim/src/org/lcsim/util/swim
diff -u -r1.6 -r1.7
--- Line.java 28 Sep 2005 22:39:57 -0000 1.6
+++ Line.java 18 Oct 2005 19:08:48 -0000 1.7
@@ -7,7 +7,7 @@
/**
* A straight line
* @author tonyj
- * @version $Id: Line.java,v 1.6 2005/09/28 22:39:57 mcharles Exp $
+ * @version $Id: Line.java,v 1.7 2005/10/18 19:08:48 mcharles Exp $
*/
public class Line implements Trajectory
{
@@ -251,4 +251,150 @@
Hep3Vector lineDirection = new BasicHep3Vector(cosLambda*cosPhi, cosLambda*sinPhi, sinLambda);
return lineDirection;
}
+
+ /**
+ * Return the direction of the line. Marked Deprecated since this should
+ * really be replaced by a general method for the Trajectory interface,
+ * perhaps getDirection(double alpha) for distance parameter alpha.
+ */
+ @Deprecated public Hep3Vector getDirection()
+ {
+ return this.getUnitVector();
+ }
+
+ /**
+ * Routine to find where two lines meet.
+ * This will probably be replaced by a more general vertexing solution
+ * at some point, but for now it works.
+ * @param
+ */
+ @Deprecated static public double[] getPOCAOfLines(Line line1, Line line2)
+ {
+ double[] output = null;
+ try {
+ output = line1.findTrackToTrackPOCA(line2);
+ return output;
+ } catch (ParallelLinesException x) {
+ // Er -- hm.
+ return null;
+ }
+ }
+
+ /**
+ * Internal utility routine to find where two lines meet.
+ * This will probably be replaced by a more general vertexing solution
+ * at some point, but for now it works.
+ */
+ @Deprecated private double[] findTrackToTrackPOCA(Line otherLine) throws ParallelLinesException
+ {
+ Line line1 = this;
+ Line line2 = otherLine;
+
+ // Unit vectors of the two lines:
+ Hep3Vector v1 = line1.getUnitVector();
+ Hep3Vector v2 = line2.getUnitVector();
+ // Origin points on the two lines:
+ Hep3Vector x1 = line1.origin;
+ Hep3Vector x2 = line2.origin;
+
+ // Find the common perpendicular:
+ Hep3Vector perp = VecOp.cross(v1, v2);
+ if (perp.magnitude() == 0.0) {
+ // Lines are parallel!
+ throw new ParallelLinesException("Parallel lines");
+ } else {
+ // Let the origin point along the lines be x1 and x2.
+ // Let the unit vectors along the lines be v1 and v2.
+ // Suppose that the points of closest approach are at:
+ // x1' = x1 + a(v1)
+ // x2' = x2 + b(v2)
+ // Construct a vector p1 which is perpendicular to perp and to v1:
+ // p1 = v1 x perp = v1 x (v1 x v2) = v1(v1.v2) - v2(v1.v1)
+ // p2 = v2 x perp = v2 x (v1 x v2) = v1(v2.v2) - v2(v1.v2)
+ // Then v1.p1 = (v1.v1)(v1.v2) - (v1.v2)(v1.v1) = 0
+ // v1.p2 = (v1.v1)(v2.v2) - (v1.v2)(v1.v2) = L
+ // v2.p1 = (v2.v1)(v1.v2) - (v2.v2)(v1.v1) = -L
+ // v2.p2 = (v2.v1)(v2.v2) - (v2.v2)(v1.v2) = 0
+ // Thus,
+ // (x2'-x1').p1 = (x2 + b(v2) - x1 -a(v1)).p1
+ // = x2.p1 + b(v2.p1) - x1.p1 - a(v1.p1)
+ // = x2.p1 - bL - x1.p1
+ // and similarly,
+ // (x2'-x1').p2 = x2.p2 -x1.p2 - aL
+ // But p1 and p2 are perpendicular to perp, and (x2'-x1') is parallel
+ // to perp for x1' and x2' the closest points. So
+ // (x2'-x1').p1 = x2.p1 - bL - x1.p1 = 0
+ // (x2'-x1').p2 = x2.p2 - x1.p2 - aL = 0
+ // Hence,
+ // -bL = -x2.p1 + x1.p1 => b = (x2-x1).p1 / L
+ // -aL = -x2.p2 + x1.p2 => a = (x2-x1).p2 / L
+ // and so the POCA is at
+ // (x1' + x2')/2 = (x1 + a(v1) + x2 + b(v2)) / 2
+ // = (x1 + [(x2-x1).p2/L]v1 + x2 + [(x2-x1).p1/L]v2) / 2
+
+ // p1 = v1 x perp (and similarly p2)
+ Hep3Vector p1 = VecOp.cross(v1, perp);
+ Hep3Vector p2 = VecOp.cross(v2, perp);
+
+ // L = (v1.v1)(v2.v2) - (v1.v2)(v1.v2) = 1 - (v1.v2)^2
+ double L = 1.0 - (VecOp.dot(v1,v2) * VecOp.dot(v1,v2));
+
+ // a = (x2-x1).p2 / L
+ // b = (x2-x1).p1 / L
+ double a = VecOp.dot(VecOp.sub(x2,x1), p2) / L;
+ double b = VecOp.dot(VecOp.sub(x2,x1), p1) / L;
+
+ double[] output = new double[2];
+ output[0] = a;
+ output[1] = b;
+ return output;
+ }
+ }
+
+ @Deprecated protected class ParallelLinesException extends Exception {
+ public ParallelLinesException(String m) { super(m); }
+ }
+
+ /**
+ * Defines a line in space. The line is defined by a point and a
+ * direction vector. The magnitude of the direction vector is ignored.
+ * @param origin A point on the line
+ * @param dir The signed direction of the line
+ */
+ public Line(Hep3Vector origin, Hep3Vector dir)
+ {
+ // The easy bit:
+ this.origin = origin;
+
+ // The hard bit: find the direction as (cosPhi, sinPhi, cosLambda, sinLambda)
+ double normalization = dir.magnitude();
+ Hep3Vector dirNormalized = hep.physics.vec.VecOp.unit(dir);
+
+ // cosPhi, sinPhi, sinLambda are fairly easy.
+ // It helps that the signed unit vector is (cos(lambda)*cos(phi), cos(lambda)*sin(phi), sin(lambda)).
+ double phi = hep.physics.vec.VecOp.phi(dirNormalized);
+ cosPhi = Math.cos(phi);
+ sinPhi = Math.sin(phi);
+ sinLambda = dirNormalized.z();
+
+ // cosLambda is harder. First find the absolute value:
+ double cosLambdaSquared = (dirNormalized.x()*dirNormalized.x() + dirNormalized.y()*dirNormalized.y());
+ if (cosLambdaSquared == 0.0) {
+ // Easy case
+ cosLambda = 0.0;
+ } else {
+ // Need to tease it apart using either X or Y. Choose the one with the bigger lever arm...
+ boolean useX = (Math.abs(dir.x()) > Math.abs(dir.y()));
+
+ // We will use one of the relations:
+ // vx = cosLambda * cosPhi
+ // vy = cosLambda * sinPhi
+ // where the unit direction vector is (vx, vy, vz)
+ if (useX) {
+ cosLambda = dirNormalized.x() / cosPhi;
+ } else {
+ cosLambda = dirNormalized.y() / sinPhi;
+ }
+ }
+ }
}