 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
| WTrack.java | +95 | -105 | 1.6 -> 1.7 |
Reorganize for clarity
diff -u -r1.6 -r1.7 --- WTrack.java 9 Feb 2013 00:32:06 -0000 1.6 +++ WTrack.java 12 Jul 2013 20:56:21 -0000 1.7 @@ -19,76 +19,51 @@
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];
public HelicalTrackFit _htf = null;
private double _bfield;
- private double _bfield_constant; - int _q = 0;
+ private double _a;
static int max_iterations_intercept = 10;
static double epsilon_intercept = 1e-4;
- public WTrack(double [] params, double bfield, int q) {
- _bfield = bfield;
- _q = q;
- _bfield_constant = -1*Constants.fieldConversion*_bfield*_q;
- _parameters = params;
- }
-
public WTrack(WTrack trk) {
_bfield = trk._bfield;
- _q = trk.getCharge(); - _bfield_constant = trk._bfield_constant; - _parameters = trk.getParameters();
+ _a = trk._a; + _parameters = trk._parameters;
_htf = trk._htf;
_debug = trk._debug;
}
public WTrack(HelicalTrackFit track, double bfield) {
- this.initWithTrack(track, bfield, false);
+ initWithTrack(track, bfield, false);
}
public WTrack(HelicalTrackFit track, double bfield, boolean flip) {
- this.initWithTrack(track, bfield, flip);
+ initWithTrack(track, bfield, flip);
}
- public void initWithTrack(HelicalTrackFit track, double bfield, boolean flip) {
- this._htf = track;
- double signR = flip ? -1*Math.signum(track.R()): Math.signum(track.R());
- _bfield = bfield;
- _q = (int) signR;
- _bfield_constant = -1*Constants.fieldConversion*_bfield*_q;
-
- double a = _bfield_constant;
- double c = 1/(2*Math.abs(track.R())*signR); //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
+ public void initWithTrack(HelicalTrackFit track, double bfield, boolean flip) {
+ _htf = track;
+ _bfield = flip ? -1.0*bfield : bfield; // flip if needed
+ _a = -1*Constants.fieldConversion*_bfield*Math.signum(track.R());
+ double p = track.p(Math.abs(_bfield));
+ double theta = Math.PI/2.0 - Math.atan(track.slope());
+ double phi = track.phi0();
+ _parameters[0] = p*Math.cos(phi)*Math.sin(theta);
+ _parameters[1] = p*Math.sin(phi)*Math.sin(theta);
+ _parameters[2] = p*Math.cos(theta);
+ _parameters[3] = Math.sqrt(_parameters[0]*_parameters[0]+_parameters[1]*_parameters[1]+_parameters[2]*_parameters[2]);
+ _parameters[4] = -1*track.dca()*Math.sin(phi); //x0
+ _parameters[5] = track.dca()*Math.cos(phi); //y0
_parameters[6] = track.z0(); //z0
-
if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": WTrack initialized from HelicalTrackFit with params"); - System.out.println(track.toString()); - System.out.println(this.getClass().getSimpleName() + ": B-field " + _bfield + " q " + _q + " a " + a); - System.out.println(this.getClass().getSimpleName() + ": R " + track.R() + " c " + c + " phi0 " + phi0 + " pt " + pt); - double pt_tmp = Math.sqrt(_parameters[0]*_parameters[0] + _parameters[1]*_parameters[1]); - System.out.println(this.getClass().getSimpleName() + ": P0x " + _parameters[0] + " P0y " + _parameters[1] + " -> pt " + pt_tmp);
+ System.out.printf("%s: WTrack initialized (p=%f,bfield=%f,theta=%f,phi=%f) from HelicalTrackFit:\n%s:%s\n",this.getClass().getSimpleName(),
+ p,_bfield,theta,phi,
+ this.getClass().getSimpleName(),this.toString());
}
- }
- public void setTrackParameters(double [] params) {
+ }
+ public void setTrackParameters(double [] params) {
_parameters = params;
}
@@ -96,19 +71,19 @@
return _parameters;
}
- private void check() {
- if(getCharge()==0) {
- System.out.println(this.getClass().getSimpleName() + ": no track was set!");
- System.exit(1);
- }
+ private boolean goingForward() {
+ // assuming the track should go in the x-direction --> not very general -> FIX THIS!?
+ return getP0().x()>0 ? true : false;
}
- public double getBfieldConstant() {
- return this._bfield_constant;
+
+ public double a() {
+ return _a;
+
}
private int getCharge() {
- return _q;
+ return (int) Math.signum(_htf.R());
}
public Hep3Vector getP0() {
@@ -139,8 +114,9 @@
- private Hep3Vector getMomentumOnHelix(WTrack track, double s) {
- double a = track.getBfieldConstant();
+ private Hep3Vector getMomentumOnHelix(double s) {
+ WTrack track = this;
+ double a = track.a();
Hep3Vector p0 = track.getP0();
double rho = a / p0.magnitude();
double px = p0.x()*Math.cos(rho*s) - p0.y()*Math.sin(rho*s);
@@ -149,8 +125,9 @@
return (new BasicHep3Vector(px,py,pz));
}
- private Hep3Vector getPointOnHelix(WTrack track, double s) {
- double a = track.getBfieldConstant();
+ private Hep3Vector getPointOnHelix(double s) {
+ WTrack track = this;
+ double a = track.a();
Hep3Vector p0 = track.getP0();
Hep3Vector x0 = track.getX0();
double rho = a / p0.magnitude();
@@ -160,13 +137,13 @@
return (new BasicHep3Vector(x,y,z));
}
- private double getPathLengthToPlaneApprox(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
+ private double getPathLengthToPlaneApprox(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
*/
- - double a = track.getBfieldConstant();
+ WTrack track = this; + double a = track.a();
Hep3Vector p0 = track.getP0();
Hep3Vector x0 = track.getX0();
double p = p0.magnitude();
@@ -177,8 +154,8 @@
double t = B*B-4*A*C;
if(t<0) {
System.out.println(" getPathLengthToPlaneApprox ERROR t is negative (" + t + ")!" );
- System.out.println(" p " + p + " rho " + rho + " A " + A + " B " + B + " C " + C);
- System.out.println(" track params\n" + track.toString());
+ System.out.println(" p " + p + " rho " + rho + " a " + a + " A " + A + " B " + B + " C " + C);
+ System.out.println(" track params: " + track.paramsToString());
System.out.println(" xp " + xp.toString());
System.out.println(" eta " + eta.toString());
System.out.println(" h " + h.toString());
@@ -201,7 +178,7 @@
// }
if(_debug) {
System.out.println(" getPathLengthToPlaneApprox ");
- System.out.println(" " + track.toString());
+ System.out.println(" " + track.paramsToString());
System.out.println(" xp " + xp.toString());
System.out.println(" eta " + eta.toString());
System.out.println(" h " + h.toString());
@@ -213,13 +190,13 @@
}
- private Hep3Vector getPointOnHelix(WTrack track, double s, Hep3Vector h) {
+ private Hep3Vector getPointOnHelix(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();
+ WTrack track = this; + double a = track.a();
Hep3Vector p0 = track.getP0();
double p = p0.magnitude();
Hep3Vector x0 = track.getX0();
@@ -235,13 +212,13 @@
return x;
}
- private Hep3Vector getMomentumOnHelix(WTrack track, double s, Hep3Vector h) {
+ private Hep3Vector getMomentumOnHelix(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();
+ WTrack track = this; + double a = track.a();
Hep3Vector p0 = track.getP0();
double rho = a / p0.magnitude();
double srho = s*rho;
@@ -254,7 +231,7 @@
}
- private double[] getHelixParametersAtPathLength(WTrack track, double s, Hep3Vector h) {
+ private double[] getHelixParametersAtPathLength(double s, Hep3Vector h) {
/*
* Calculate the exact position of the new helix parameters at path length s
* in an arbitrarily oriented, constant magnetic field
@@ -264,11 +241,11 @@
*/
// Find track parameters at that path length
- Hep3Vector p = getMomentumOnHelix(track, s, h); - Hep3Vector x = getPointOnHelix(track,s, h);
+ Hep3Vector p = getMomentumOnHelix(s, h); + Hep3Vector x = getPointOnHelix(s, h);
- Hep3Vector p_tmp = getMomentumOnHelix(track, s); - Hep3Vector x_tmp = getPointOnHelix(track, s);
+ Hep3Vector p_tmp = getMomentumOnHelix(s); + Hep3Vector x_tmp = getPointOnHelix(s);
if(_debug) {
System.out.println(" point on helix at s");
@@ -282,7 +259,7 @@
pars[0] = p.x();
pars[1] = p.y();
pars[2] = p.z();
- pars[3] = track.getParameters()[3]; //E is unchanged
+ pars[3] = getParameters()[3]; //E is unchanged
pars[4] = x.x();
pars[5] = x.y();
pars[6] = x.z();
@@ -291,9 +268,6 @@
public Hep3Vector getHelixAndPlaneIntercept(Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
- return this.getHelixAndPlaneIntercept(this, xp, eta, h);
- }
- public Hep3Vector getHelixAndPlaneIntercept(WTrack track, Hep3Vector xp, Hep3Vector eta, Hep3Vector h) {
/*
* Find the interception point between the helix and plane
@@ -304,40 +278,56 @@
int iteration = 1;
- double delta_s = 9999; - double s_prev = 999999;
+ double s_total = 0.; + double step = 9999999.9;
//List<WTrack> tracks = new ArrayList<WTrack>();
- WTrack trk = new WTrack(track);
- while(iteration<=max_iterations_intercept && Math.abs(s_prev)>epsilon_intercept) {
- if(_debug) System.out.printf("%s: Iteration %d: s_pref =%.3f delta_s=%.3f\ntrk params: [%s] \n",this.getClass().getSimpleName(),iteration,s_prev,delta_s,paramsToString());
-
- // Start by approximating the path length to the point
- double s = getPathLengthToPlaneApprox(trk, xp, eta, h);
-
- if(_debug) System.out.printf("%s: Iteration %d: path length step s=%.3f\n",this.getClass().getSimpleName(),iteration,s);
+ WTrack trk = this;
+ while(iteration<=max_iterations_intercept && Math.abs(step)>epsilon_intercept) {
- // Find the track parameters at this point - double[] params = getHelixParametersAtPathLength(trk, s, h); - // update the track parameters - trk.setTrackParameters(params);
+ if(_debug) {
+ System.out.printf("%s: Iteration %d\n", this.getClass().getSimpleName(),iteration);
+ System.out.printf("%s: s_total %f prev_step %.3f current trk params: %s \n",
+ this.getClass().getSimpleName(),s_total,step,trk.paramsToString());
+ }
- if(_debug) System.out.printf("%s: Iteration %d: updated trak params: [%s]\n",this.getClass().getSimpleName(),iteration,trk.paramsToString());
+ // check that the track is not looping
- //tracks.add(trk); - iteration++; - delta_s = Math.abs(s-s_prev); - s_prev = s;
+ if(trk.goingForward()) {
+
+
+ // Start by approximating the path length to the point
+ step = getPathLengthToPlaneApprox(xp, eta, h);
+
+ if(_debug) System.out.printf("%s: path length step s=%.3f\n",this.getClass().getSimpleName(),step);
+
+ // Find the track parameters at this point
+ double[] params = getHelixParametersAtPathLength(step, h);
+
+ // update the track parameters
+ trk.setTrackParameters(params);
+
+ if(_debug) System.out.printf("%s: updated track params: [%s]\n",this.getClass().getSimpleName(),trk.paramsToString());
+
+ //tracks.add(trk);
+ iteration++;
+ s_total += step;
+
+ //Save distance between point and estimate
+ //Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
- //Save distance between point and estimate - //Hep3Vector dpoint = VecOp.sub(xp, trk.getX0());
+ } else {
+ //if(_debug)
+ System.out.printf("%s: this track started to go backwards?! params [%s]\n",this.getClass().getSimpleName(),trk.toString());
+ return null;
+ }
}
- if(_debug) System.out.printf("%s: final s=%.3f and ds=%.3f after %d iterations with track params: [%s]\n",this.getClass().getSimpleName(),s_prev,delta_s,iteration,trk.paramsToString());
+ if(_debug) System.out.printf("%s: final total_s=%f with final step %f after %d iterations gave track params: %s\n",
+ this.getClass().getSimpleName(),s_total,step,iteration,trk.paramsToString());
- Hep3Vector point = trk.getX0(); - return point;
+ return trk.getX0();
}
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