hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N WTrackUtils.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ WTrackUtils.java 9 Oct 2012 01:20:02 -0000 1.1
@@ -0,0 +1,244 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.users.phansson;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.List;
+
+/**
+ *
+ * @author phansson
+ */
+public class WTrackUtils {
+ boolean _debug = false;
+ List<Double> _delta_point = new ArrayList<Double>();
+
+ public WTrackUtils() {
+ }
+
+ public WTrackUtils(boolean debug) {
+ _debug = debug;
+ }
+
+ public void setDebug(boolean debug) {
+ _debug = debug;
+ }
+
+ public List<Double> getDeltas() {
+ return _delta_point;
+ }
+
+ private Hep3Vector getMomentumOnHelix(WTrack track, double s) {
+ double a = track.getBfieldConstant();
+ Hep3Vector p0 = track.getP0();
+ double rho = a / p0.magnitude();
+ double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s);
+ double py = p0.y()*Math.cos(rho*s) + p0.x()*Math.sin(rho*s);
+ double pz = p0.z();
+ return (new BasicHep3Vector(px,py,pz));
+ }
+
+ private Hep3Vector getPointOnHelix(WTrack track, double s) {
+ double a = track.getBfieldConstant();
+ Hep3Vector p0 = track.getP0();
+ Hep3Vector x0 = track.getX0();
+ double rho = a / p0.magnitude();
+ double x = x0.x() + p0.x()/a*Math.sin(rho*s) - p0.y()/a*(1-Math.cos(rho*s));
+ double y = x0.y() + p0.y()/a*Math.sin(rho*s) + p0.x()/a*(1-Math.cos(rho*s));
+ double z = x0.z() + p0.z()/p0.magnitude()*s;
+ return (new BasicHep3Vector(x,y,z));
+ }
+
+ public double getPathLengthToPlaneApprox(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+ /*
+ * Find the approximate path length to the point xp
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": getPathLengthToPlaneApprox ");
+ System.out.println(this.getClass().getSimpleName() + ": " + track.toString());
+ System.out.println(this.getClass().getSimpleName() + ": xp " + xp.toString());
+ System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
+ System.out.println(this.getClass().getSimpleName() + ": h " + h.toString());
+ }
+ double a = track.getBfieldConstant();
+ Hep3Vector p0 = track.getP0();
+ Hep3Vector x0 = track.getX0();
+ double p = p0.magnitude();
+ double rho = a / p0.magnitude();
+ double A = VecOp.dot(eta,VecOp.cross(p0, h))/p*0.5*rho;
+ double B = VecOp.dot(p0,eta)/p;
+ double C = VecOp.dot(VecOp.sub(x0,xp),eta);
+ double root1 = (-B + Math.sqrt(B*B-4*A*C)) /(2*A);
+ double root2 = (-B - Math.sqrt(B*B-4*A*C)) /(2*A);
+ double root = -99999999;
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": p " + p + " rho " + rho + " A " + A + " B " + B + " C " + C);
+ }
+ // choose the smallest positive solution
+ // if both negative choose the smallest negative ???
+ //if(root1==0 || root2==0) root=0;
+ root = Math.abs(root1) <= Math.abs(root2) ? root1 : root2;
+// else if(Math.signum(root1)>0 && Math.signum(root2)<0) root = root1;
+// else if(Math.signum(root2)>0 && Math.signum(root1)<0) root = root2;
+// else if(Math.signum(root1)>0 && Math.signum(root2)>0) root = root1 > root2 ? root2 : root1;
+// else if(Math.signum(root1)<0 && Math.signum(root2)<0) root = root1 < root2 ? root2 : root1;
+// else {
+// System.out.println(this.getClass().getSimpleName() + ": I should never get here! (root1=" + root1 + " root2=" + root2+")");
+// System.exit(1);
+// }
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": root1 " + root1 + " root2 " + root2 + " -> root " + root);
+ }
+ return root;
+
+ }
+
+
+ private Hep3Vector getPointOnHelix(WTrack track, double s, Hep3Vector h) {
+ /*
+ * Get point on helix at path lenght s
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+
+ double a = track.getBfieldConstant();
+ Hep3Vector p0 = track.getP0();
+ double p = p0.magnitude();
+ Hep3Vector x0 = track.getX0();
+ double rho = a / p0.magnitude();
+ double srho = s*rho;
+ Hep3Vector a1 = VecOp.mult(1/a*Math.sin(srho), p0);
+ Hep3Vector a2 = VecOp.mult(1/a*(1-Math.cos(srho)),VecOp.cross(p0,h));
+ Hep3Vector a3 = VecOp.mult(VecOp.dot(p0, h)/p,h);
+ Hep3Vector a4 = VecOp.mult(s-Math.sin(srho)/rho, a3);
+ Hep3Vector x = VecOp.add(x0,a1);
+ x = VecOp.sub(x,a2);
+ x = VecOp.add(x,a4);
+ return x;
+ }
+
+ private Hep3Vector getMomentumOnHelix(WTrack track, double s, Hep3Vector h) {
+ /*
+ * Get point on helix at path lenght s
+ * in arbitrary oriented, constant magnetic field with unit vector h
+ */
+
+ double a = track.getBfieldConstant();
+ Hep3Vector p0 = track.getP0();
+ double rho = a / p0.magnitude();
+ double srho = s*rho;
+ Hep3Vector a1 = VecOp.mult(Math.cos(srho), p0);
+ Hep3Vector a2 = VecOp.cross(p0, VecOp.mult(Math.sin(srho),h));
+ Hep3Vector a3 = VecOp.mult(VecOp.dot(p0,h),VecOp.mult(1-Math.cos(srho),h));
+ Hep3Vector p = VecOp.sub(a1,a2);
+ p = VecOp.add(p,a3);
+ return p;
+ }
+
+
+ public double[] getHelixParametersAtPathLength(WTrack track, double s, Hep3Vector h) {
+ /*
+ * Calculate the exact position of the new helix parameters at path length s
+ * in an arbitrarily oriented, constant magnetic field
+ *
+ * point xp is the point
+ * h is a unit vector in the direction of the magnetic field
+ */
+
+ // Find track parameters at that path length
+ Hep3Vector p = this.getMomentumOnHelix(track, s, h);
+ Hep3Vector x = this.getPointOnHelix(track, s, h);
+
+ Hep3Vector p_tmp = this.getMomentumOnHelix(track, s);
+ Hep3Vector x_tmp = this.getPointOnHelix(track, s);
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": point on helix at s");
+ System.out.println(this.getClass().getSimpleName() + ": p " + p.toString() + " p_tmp " + p_tmp.toString());
+ System.out.println(this.getClass().getSimpleName() + ": x " + x.toString() + " x_tmp " + x_tmp.toString());
+ }
+
+
+ //Create the new parameter array
+ double [] pars = new double[7];
+ pars[0] = p.x();
+ pars[1] = p.y();
+ pars[2] = p.z();
+ pars[3] = track.getParameters()[3]; //E is unchanged
+ pars[4] = x.x();
+ pars[5] = x.y();
+ pars[6] = x.z();
+ return pars;
+ }
+
+
+
+ public Hep3Vector getHelixAndPlaneIntercept(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+
+ /*
+ * Find the interception point between the helix and plane
+ * track: the helix in W representation
+ * xp: point on the plane
+ * eta: unit vector of the plane
+ * h: unit vector of magnetic field
+ */
+
+ _delta_point.clear();
+
+ int max_it = 10;
+ double epsilon = 1e-4;
+ int iteration = 1;
+ double delta_s = 9999;
+ double s_prev = 999999;
+ List<WTrack> tracks = new ArrayList<WTrack>();
+ WTrack trk = new WTrack(track);
+ while(iteration<=max_it && Math.abs(s_prev)>epsilon) {
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": iteration " + iteration + " --- ");
+ System.out.println(this.getClass().getSimpleName() + ": s_prev " + s_prev + " delta_s " + delta_s);
+ System.out.println(this.getClass().getSimpleName() + ": " + trk.toString());
+ }
+ // Start by approximating the path length to the point
+ double s = this.getPathLengthToPlaneApprox(trk, xp, eta, h);
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": Approx. path s " + s);
+ }
+
+ // Find the track parameters at this point
+ double[] params = this.getHelixParametersAtPathLength(trk, s, h);
+ // update the track parameters
+ trk.setTrackParameters(params);
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": updated params " + trk.toString());
+ }
+
+ tracks.add(trk);
+ iteration++;
+ delta_s = Math.abs(s-s_prev);
+ s_prev = s;
+
+ //Save distance between point and estimate
+ Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
+ _delta_point.add(Math.abs(VecOp.dot(eta, dpoint)));
+
+
+ }
+
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": Final s " + s_prev + " delta_s " + delta_s);
+ System.out.println(this.getClass().getSimpleName() + ": Final track params " + trk.toString());
+ }
+
+ Hep3Vector point = trk.getX0();
+ return point;
+
+ }
+
+}
hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N WTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ WTrack.java 9 Oct 2012 01:20:02 -0000 1.1
@@ -0,0 +1,121 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.users.phansson;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.constants.Constants;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+
+/**
+ *
+ * @author phansson
+ */
+public class WTrack {
+
+ private boolean _debug = false;
+ //public enum PARAM{P0X(0),P0Y(1),P0Z(2),E(3),X0(4),Y0(5),Z0(6);}
+ public enum PARAM{TEST;}
+ private double[] _parameters = new double[7];
+ //HelicalTrackFit _canonical_track = null;
+ double _bfield;
+ double _bfield_constant;
+
+ int _q = 0;
+
+ public WTrack(double [] params, double bfield, int q) {
+ _bfield = bfield;
+ _q = q;
+ _bfield_constant = Constants.fieldConversion*_bfield*_q;
+ _parameters = params;
+ }
+
+ public WTrack(WTrack trk) {
+ _bfield = trk._bfield;
+ _q = trk.getCharge();
+ _bfield_constant = Constants.fieldConversion*_bfield*_q;
+ _parameters = trk.getParameters();
+ }
+
+
+ public WTrack(HelicalTrackFit track, double bfield) {
+ _bfield = bfield;
+ _q = (int) Math.signum(track.R());
+ _bfield_constant = Constants.fieldConversion*_bfield*_q;
+
+ double a = _bfield_constant;
+ double c = 1/(2*track.R()); //note the different definition than in canonical track
+ double phi0 = track.phi0();
+ double lambda = track.slope(); // lambda = 1/tan(theta)
+ double pt = track.pT(_bfield);
+ double m = 0.0;
+ double E0 = Math.sqrt(pt*pt*(1+lambda*lambda) + m*m);
+ double x0 = -1*track.dca()*Math.sin(phi0);
+ double y0 = track.dca()*Math.cos(phi0);
+ double z0 = track.z0();
+
+ _parameters[0] = a/(2*c)*Math.cos(phi0); //P0x
+ _parameters[1] = a/(2*c)*Math.sin(phi0); //P0y
+ _parameters[2] = a/(2*c)*lambda; //P0z
+ _parameters[3] = Math.sqrt(pt*pt*(1+lambda*lambda) + m*m); //E0
+ _parameters[4] = -1*track.dca()*Math.sin(phi0); //x0
+ _parameters[5] = track.dca()*Math.cos(phi0); //y0
+ _parameters[6] = track.z0(); //z0
+
+ double pt_tmp = Math.sqrt(_parameters[0]*_parameters[0] + _parameters[1]*_parameters[1]);
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": a " + a + " R " + track.R() + " c " + c + " phi0 " + phi0 + " pt " + pt);
+ System.out.println(this.getClass().getSimpleName() + ": P0x " + _parameters[0] + " P0y " + _parameters[1] + " -> pt " + pt_tmp);
+ }
+ }
+
+ public void setTrackParameters(double [] params) {
+ _parameters = params;
+ }
+
+ public double[] getParameters() {
+ return _parameters;
+ }
+
+ private void check() {
+ if(getCharge()==0) {
+ System.out.println(this.getClass().getSimpleName() + ": no track was set!");
+ System.exit(1);
+ }
+ }
+
+ public double getBfieldConstant() {
+ return this._bfield_constant;
+ }
+
+ private int getCharge() {
+ return _q;
+ }
+
+ public Hep3Vector getP0() {
+ return ( new BasicHep3Vector(_parameters[0],_parameters[1],_parameters[2]));
+ }
+
+ public Hep3Vector getX0() {
+ return ( new BasicHep3Vector(_parameters[4],_parameters[5],_parameters[6]));
+ }
+
+ public String toString() {
+
+ String str = "WTrack params [";
+ for(int i=0;i<7;++i) str += _parameters[i] + ", ";
+ str += "]";
+ return str;
+ }
+
+
+
+
+
+
+
+
+}