lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -u -r1.4 -r1.5
--- FitterTestDriver.java 10 Sep 2006 11:47:39 -0000 1.4
+++ FitterTestDriver.java 12 Sep 2006 13:28:23 -0000 1.5
@@ -12,11 +12,13 @@
import org.lcsim.event.MCParticle;
import org.lcsim.spacegeom.CartesianPoint;
import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpaceVector;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
import static java.lang.Math.abs;
-
+import static java.lang.Math.sqrt;
class FitterTestDriver extends Driver {
private AIDA aida = AIDA.defaultInstance();
private SpacePoint referencePoint;
@@ -42,68 +44,80 @@
continue;
List<Track> unsmearedInput = new ArrayList<Track>();
List<Track> smearedInput = new ArrayList<Track>();
+ origin = new SpacePoint(part.getEndPoint());
for (MCParticle daughter : part.getDaughters()) {
Track unsmearedTrack = tf.getTrack(
- new SpacePoint(daughter.getMomentum()),
+ new SpaceVector(daughter.getMomentum()),
new SpacePoint(daughter.getOrigin()),
- referencePoint,
+ origin,
(int)daughter.getCharge(),
rand,
false);
- System.err.println("Particle: charge " + daughter.getCharge());
- System.err.println("Position:" + new SpacePoint(daughter.getOrigin()));
- System.err.println("Momentum:" + new SpacePoint(daughter.getMomentum()));
- SpacePoint unsmearedPosition = NewTrack.Parameters2Position(unsmearedTrack.getParameters());
+// System.err.println("Particle: charge " + daughter.getCharge());
+// System.err.println("Position:" + new SpacePoint(daughter.getOrigin()));
+// System.err.println("Momentum:" + new SpacePoint(daughter.getMomentum()));
+ SpacePoint unsmearedPosition = EMap.Parameters2Position(unsmearedTrack.getParameters(), referencePoint);
SpacePoint unsmearedMomentum = EMap.Parameters2Momentum(unsmearedTrack.getParameters());
- System.err.println("unsmeared Position: " + unsmearedPosition);
- System.err.println("unsmeared Momentum: " + unsmearedMomentum);
+// System.err.println("unsmeared Position: " + unsmearedPosition);
+// System.err.println("unsmeared Momentum: " + unsmearedMomentum);
unsmearedInput.add(unsmearedTrack);
- aida.cloud1D("unsmeared_pt").fill(unsmearedTrack.getPt());
- aida.cloud1D("unsmeared_diffx").fill(unsmearedPosition.x() - daughter.getOriginX());
- aida.cloud1D("unsmeared_diffy").fill(unsmearedPosition.y() - daughter.getOriginY());
- aida.cloud1D("unsmeared_diffz").fill(unsmearedPosition.z() - daughter.getOriginZ());
- aida.cloud1D("unsmeared_diffpx").fill(unsmearedMomentum.x() - daughter.getPX());
- aida.cloud1D("unsmeared_diffpy").fill(unsmearedMomentum.y() - daughter.getPY());
- aida.cloud1D("unsmeared_diffpz").fill(unsmearedMomentum.z() - daughter.getPZ());
+ aida.cloud1D("track_unsmeared_pt").fill(unsmearedTrack.getPt());
+ aida.cloud1D("track_unsmeared_diffx").fill(unsmearedPosition.x() - daughter.getOriginX());
+ aida.cloud1D("track_unsmeared_diffy").fill(unsmearedPosition.y() - daughter.getOriginY());
+ aida.cloud1D("track_unsmeared_diffz").fill(unsmearedPosition.z() - daughter.getOriginZ());
+ aida.cloud1D("track_unsmeared_diffpx").fill(unsmearedMomentum.x() - daughter.getPX());
+ aida.cloud1D("track_unsmeared_diffpy").fill(unsmearedMomentum.y() - daughter.getPY());
+ aida.cloud1D("track_unsmeared_diffpz").fill(unsmearedMomentum.z() - daughter.getPZ());
Track smearedTrack = tf.getTrack(
- new SpacePoint(daughter.getMomentum()),
+ new SpaceVector(daughter.getMomentum()),
new SpacePoint(daughter.getOrigin()),
- referencePoint,
+ origin,
(int)daughter.getCharge(),
rand,
true);
- SpacePoint smearedPosition = NewTrack.Parameters2Position(smearedTrack.getParameters());
+ SpacePoint smearedPosition = EMap.Parameters2Position(smearedTrack.getParameters(), referencePoint);
SpacePoint smearedMomentum = EMap.Parameters2Momentum(smearedTrack.getParameters());
- System.err.println("smeared Position: " + smearedPosition);
- System.err.println("smeared Momentum: " + smearedMomentum);
+// System.err.println("smeared Position: " + smearedPosition);
+// System.err.println("smeared Momentum: " + smearedMomentum);
smearedInput.add(smearedTrack);
- aida.cloud1D("smeared_pt").fill(smearedTrack.getPt());
- aida.cloud1D("smeared_diffx").fill(smearedPosition.x() - daughter.getOriginX());
- aida.cloud1D("smeared_diffy").fill(smearedPosition.y() - daughter.getOriginY());
- aida.cloud1D("smeared_diffz").fill(smearedPosition.z() - daughter.getOriginZ());
- aida.cloud1D("smeared_diffpx").fill(smearedMomentum.x() - daughter.getPX());
- aida.cloud1D("smeared_diffpy").fill(smearedMomentum.y() - daughter.getPY());
- aida.cloud1D("smeared_diffpz").fill(smearedMomentum.z() - daughter.getPZ());
- aida.cloud1D("diff_d0").fill(unsmearedTrack.getParameters().getValues()[0]-smearedTrack.getParameters().getValues()[0]);
- aida.cloud1D("diff_phi0").fill(unsmearedTrack.getParameters().getValues()[1]-smearedTrack.getParameters().getValues()[1]);
- aida.cloud1D("diff_omega").fill(unsmearedTrack.getParameters().getValues()[2]-smearedTrack.getParameters().getValues()[2]);
- aida.cloud1D("diff_z").fill(unsmearedTrack.getParameters().getValues()[3]-smearedTrack.getParameters().getValues()[3]);
- aida.cloud1D("diff_tanLambda").fill(unsmearedTrack.getParameters().getValues()[4]-smearedTrack.getParameters().getValues()[4]);
+ aida.cloud1D("track_smeared_pt").fill(smearedTrack.getPt());
+ aida.cloud1D("track_smeared_diffx").fill(smearedPosition.x() - daughter.getOriginX());
+ aida.cloud1D("track_smeared_diffy").fill(smearedPosition.y() - daughter.getOriginY());
+ aida.cloud1D("track_smeared_diffz").fill(smearedPosition.z() - daughter.getOriginZ());
+ aida.cloud1D("track_smeared_diffpx").fill(smearedMomentum.x() - daughter.getPX());
+ aida.cloud1D("track_smeared_diffpy").fill(smearedMomentum.y() - daughter.getPY());
+ aida.cloud1D("track_smeared_diffpz").fill(smearedMomentum.z() - daughter.getPZ());
+ aida.cloud1D("track_diff_d0").fill(unsmearedTrack.getParameters().getValues()[0]-smearedTrack.getParameters().getValues()[0]);
+ aida.cloud1D("track_diff_phi0").fill(unsmearedTrack.getParameters().getValues()[1]-smearedTrack.getParameters().getValues()[1]);
+ aida.cloud1D("track_diff_omega").fill(unsmearedTrack.getParameters().getValues()[2]-smearedTrack.getParameters().getValues()[2]);
+ aida.cloud1D("track_diff_z").fill(unsmearedTrack.getParameters().getValues()[3]-smearedTrack.getParameters().getValues()[3]);
+ aida.cloud1D("track_diff_tanLambda").fill(unsmearedTrack.getParameters().getValues()[4]-smearedTrack.getParameters().getValues()[4]);
}
Fitter fitter = new Fitter(event);
Vertex newVtx1 = fitter.fit(unsmearedInput, origin);
Fitter fitter2 = new Fitter(event);
Vertex newVtx2 = fitter2.fit(smearedInput, origin);
+ double sigmaX = sqrt(newVtx2.getSpatialCovarianceMatrix().get(0, 0));
+ double sigmaY = sqrt(newVtx2.getSpatialCovarianceMatrix().get(1, 1));
+ double sigmaZ = sqrt(newVtx2.getSpatialCovarianceMatrix().get(2, 2));
+ double deltaX = newVtx2.origin().x() - part.getOrigin().x();
+ double deltaY = newVtx2.origin().y() - part.getOrigin().y();
+ double deltaZ = newVtx2.origin().z() - part.getOrigin().z();
aida.cloud1D("decayLength").fill(new SpacePoint(part.getEndPoint()).magnitude());
- aida.cloud1D("unsmeared_deltaX").fill(newVtx1.origin().x() - part.getOrigin().x());
- aida.cloud1D("unsmeared_deltaY").fill(newVtx1.origin().y() - part.getOrigin().y());
- aida.cloud1D("unsmeared_deltaZ").fill(newVtx1.origin().z() - part.getOrigin().z());
- aida.cloud1D("smeared_deltaX").fill(newVtx2.origin().x() - part.getOrigin().x());
- aida.cloud1D("smeared_deltaY").fill(newVtx2.origin().y() - part.getOrigin().y());
- aida.cloud1D("smeared_deltaZ").fill(newVtx2.origin().z() - part.getOrigin().z());
+ aida.cloud1D("vtx_smeared_deltaX").fill(deltaX);
+ aida.cloud1D("vtx_smeared_deltaY").fill(deltaY);
+ aida.cloud1D("vtx_smeared_deltaZ").fill(deltaZ);
+ aida.cloud1D("vtx_sigma_x").fill(sigmaX);
+ aida.cloud1D("vtx_sigma_y").fill(sigmaY);
+ aida.cloud1D("vtx_sigma_z").fill(sigmaZ);
+ aida.cloud1D("vtx_pull_x").fill(deltaX/sigmaX);
+ aida.cloud1D("vtx_pull_y").fill(deltaY/sigmaY);
+ aida.cloud1D("vtx_pull_z").fill(deltaZ/sigmaZ);
+ aida.cloud1D("vtx_unsmeared_deltaX").fill(newVtx1.origin().x() - part.getOrigin().x());
+ aida.cloud1D("vtx_unsmeared_deltaY").fill(newVtx1.origin().y() - part.getOrigin().y());
+ aida.cloud1D("vtx_unsmeared_deltaZ").fill(newVtx1.origin().z() - part.getOrigin().z());
}
}
-
}
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.5 -r1.6
--- Helix.java 10 Sep 2006 11:47:35 -0000 1.5
+++ Helix.java 12 Sep 2006 13:28:23 -0000 1.6
@@ -17,187 +17,239 @@
import org.lcsim.spacegeom.SpaceVector;
/**
- * 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).
+ * 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).
* <p>
- * For more info on swimming see <a href="doc-files/transport.pdf">this paper</a> by Paul Avery.
+ * 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.5 2006/09/10 11:47:35 jstrube Exp $
+ * @version $Id: Helix.java,v 1.6 2006/09/12 13:28:23 jstrube Exp $
*/
public class Helix implements Trajectory {
- private 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.
- private 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 The <em>signed</em> radius of curvature of the helix. The conventions is such that for
- * <em>positive</em> radii, the direction is <em>clockwise</em>.
- * @param phi The polar angle of the helix <em>momentum</em> in x-y plane w/rt positive x-axis at the origin
- * @param lambda The dip angle w/rt positive part of the x-y plane
- */
- 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);
- cosLambda = cos(lambda);
- sinLambda = sin(lambda);
- xCenter = origin.x() + radius * sinPhi;
- 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) {
- 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() {
- 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) {
- double phiToCenter = atan2(yCenter, xCenter);
- double radiusOfCenter = sqrt(xCenter * xCenter + yCenter * yCenter);
- // 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;
- while (result < 0)
- result += Math.abs(radius / cosLambda) * 2 * Math.PI;
- return result;
- }
-
- 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) {
- double angle = phi + alpha * rho;
- return new CartesianVector(cos(angle), sin(angle), sinLambda / cosLambda);
- }
-
- 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 of closest approach
- */
- // 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) {
- 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
- double zPos = xDiff.z();
- // these are two mutually exclusive cases and two while loops
- // may not be the best way to express this
- while (zPos > abs(radius * tanLambda * Math.PI / 4)) {
- zPos -= abs(radius * tanLambda * Math.PI / 4);
- }
- while (zPos < -abs(radius * tanLambda * Math.PI / 4)) {
- 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;
- 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;
- }
-
- /**
- * 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);
- 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 -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;
- }
+ 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
+ * The <em>signed</em> radius of curvature of the helix. The
+ * conventions is such that for <em>positive</em> radii, the
+ * direction is <em>clockwise</em>.
+ * @param phi
+ * The polar angle of the helix <em>momentum</em> in x-y plane
+ * w/rt positive x-axis at the origin
+ * @param lambda
+ * The dip angle w/rt positive part of the x-y plane
+ */
+ 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);
+ cosLambda = cos(lambda);
+ sinLambda = sin(lambda);
+ xCenter = origin.x() + radius * sinPhi;
+ 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) {
+ 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() {
+ 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) {
+ double phiToCenter = atan2(yCenter, xCenter);
+ double radiusOfCenter = sqrt(xCenter * xCenter + yCenter * yCenter);
+ // 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;
+ while (result < 0)
+ result += Math.abs(radius / cosLambda) * 2 * Math.PI;
+ return result;
+ }
+
+ 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) {
+ double angle = phi + alpha * rho;
+ return new CartesianVector(cos(angle), sin(angle), sinLambda
+ / cosLambda);
+ }
+
+ 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
+ * of closest approach
+ */
+ // 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) {
+ 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
+ double zPos = xDiff.z();
+ // these are two mutually exclusive cases and two while loops
+ // may not be the best way to express this
+ while (zPos > abs(radius * tanLambda * Math.PI / 4)) {
+ zPos -= abs(radius * tanLambda * Math.PI / 4);
+ }
+ while (zPos < -abs(radius * tanLambda * Math.PI / 4)) {
+ 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;
+ 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;
+ }
+
+ /**
+ * 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;
+ }
}
lcsim/test/org/lcsim/contrib/JanStrube/tracking
diff -N HelixTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixTest.java 12 Sep 2006 13:28:24 -0000 1.1
@@ -0,0 +1,356 @@
+package org.lcsim.contrib.JanStrube.tracking;
+
+import static java.lang.Math.abs;
+import static java.lang.Math.cos;
+import static java.lang.Math.sin;
+import hep.aida.ICloud2D;
+import junit.framework.*;
+import java.io.IOException;
+
+import org.lcsim.spacegeom.CartesianPoint;
+import org.lcsim.spacegeom.CartesianVector;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.aida.AIDA;
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+/**
+ *
+ * @author tonyj
+ * @version $Id: HelixTest.java,v 1.1 2006/09/12 13:28:24 jstrube Exp $
+ */
+public class HelixTest extends TestCase
+{
+
+ public HelixTest(String testName)
+ {
+ super(testName);
+ }
+
+ public static Test suite()
+ {
+ return new TestSuite(HelixTest.class);
+ }
+
+ public void testCircle()
+ {
+ SpaceVector origin = new CartesianVector(1,0,0);
+ double radius = -1;
+ double phi = Math.PI/2;
+ double lambda = 0;
+ Helix circle = new Helix(origin,radius,phi,lambda);
+
+ assertEquals(origin, circle.getPointAtLength(0));
+ assertEquals(origin, circle.getPointAtLength(radius*Math.PI*2));
+ assertEquals(new CartesianVector(-1,0,0), circle.getPointAtLength(Math.abs(radius*Math.PI)));
+ assertEquals(new CartesianVector(0,1,0), circle.getPointAtLength(Math.abs(radius*Math.PI/2)));
+
+ assertTrue(Double.isInfinite(circle.getDistanceToZPlane(1)));
+ assertTrue(Double.isNaN(circle.getDistanceToInfiniteCylinder(2)));
+ assertTrue(Double.isNaN(circle.getDistanceToInfiniteCylinder(0)));
+ assertEquals(0,circle.getDistanceToInfiniteCylinder(1),1e-14);
+ }
+
+ public void testCircle2()
+ {
+ SpacePoint origin = new SpacePoint();
+ double radius = 1;
+ double phi = Math.PI/2;
+ double lambda = 0;
+ Helix circle = new Helix(origin,radius,phi,lambda);
+
+ assertEquals(origin, circle.getPointAtLength(0));
+ assertEquals(origin, circle.getPointAtLength(radius*Math.PI*2));
+ assertEquals(new CartesianVector(2,0,0), circle.getPointAtLength(radius*Math.PI));
+ assertEquals(new CartesianVector(1,1,0), circle.getPointAtLength(radius*Math.PI/2));
+
+ assertTrue(Double.isInfinite(circle.getDistanceToZPlane(1)));
+ assertTrue(Double.isNaN(circle.getDistanceToInfiniteCylinder(3*radius)));
+ assertEquals(radius*Math.PI,circle.getDistanceToInfiniteCylinder(2*radius),1e-14);
+ assertEquals(0,circle.getDistanceToInfiniteCylinder(0),1e-14);
+ }
+
+ public void testHelix()
+ {
+ SpacePoint origin = new SpacePoint();
+ double radius = 1;
+ double phi = Math.PI/2;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin,radius,phi,lambda);
+
+ assertEquals(origin, helix.getPointAtLength(0));
+ double d = radius*Math.PI*2*Math.sqrt(2);
+ assertEquals(new CartesianVector(0,0,Math.PI*2), helix.getPointAtLength(d));
+ assertEquals(new CartesianVector(2,0,Math.PI), helix.getPointAtLength(d/2));
+
+ assertEquals(d/2,helix.getDistanceToZPlane(Math.PI));
+ assertEquals(d,helix.getDistanceToZPlane(Math.PI*2));
+ assertTrue(Double.isNaN(helix.getDistanceToInfiniteCylinder(3*radius)));
+ assertEquals(d/2,helix.getDistanceToInfiniteCylinder(2*radius),1e-14);
+ assertEquals(1.4809609793861218,helix.getDistanceToInfiniteCylinder(radius),1e-14);
+ assertEquals(0,helix.getDistanceToInfiniteCylinder(0),1e-14);
+ }
+
+ public void testHelix2()
+ {
+ SpacePoint origin = new SpacePoint();
+ double radius = -1;
+ double phi = Math.PI/2;
+ double lambda = -Math.PI/4;
+ Helix helix = new Helix(origin,radius,phi,lambda);
+
+ assertEquals(origin, helix.getPointAtLength(0));
+ double d = Math.PI*2*Math.sqrt(2);
+ assertEquals(new CartesianVector(0,0,-Math.PI*2), helix.getPointAtLength(d));
+ assertEquals(new CartesianVector(-2,0,-Math.PI), helix.getPointAtLength(d/2));
+
+ assertEquals(d/2,helix.getDistanceToZPlane(-Math.PI));
+ assertEquals(d,helix.getDistanceToZPlane(-Math.PI*2));
+ assertTrue(Double.isNaN(helix.getDistanceToInfiniteCylinder(3)));
+ assertEquals(d/2,helix.getDistanceToInfiniteCylinder(2),1e-14);
+ assertEquals(1.4809609793861218,helix.getDistanceToInfiniteCylinder(1),1e-14);
+ assertEquals(0,helix.getDistanceToInfiniteCylinder(0),1e-14);
+ }
+
+ public void testHelix3() throws IOException
+ {
+ SpacePoint origin = new SpacePoint();
+ double radius = 1;
+ double phi = Math.PI/2;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin,radius,phi,lambda);
+
+ double d = radius*Math.PI*2*Math.sqrt(2);
+
+ AIDA aida = AIDA.defaultInstance();
+ ICloud2D xy = aida.cloud2D("xy");
+ ICloud2D rz = aida.cloud2D("rz");
+ for (int i=0; i<100; i++)
+ {
+ double alpha = i*d/100;
+ SpacePoint point = helix.getPointAtLength(alpha);
+ xy.fill(point.x(),point.y());
+ double r = Math.sqrt(point.x()*point.x()+point.y()*point.y());
+ rz.fill(r,point.z());
+ }
+ assertEquals(1,xy.meanX(),1e-14);
+ assertEquals(0,xy.meanY(),1e-14);
+ assertEquals(Math.sqrt(2)/2,xy.rmsX(),1e-14);
+ assertEquals(Math.sqrt(2)/2,xy.rmsY(),1e-14);
+ aida.saveAs("helix.aida");
+ }
+
+ public void testGetDistanceFromHelixToPoint() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = 1.7;
+ double phi = Math.PI/4;
+ 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) {
+ SpacePoint pointOnHelix = helix.getPointAtLength(i*.27);
+ Hep3Vector xDiff = VecOp.sub(helix.origin, pointOnHelix);
+ double rho = -cos(lambda)/helix.getRadius();
+ if (abs(xDiff.z()) > abs(helix.getRadius()*Math.tan(lambda)*Math.PI)) {
+ System.out.printf("NOW:%f %f\n%s\n%s",abs(xDiff.z()), abs(helix.getRadius()*Math.tan(lambda)*Math.PI),helix.origin,pointOnHelix);
+ }
+
+ System.out.printf("%f\t%f\n", i*.27, Math.PI*radius/cos(lambda));
+ // System.err.printf("%.3f, found: %s\n", i*.27, helix.getDistanceToPoint(pointOnHelix));
+ assertEquals(helix.getTrackLengthToPoint(pointOnHelix), i*.27, 1e-10);
+ }
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+ double offset = helix.getTrackLengthToPoint(new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()));
+ for (int i=0; i<5; ++i) {
+ SpacePoint newPoint = new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ assertEquals(helix.getTrackLengthToPoint(newPoint), offset+i*period, 1e-10);
+ }
+ }
+
+ // same as before, neg. alpha
+ public void testGetDistanceFromNegHelixToPoint() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = 1.7;
+ double phi = Math.PI/4;
+ 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) {
+ SpacePoint pointOnHelix = helix.getPointAtLength(-i*.27);
+// System.err.printf("%.3f, found: %s\n", -i*.27, helix.getDistanceToPoint(pointOnHelix));
+ assertEquals(helix.getTrackLengthToPoint(pointOnHelix), -i*.27, 1e-10);
+ }
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+ double offset = helix.getTrackLengthToPoint(new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()));
+ for (int i=0; i<5; ++i) {
+ SpacePoint newPoint = new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ assertEquals(helix.getTrackLengthToPoint(newPoint), offset+i*period, 1e-10);
+ }
+ }
+
+ // same as before, neg. radius
+ public void testGetDistanceFromNegHelixRadiusToPoint() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = -1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+// System.out.println(period);
+ // test "random" points on helix
+ for (int i=0; i<500; ++i) {
+ SpacePoint pointOnHelix = helix.getPointAtLength(-i*.27);
+// System.out.println(pointOnHelix.z());
+// System.err.printf("%.3f, found: %s\n", -i*.27, helix.getDistanceToPoint(pointOnHelix));
+ assertEquals(helix.getTrackLengthToPoint(pointOnHelix), -i*.27, 1e-10);
+ }
+
+ // test that the distance to points has the right period along z
+ double offset = helix.getTrackLengthToPoint(new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()));
+ for (int i=0; i<5; ++i) {
+ SpacePoint newPoint = new CartesianPoint(origin.x()+0.333, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ assertEquals(helix.getTrackLengthToPoint(newPoint), offset+i*period, 1e-10);
+ }
+ }
+
+ public void testGetSignedClosestDifferenceToPoint() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = 1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+
+ final double offset1 = 0.333;
+ final double offset2 = 1.23;
+ final double offset3 = .77;
+ double point1 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x()+offset1, origin.y(), origin.z()));
+ double point2 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x(), origin.y()+offset2, origin.z()));
+ double point3 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x()+offset3, origin.y()+offset3, origin.z()));
+
+ for (int i=0; i<5; ++i) {
+ SpacePoint newPoint1 = new CartesianPoint(origin.x()+offset1, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ SpacePoint newPoint2 = new CartesianPoint(origin.x(), origin.y()+offset2, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ SpacePoint newPoint3 = new CartesianPoint(origin.x()+offset3, origin.y()+offset3, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+// System.out.println("Point " + i + newPoint.toString());
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint1), point1, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint2), point2, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint3), point3, 1e-10);
+ }
+ // move along the center and check for constant distance
+ SpacePoint originCenter = new CartesianPoint(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), origin.z());
+ double dist0 = helix.getSignedClosestDifferenceToPoint(originCenter);
+ for (int i=-250; i<250; ++i) {
+ SpacePoint pointOnCenter = new CartesianPoint(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), i/2);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(pointOnCenter), dist0, 1e-10);
+ }
+ }
+
+ // same as before, neg. radius
+ public void testGetSignedClosestDifferenceToPointNeg() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = -1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix = new Helix(origin, radius, phi, lambda);
+
+ // test that the distance to points has the right period along z
+ double period = 2*Math.PI*radius/Math.cos(lambda);
+
+ final double offset1 = 0.333;
+ final double offset2 = 1.23;
+ final double offset3 = .77;
+ double point1 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x()+offset1, origin.y(), origin.z()));
+ double point2 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x(), origin.y()+offset2, origin.z()));
+ double point3 = helix.getSignedClosestDifferenceToPoint(new CartesianPoint(origin.x()+offset3, origin.y()+offset3, origin.z()));
+
+ for (int i=0; i<5; ++i) {
+ SpacePoint newPoint1 = new CartesianPoint(origin.x()+offset1, origin.y(), origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ SpacePoint newPoint2 = new CartesianPoint(origin.x(), origin.y()+offset2, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+ SpacePoint newPoint3 = new CartesianPoint(origin.x()+offset3, origin.y()+offset3, origin.z()+i*radius*Math.tan(lambda)*Math.PI*2);
+// System.out.println("Point " + i + newPoint.toString());
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint1), point1, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint2), point2, 1e-10);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(newPoint3), point3, 1e-10);
+ }
+ // move along the center and check for constant distance
+ SpacePoint originCenter = new CartesianPoint(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), origin.z());
+ double dist0 = helix.getSignedClosestDifferenceToPoint(originCenter);
+ for (int i=-250; i<250; ++i) {
+ SpacePoint pointOnCenter = new CartesianPoint(origin.x()-radius*cos(phi), origin.y()+radius*sin(phi), i/2);
+ assertEquals(helix.getSignedClosestDifferenceToPoint(pointOnCenter), dist0, 1e-10);
+ }
+ }
+
+ public void testGetTangentAtDistance() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = -1.7;
+ double phi = Math.PI/4;
+ double lambda = Math.PI/4;
+ Helix helix1 = new Helix(origin, radius, phi, lambda);
+ Helix helix2 = new Helix(origin, -radius, phi, lambda);
+
+ for (int i=0; i<25; ++i) {
+ SpacePoint tangentAtThisPiece_1 = helix1.getUnitTangentAtLength(i*.37);
+ SpacePoint tangentAtThisPiece_2 = helix2.getUnitTangentAtLength(i*.37);
+ // momentum conservation
+ assertEquals(tangentAtThisPiece_1.magnitude(), helix1.getUnitTangentAtLength((i+1)*.37).magnitude(), 1e-10);
+ assertEquals(tangentAtThisPiece_2.magnitude(), helix2.getUnitTangentAtLength(i*.37).magnitude());
+ // both positively and negatively charged tracks must be moving forward
+ assertEquals(tangentAtThisPiece_1.z(), tangentAtThisPiece_2.z());
+ }
+ }
+
+ public void testCurvatureSign() {
+ SpacePoint origin = new CartesianPoint(3, 6, 9);
+ double radius = 1.5;
+ double phi = Math.PI/2;
+ double lambda = 1.3;
+ for (int iPhi=0; iPhi<50; ++iPhi) {
+ Helix helix1 = new Helix(origin, radius, iPhi*phi, lambda);
+ Helix helix2 = new Helix(origin, -radius, iPhi*phi, lambda);
+ Helix helix5 = new Helix(origin, radius, Math.PI-iPhi*phi, lambda);
+ Helix helix6 = new Helix(origin, -radius, Math.PI-iPhi*phi, lambda);
+ assertEquals(helix1.yCenter, helix6.yCenter, 1e-10);
+ assertEquals(helix1.xCenter, helix5.xCenter, 1e-10);
+ assertEquals(helix2.xCenter, helix6.xCenter, 1e-10);
+ assertEquals(helix2.yCenter, helix5.yCenter, 1e-10);
+ }
+ Helix helix1 = new Helix(origin, radius, phi, lambda);
+ Helix helix2 = new Helix(origin, -radius, phi, lambda);
+
+ SpacePoint h1 = helix1.getPointAtLength(Math.PI/2*radius/cos(lambda));
+ assertEquals(h1.x(), 4.5, 1e-14);
+ assertEquals(h1.y(), 7.5, 1e-14);
+ assertTrue(h1.z() > origin.z());
+ SpacePoint h2 = helix2.getPointAtLength(Math.PI/2*radius/cos(lambda));
+ assertEquals(h2.x(), 1.5, 1e-14);
+ assertEquals(h2.y(), 7.5, 1e-14);
+ assertTrue(h2.z() > origin.z());
+ SpacePoint h3 = helix1.getPointAtLength(-Math.PI/2*radius/cos(lambda));
+ assertEquals(h3.x(), 4.5, 1e-14);
+ assertEquals(h3.y(), 4.5, 1e-14);
+ assertTrue(h3.z() < origin.z());
+ SpacePoint h4 = helix2.getPointAtLength(-Math.PI/2*radius/cos(lambda));
+ assertEquals(h4.x(), 1.5, 1e-14);
+ assertEquals(h4.y(), 4.5, 1e-14);
+ assertTrue(h4.z() < origin.z());
+ }
+
+
+ private void assertEquals(SpacePoint v1, SpacePoint v2)
+ {
+ assertEquals(v1.x(),v2.x(), 1e-14);
+ assertEquals(v1.y(),v2.y(), 1e-14);
+ assertEquals(v1.z(),v2.z(), 1e-14);
+ }
+}