Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
WTrackUtils.java | +244 | added 1.1 | |
WTrack.java | +121 | added 1.1 | |
+365 |
Simple helix class with utils based on W representation.
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; + + } + +}
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; + } + + + + + + + + +}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1