lcsim/src/org/lcsim/recon/vertexing/zvtop4
diff -N ZvSwimmer.java
--- ZvSwimmer.java 23 Aug 2005 05:41:08 -0000 1.13
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,274 +0,0 @@
-/**
- *
- */
-package org.lcsim.recon.vertexing.zvtop4;
-
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-
-import org.lcsim.event.Track;
-import org.lcsim.geometry.Detector;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.util.swim.HelixSwim;
-
-import Jama.Matrix;
-
-import static org.lcsim.recon.vertexing.zvtop4.ZvUtil.getUnitTangent;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.cross;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.unit;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.dot;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.subtract;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.multiply;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.add;
-import static org.lcsim.recon.vertexing.zvtop4.VectorArithmetic.magnitude;
-import static org.lcsim.recon.vertexing.zvtop4.ZvUtil.fieldConversion;
-
-/**
- * @author jstrube
- * @version
- */
-public class ZvSwimmer extends HelixSwim {
-
- private static double maxSinLambda = 0.99;
- private static double stepMin = 0.005;
- private static int nMaxTry = 7;
- private static double stepSize = 0.2;
-
-
- static Detector detector;
- /**
- * @param B
- */
- // FIXME all the constructors can leave the class in an inconsistent state
- public ZvSwimmer(double B) {
- super(B);
- // TODO Auto-generated constructor stub
- }
-
- /**
- * @param B
- * @param rhoMin
- * @param rhoMax
- * @param zMax
- */
- public ZvSwimmer(double B, double rhoMin, double rhoMax, double zMax, Detector det) {
- super(B, rhoMin, rhoMax, zMax);
- // TODO Auto-generated constructor stub
- detector = det;
- }
-
- /**
- * @param B
- * @param p
- * @param r0
- * @param iq
- */
- public ZvSwimmer(double B, double[] p, double[] r0, int iq) {
- super(B, p, r0, iq);
- // TODO Auto-generated constructor stub
- }
-
-
- /**
- * Swims the track to a POCA at a line parallel to the target position.
- * @param track the Track to swim
- * @param targetPosition the point in space where to swim the track
- * @return a ZvFitStatus object to determine if the action was successful
- */
- public ZvFitStatus swimTo(Track track, SpacePoint targetPosition) {
- Hep3Vector orgBField = new BasicHep3Vector(detector.getFieldMap().getField(track.getReferencePoint()));
- Hep3Vector orgTangent = getUnitTangent(track.getTrackParameters());
- Hep3Vector orgNormal = cross(orgTangent, orgBField);
- ZvFitStatus status = ZvFitStatus.FAILED_ON_ENTRY;
- Track.Parameters orgParameters = new Track.Parameters(track.getTrackParameters());
- Matrix orgErrorMatrix = new Matrix(track.getErrorMatrix());
-
- // FIXME detector radius
- double radius = 100;
- // set and check if target position is outside the detector
- if (targetPosition.rxyz()> radius)
- return ZvFitStatus.POS_OUT_OF_VOLUME;
-
- // get track parameters and error matrix
- double[] hlxPar = track.getTrackParameters();
- double[][] dHlxPar = track.getErrorMatrix();
-
- double charge = track.getCharge();
- double sgnq = charge / Math.abs(charge);
-
- // FIXME cast from double to BasicHep3Vector
- double pTot = new BasicHep3Vector(track.getMomentum()).magnitude();
-
- // TODO
-// ZvHelixVec hlxVec = new ZvHelixVec(hlxPar);
-
- // current position vector
- // FIXME
- SpacePoint vCurPos;
- {
- double[] x = track.getReferencePoint();
- vCurPos = new CartesianPoint(x[0], x[1], x[2]);
- }
- SpacePoint orgPosition = new SpacePoint(vCurPos);
-
- Hep3Vector tangentUnit = getUnitTangent(hlxPar);
- Hep3Vector magneticFieldUnit = new BasicHep3Vector(detector.getFieldMap().getField(vCurPos.getCartesianArray()));
- Hep3Vector normalUnit = unit(cross(tangentUnit, magneticFieldUnit));
-
- Hep3Vector vbdUnit = new BasicHep3Vector(0, 0, 1);
-
-
- // check for cos(Lambda)
- if (Math.abs(tangentUnit.z()) > maxSinLambda)
- return ZvFitStatus.TOO_BIG_DIP_ANGLE;
-
- // master loop
- double f = 0.;
- double sumArc = 0.;
- int iter = 0;
-
- do {
- // counter
- if (++iter > nMaxTry)
- return ZvFitStatus.TOO_MANY_ITERATIONS;
-
- // now calculate travel distance along track
- double pb = dot(tangentUnit, vbdUnit);
- Hep3Vector vu = subtract(targetPosition, vCurPos);
-
- // distance along p direction to pt of closest approach
- f = (dot(vu, tangentUnit) - pb * dot(vu, vbdUnit)) / (1. - pb * pb);
-
- // preparation for stepping loop
- double fLeft = f;
- double fSign = f / Math.abs(f);
- while (Math.abs(fLeft) > 0.) {
- double arcLen;
- if (Math.abs(fLeft) > stepSize)
- arcLen = fSign * stepSize;
- else
- arcLen = fLeft;
-
- // now do a step
-
- vCurPos = leapToNewPosition(track, arcLen, vCurPos, tangentUnit, pTot, sgnq);
- tangentUnit = globalpUnit;
- // check for containment
- if (Math.abs(globalpUnit.z()) > maxSinLambda)
- return ZvFitStatus.TOO_BIG_DIP_ANGLE;
-
- // final accounting
- sumArc += arcLen;
- fLeft -= arcLen;
- } // stepping loop
- } while (Math.abs(f) > stepMin); // master loop
-
- // recalculate helix parameters
-// int phi0 = Track.ParameterName.phi0.ordinal();
-// int tanLambda = Track.ParameterName.s.ordinal();
-// hlxPar[phi0] = Math.atan2(tangentUnit.y(), tangentUnit.x());
-// if (hlxPar[phi0] < 0)
-// hlxPar[phi0] += 2. * Math.PI; // phi
-// // FIXME these parameters are screwed up for sure, ZvSwimmer is to be obsoleted anyway
- double cosL = Math.sqrt(1. - tangentUnit.z() * tangentUnit.z());
-// hlxPar[2] = 1. / (pTot * cosL); // kappa
-// hlxPar[tanLambda] = tangentUnit.z() / cosL; // tan(lambda)
-// hlxPar[4] = vCurPos.x(); // x
-// hlxPar[5] = vCurPos.y(); // y
-// hlxPar[6] = vCurPos.z(); // z
-
- // check current position
- if ( vCurPos.rxy() > maxRadius || vCurPos.z() > maxZ )
- return ZvFitStatus.POS_OUT_OF_VOLUME;
-
- // calculate doca reference point on target line
- Hep3Vector vu = subtract(vCurPos, targetPosition);
- double lambda = dot(vu, tangentUnit) / dot(vbdUnit, tangentUnit);
- SpacePoint vxRef = add(targetPosition, multiply(vbdUnit, lambda));
-
- // now swim the error matrix
- if (Math.abs(sumArc) > 0.) {
- double sign = -sumArc / Math.abs(sumArc);
- double fsum = Math.abs(sumArc);
- Hep3Vector dVec = subtract(vCurPos, orgPosition);
- double dis2 = Math.sqrt(dVec.x() * dVec.x() + dVec.y() * dVec.y());
-
- double bFieldStrength = magnitude(detector.getFieldMap().getField(vCurPos.getCartesianArray()));
- magneticFieldUnit = unit(detector.getFieldMap().getField(vCurPos.getCartesianArray()));
-
- normalUnit = unit(cross(tangentUnit, magneticFieldUnit));
-
- // build the transformation matrix
- Matrix A = new Matrix(5, 5, 0.);
- A.set(0, 0, 1.);
- A.set(0, 1, fsum * fieldConversion * bFieldStrength * cosL * sgnq * sign);
- A.set(1, 1, 1.);
- A.set(2, 2, 1.);
- A.set(3, 0, dVec.x() * normalUnit.y() - dVec.y() * normalUnit.x());
- A.set(3, 1, sgnq * 0.5 * dis2 * fieldConversion * bFieldStrength
- * (tangentUnit.y() * normalUnit.x() - tangentUnit.x() * normalUnit.y()) / cosL);
- A.set(3, 2, -sign * fsum * cosL * normalUnit.z());
- A.set(3, 3, dot(normalUnit, orgNormal));
- A.set(3, 4, dot(normalUnit, orgBField));
- A.set(4, 0, dVec.x() * magneticFieldUnit.y() - dVec.y() * magneticFieldUnit.x());
- A.set(4, 1, sgnq * 0.5 * dis2 * fieldConversion * bFieldStrength
- * (tangentUnit.y() * magneticFieldUnit.x() - tangentUnit.x() * magneticFieldUnit.y()) / cosL);
- A.set(4, 2, -sign * fsum * cosL * magneticFieldUnit.z());
- A.set(4, 3, dot(magneticFieldUnit, orgNormal));
- A.set(4, 4, dot(magneticFieldUnit, orgBField));
-
- dHlxPar = A.times(orgErrorMatrix.times(A.transpose())).getArray();
- }
-
- // record swum distance
-// status.setDist(sumArc);
-//
-// // set up the resulting track
-// swumTrack = track.copy();
-// swumTrack.setZvTrackParam(hlxPar);
-// swumTrack.setZvTrackError(dHlxPar);
-// swumTrack.setSwimStatus(new ZvSwimStatus(status));
- return ZvFitStatus.OK;
- }
- // FIXME bad example of global variable
- static Hep3Vector globalpUnit;
- // TODO naming of vectors is obscure
- private static SpacePoint leapToNewPosition(Track track
- , double arcLen
- , SpacePoint xPosition
- , Hep3Vector tan
- , double pTot
- , double charge) {
- Hep3Vector oldTangent = getUnitTangent(track.getTrackParameters());
-
- Hep3Vector tangentUnit = unit(tan);
- Hep3Vector bField = unit(detector.getFieldMap().getField(xPosition.getCartesianArray()));
- // is already unit vec
- Hep3Vector parallelUnit = cross(tangentUnit, bField);
- Hep3Vector normalFieldVector = unit(cross(bField, parallelUnit));
- double bFieldMag = bField.magnitude();
- double cosPit = dot(tangentUnit, normalFieldVector);
- double sinPit = dot(tangentUnit, bField);
- double tanPit = sinPit / cosPit;
- double omega = charge/Math.abs(charge)*bFieldMag/ ( 333.567 * pTot*cosPit );
-
- double alpha = omega * cosPit * arcLen;
-
- Hep3Vector delPos = multiply(add(multiply(normalFieldVector, Math.sin(alpha)),
- add(multiply(parallelUnit, 1.-Math.cos(alpha)), multiply(bField, alpha*tanPit))), 1./omega);
-
- globalpUnit = multiply(add(multiply(normalFieldVector, Math.cos(alpha)),
- add(multiply(parallelUnit, Math.sin(alpha)), multiply(bField, tanPit))), cosPit);
-
- return add(xPosition, delPos);
- }
-
- public void setDetector(Detector d) {
- detector = d;
- }
-
- private static double maxRadius = 999999.; // cm
- private static double maxZ = 999999.; // cm
-}
