Commit in hps-java/src/main/java/org/lcsim/hps on MAIN | |||
recon/tracking/HelicalTrackHitResidualsDriver.java | +21 | -18 | 1.1 -> 1.2 |
/DataTrackerFakeHitDriver.java | +7 | -30 | 1.6 -> 1.7 |
/TrackUtils.java | +41 | -94 | 1.11 -> 1.12 |
users/phansson/StripMPAlignmentInput.java | +15 | -25 | 1.1 -> 1.2 |
/MPAlignmentInputCalculator.java | -2 | 1.1 -> 1.2 | |
/WTrack.java | +223 | -8 | 1.4 -> 1.5 |
/ModuleMPAlignmentInput.java | +2 | -1 | 1.1 -> 1.2 |
+309 | -178 |
Cleaning up utils for wtrack, improving usability and reusability by making helixintercept and a few other methods static.
diff -u -r1.1 -r1.2 --- HelicalTrackHitResidualsDriver.java 19 Nov 2012 21:53:42 -0000 1.1 +++ HelicalTrackHitResidualsDriver.java 27 Nov 2012 02:17:29 -0000 1.2 @@ -34,32 +34,32 @@
public class HelicalTrackHitResidualsDriver extends Driver { private AIDA aida = AIDA.defaultInstance();
- int totalTracks=0; - int totalTracksProcessed=0;
+ private int totalTracks=0; + private int totalTracksProcessed=0;
private String outputPlotFileName=""; private String trackCollectionName="MatchedTracks"; private boolean hideFrame = false; private boolean _debug = false;
- private TrackUtils _trackUtils = new TrackUtils();
+ private boolean _includeMS = true;
private AIDAFrame plotterFrame; // ICloud1D[] _h_resz_track_top = new ICloud1D[5]; // ICloud1D[] _h_resz_track_bottom = new ICloud1D[5]; // ICloud1D[] _h_resy_track_top = new ICloud1D[5]; // ICloud1D[] _h_resy_track_bottom = new ICloud1D[5];
- IHistogram1D[] _h_resz_track_top = new IHistogram1D[5]; - IHistogram1D[] _h_resz_track_bottom = new IHistogram1D[5]; - IHistogram1D[] _h_resy_track_top = new IHistogram1D[5]; - IHistogram1D[] _h_resy_track_bottom = new IHistogram1D[5]; - IDataPointSet dps_hth_y_b; - IDataPointSet dps_hth_y_t; - IDataPointSet dps_hth_z_b; - IDataPointSet dps_hth_z_t; - IPlotter _plotter_resz_top; - IPlotter _plotter_resy_top; - IPlotter _plotter_resz_bottom; - IPlotter _plotter_resy_bottom; - IPlotter _plotter_mean_res;
+ private IHistogram1D[] _h_resz_track_top = new IHistogram1D[5]; + private IHistogram1D[] _h_resz_track_bottom = new IHistogram1D[5]; + private IHistogram1D[] _h_resy_track_top = new IHistogram1D[5]; + private IHistogram1D[] _h_resy_track_bottom = new IHistogram1D[5]; + private IDataPointSet dps_hth_y_b; + private IDataPointSet dps_hth_y_t; + private IDataPointSet dps_hth_z_b; + private IDataPointSet dps_hth_z_t; + private IPlotter _plotter_resz_top; + private IPlotter _plotter_resy_top; + private IPlotter _plotter_resz_bottom; + private IPlotter _plotter_resy_bottom; + private IPlotter _plotter_mean_res;
@@ -80,6 +80,10 @@
hideFrame = hide; }
+ public void setIncludeMS(boolean inc) { + this._includeMS = inc; + } +
public HelicalTrackHitResidualsDriver() { }
@@ -118,8 +122,7 @@
List<TrackerHit> hitsOnTrack = track.getTrackerHits(); for(TrackerHit hit : hitsOnTrack) { HelicalTrackHit hth = (HelicalTrackHit) hit;
- Map<String,Double> res_track = _trackUtils.calculateTrackHitResidual(hth, trk, false); - //Map<String,Double> res_trackgen = _trackUtils.calculateTrackHitResidual(hth, trkgen, false);
+ Map<String,Double> res_track = TrackUtils.calculateTrackHitResidual(hth, trk, this._includeMS);
boolean isTop = false; if(SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)hth.getRawHits().get(0)).getDetectorElement())) { isTop = true;
diff -u -r1.6 -r1.7 --- DataTrackerFakeHitDriver.java 19 Nov 2012 22:11:11 -0000 1.6 +++ DataTrackerFakeHitDriver.java 27 Nov 2012 02:17:29 -0000 1.7 @@ -10,11 +10,7 @@
import hep.physics.vec.Hep3Vector; import hep.physics.vec.VecOp; import java.util.*;
-import org.lcsim.detector.IDetectorElement; -import org.lcsim.detector.ITransform3D; -import org.lcsim.detector.ITranslation3D; -import org.lcsim.detector.Transform3D; -import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.*;
import org.lcsim.detector.identifier.IIdentifier; import org.lcsim.detector.solids.Box; import org.lcsim.detector.solids.Inside;
@@ -33,7 +29,6 @@
import org.lcsim.geometry.Detector; import org.lcsim.geometry.subdetector.BarrelEndcapFlag; import org.lcsim.hps.users.phansson.WTrack;
-import org.lcsim.hps.users.phansson.WTrackUtils;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit; import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D; import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
@@ -53,7 +48,6 @@
Hep3Matrix detToTrk; SvtTrackExtrapolator extrapolator = new SvtTrackExtrapolator(); Hep3Vector _bfield;
- WTrackUtils wutils = new WTrackUtils();
TrackerHitType trackerType = new TrackerHitType(TrackerHitType.CoordinateSystem.GLOBAL, TrackerHitType.MeasurementType.STRIP_1D); CoordinateSystem coordinate_system = trackerType.getCoordinateSystem(); private HitIdentifier _ID = new HitIdentifier();
@@ -161,7 +155,6 @@
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.); _bfield = new BasicHep3Vector(0,0,detector.getFieldMap().getField(IP).y()); detToTrk = trackerhitutils.detToTrackRotationMatrix();
- this.wutils.setDebug(this.debug);
// Process detectors specified by path, otherwise process entire
@@ -419,7 +412,7 @@
if(debug) System.out.println(this.getClass().getSimpleName() + ": eta " + eta.toString());
- Hep3Vector position = wutils.getHelixAndPlaneIntercept(wtrack, _orgloc, eta, h);
+ Hep3Vector position = wtrack.getHelixAndPlaneIntercept(_orgloc, eta, h);
if(debug) System.out.println(this.getClass().getSimpleName() + ": found interception point at position " + position.toString());
@@ -455,15 +448,7 @@
if(debug) System.out.println(this.getClass().getSimpleName() + ": fixed global hit position " + position.toString());
- //Fill some debug stuff - List<Double> deltas = wutils.getDeltas(); - this._delta_itercount.get(sensor).fill(deltas.size()); - IProfile1D dhisto = this._delta_histos.get(sensor); - for(int i=0;i<deltas.size();++i) { - this._prf_all_deltas.fill(i, deltas.get(i)); - dhisto.fill(i,deltas.get(i)); - } - _prf_final_deltas.fill(deltas.size(), deltas.get(deltas.size()-1));
+
//Fill dummy versions SymmetricMatrix covariance = this.getCovariance(rth_cluster, electrodes);
@@ -568,7 +553,7 @@
Hep3Vector _orgloc = this.getOrgLoc(sensor); Hep3Vector h = new BasicHep3Vector(_bfield.x(),_bfield.y(),Math.abs(_bfield.z())); h = VecOp.unit(h);
- Hep3Vector position = wutils.getHelixAndPlaneIntercept(wtrack, _orgloc, eta, h);
+ Hep3Vector position = wtrack.getHelixAndPlaneIntercept(_orgloc, eta, h);
if(debug) { HelicalTrackFit htf = wtrack._htf; List<Double> s = HelixUtils.PathToXPlane(htf, position.x(), 0, 0);
@@ -577,15 +562,7 @@
System.out.println(this.getClass().getSimpleName() + ": Path length to position " + position.toString() + ": s = " + s.get(0)); System.out.println(this.getClass().getSimpleName() + ": Difference between W and helixutils: diffpos " + posdiff.toString() + " ( " + position.toString() + " posOnHelix " + posOnHelix.toString() + " R=" + htf.R()); }
- //Fill some debug stuff - List<Double> deltas = wutils.getDeltas(); - this._delta_itercount.get(sensor).fill(deltas.size()); - IProfile1D dhisto = this._delta_histos.get(sensor); - for(int i=0;i<deltas.size();++i) { - this._prf_all_deltas.fill(i, deltas.get(i)); - dhisto.fill(i,deltas.get(i)); - } - _prf_final_deltas.fill(deltas.size(), deltas.get(deltas.size()-1));
+
return position; }
@@ -701,7 +678,7 @@
this.printDebug("p side position in lcsim coordinates: " + pSidePosition.toString()); //Hep3Vector pSideInter = wutils.getHelixAndPlaneIntercept(wtrack, pSidePosition, VecOp.unit(pSidePosition), h); Hep3Vector eta = this.getPlaneUnitVector(sensor);
- Hep3Vector pSideInter = wutils.getHelixAndPlaneIntercept(wtrack, pSidePosition, eta , h);
+ Hep3Vector pSideInter = wtrack.getHelixAndPlaneIntercept(pSidePosition, eta , h);
this.printDebug("Intersection between track and p side: " + pSideInter.toString()); // Transform back to the JLab coordinates pSideInter = VecOp.mult(VecOp.inverse(detToTrk), pSideInter);
@@ -734,7 +711,7 @@
nSidePosition = VecOp.mult(detToTrk, nSidePosition); this.printDebug("n side position in lcsim coordinates: " + nSidePosition.toString()); //Hep3Vector pSideInter = wutils.getHelixAndPlaneIntercept(wtrack, pSidePosition, VecOp.unit(pSidePosition), h);
- Hep3Vector nSideInter = wutils.getHelixAndPlaneIntercept(wtrack, nSidePosition, eta , h);
+ Hep3Vector nSideInter = wtrack.getHelixAndPlaneIntercept(nSidePosition, eta , h);
this.printDebug("Intersection between track and n side: " + nSideInter.toString()); // Transform back to the JLab coordinates nSideInter = VecOp.mult(VecOp.inverse(detToTrk), nSideInter);
diff -u -r1.11 -r1.12 --- TrackUtils.java 26 Nov 2012 18:12:01 -0000 1.11 +++ TrackUtils.java 27 Nov 2012 02:17:29 -0000 1.12 @@ -19,12 +19,11 @@
import org.lcsim.event.Track; import org.lcsim.fit.helicaltrack.*; import org.lcsim.hps.users.phansson.WTrack;
-import org.lcsim.hps.users.phansson.WTrackUtils;
/** * * @author Omar Moreno <[log in to unmask]>
- * @version $Id: TrackUtils.java,v 1.11 2012/11/26 18:12:01 phansson Exp $
+ * @version $Id: TrackUtils.java,v 1.12 2012/11/27 02:17:29 phansson Exp $
* TODO: Switch to JLab coordinates */
@@ -33,8 +32,7 @@
private boolean _debug = false; private boolean isTrackSet = false; private double[] trackParameters;
- private WTrackUtils _wutils = new WTrackUtils(); -
+
/** *
@@ -48,8 +46,7 @@
*/ public void setDebug(boolean debug){ this._debug = debug;
- _wutils.setDebug(debug); - }
+ }
/** *
@@ -262,79 +259,29 @@
return new BasicHep3Vector(x, y, z); }
- public Hep3Vector calculateIterativeHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield) { - return this.calculateIterativeHelixInterceptXPlane(helfit, strip, bfield,false); - }
- public Hep3Vector calculateIterativeHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield, boolean debug) { - if(debug) System.out.printf("%s: calculateIterativeHelixInterceptXPlane for heltrackfit\n%s \n",this.getClass().getSimpleName(),helfit.toString());
+ public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, Hep3Vector unit_vec_normal_to_plane,Hep3Vector point_on_plane, double bfield) { + /* + * Use code in WTrack to find the iterative solution to the interception + */ + boolean debug = true;
WTrack wtrack = new WTrack(helfit,bfield,true); // B-field sign is flipped so flip!
- if(debug) { - System.out.printf("%s: created Wtrack with params for heltrackfit\n%s \n",this.getClass().getSimpleName(),wtrack.toString()); - System.out.printf("%s: with associated heltrackfit \n%s \n",this.getClass().getSimpleName(),wtrack._htf.toString()); - }
+ if(debug) System.out.printf("getHelixPlaneIntercept:find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, and WTrack \n%s \n",point_on_plane.toString(),unit_vec_normal_to_plane.toString(), bfield,wtrack.toString()); + Hep3Vector intercept_point = wtrack.getHelixAndPlaneIntercept(wtrack,point_on_plane, unit_vec_normal_to_plane, new BasicHep3Vector(0,0,1)); + if(debug) System.out.printf("getHelixPlaneIntercept: found intercept point at %s\n",intercept_point.toString()); + return intercept_point; + } + + public static Hep3Vector getHelixPlaneIntercept(HelicalTrackFit helfit, HelicalTrackStrip strip, double bfield) {
Hep3Vector point_on_plane = strip.origin(); Hep3Vector unit_vec_normal_to_plane = VecOp.cross(strip.u(),strip.v());//strip.w();
- if(debug) System.out.printf("%s: find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, and track pars:\n%s \n",this.getClass().getSimpleName(),point_on_plane.toString(),unit_vec_normal_to_plane.toString(), bfield,helfit.toString()); - boolean d = _wutils.getDebug(); - if(debug) _wutils.setDebug(debug); - Hep3Vector intercept_point = _wutils.getHelixAndPlaneIntercept(wtrack, point_on_plane, unit_vec_normal_to_plane, new BasicHep3Vector(0,0,1)); - _wutils.setDebug(d); //turn off if needed
+ Hep3Vector intercept_point = getHelixPlaneIntercept(helfit, unit_vec_normal_to_plane, point_on_plane, bfield);
return intercept_point; }
+
- public double calculateHelixInterceptXPlane(HelicalTrackFit helfit, HelicalTrackStrip strip) { - /* - * Calculate the interception path length on the helix - * where the track fit and the plane (where the strip is located) intercepts - */ - boolean _DEBUG = false; - - HelicalTrackFit cfit = helfit; - if(cfit==null) { - System.out.println(this.getClass().getSimpleName() + ": no helical track fit was provided for interception calculation"); - System.exit(1); - } - if(_DEBUG) { - System.out.println(this.getClass().getSimpleName() + ": calculateHelixIntercept ---"); - System.out.println(this.getClass().getSimpleName() + ": Track curvatute " + cfit.curvature() + " phi0 " + cfit.phi0()); - } - //Get the origin of the plane - Hep3Vector org = strip.origin(); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": origin of plane " + org.toString()); - - //Get the point on the helix at the origin - double xint = org.x(); - double s = HelixUtils.PathToXPlane(cfit, xint,0,0).get(0); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": xint " + xint + " s " + s ); - - //Calculate that point in space - Hep3Vector pos_xint = HelixUtils.PointOnHelix(cfit, s); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": xint position " + pos_xint.toString()); - - //Create a vector that points from the origin to that point - Hep3Vector vec_xint = VecOp.sub(pos_xint, org); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": vec_xint " + vec_xint.toString()); - - //Get the v in global coordinates - // v tells how the unmeasured coordinate is oriented in the tracking frame - Hep3Vector v = strip.v(); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": v " + v.toString()); - - // We are basically looking for the point where v and track intercepts in the x-y plane - // Assuming that the track travels parallel to the x axis we can calculate the - // the correction to xint by scaling - // NEED TO TAKE INTO ACCOUNT THE REAL TRACK TRAJECTORY!! - double x_corr = vec_xint.y() * v.x()/v.y(); - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": vy " + vec_xint.y() + " vx/vy " + v.x()/v.y() + " -> " + "x_corr " + x_corr); - double xint_new = xint + x_corr; - if(_DEBUG) System.out.println(this.getClass().getSimpleName() + ": final xint_new " + xint_new); - - return xint_new;
- } -
/**
@@ -369,7 +316,7 @@
/** * */
- public double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){
+ public static double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){
return .5*(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1); }
@@ -377,7 +324,7 @@
/** * */
- public boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor){
+ public static boolean sensorContainsTrack(Hep3Vector trackPosition, SiSensor sensor){
boolean debug = false; ITransform3D localToGlobal = sensor.getGeometry().getLocalToGlobal();
@@ -385,8 +332,8 @@
Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid(); Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0); if(debug){
- System.out.println(this.getClass().getSimpleName() + ": Sensor: " + SvtUtils.getInstance().getDescription(sensor)); - System.out.println(this.getClass().getSimpleName() + ": Track Position: " + trackPosition.toString());
+ System.out.println("sensorContainsTrack: Sensor: " + SvtUtils.getInstance().getDescription(sensor)); + System.out.println("sensorContainsTrack: Track Position: " + trackPosition.toString());
} List<Point3D> vertices = new ArrayList<Point3D>();
@@ -399,8 +346,8 @@
//vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z())); vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z())); if(debug){
- System.out.println(this.getClass().getSimpleName() + ": Vertex 1 Position: " + vertices.get(0).toString()); - //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());
+ System.out.println("sensorContainsTrack: Vertex 1 Position: " + vertices.get(0).toString()); + //System.out.println("sensorContainsTrack: Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());
} } else if(vertex.y() > 0 && vertex.x() > 0){
@@ -408,8 +355,8 @@
//vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z())); vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z())); if(debug){
- System.out.println(this.getClass().getSimpleName() + ": Vertex 2 Position: " + vertices.get(1).toString()); - //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());
+ System.out.println("sensorContainsTrack: Vertex 2 Position: " + vertices.get(1).toString()); + //System.out.println("sensorContainsTrack: Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());
} } else if(vertex.y() > 0 && vertex.x() < 0){
@@ -417,8 +364,8 @@
//vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z())); vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z())); if(debug){
- System.out.println(this.getClass().getSimpleName() + ": Vertex 3 Position: " + vertices.get(2).toString()); - //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());
+ System.out.println("sensorContainsTrack: Vertex 3 Position: " + vertices.get(2).toString()); + //System.out.println("sensorContainsTrack: Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());
} } else if(vertex.y() < 0 && vertex.x() < 0){
@@ -426,16 +373,16 @@
//vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z())); vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z())); if(debug){
- System.out.println(this.getClass().getSimpleName() + ": Vertex 4 Position: " + vertices.get(3).toString()); - //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());
+ System.out.println("sensorContainsTrack: Vertex 4 Position: " + vertices.get(3).toString()); + //System.out.println("sensorContainsTrack: Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());
} } }
- double area1 = this.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z()); - double area2 = this.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z()); - double area3 = this.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z()); - double area4 = this.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
+ double area1 = TrackUtils.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z()); + double area2 = TrackUtils.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z()); + double area3 = TrackUtils.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z()); + double area4 = TrackUtils.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
if((area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0)) return true;
@@ -444,16 +391,16 @@
public HelicalTrackFit makeHelicalTrackFit(double[] helix_parameters) {
- SymmetricMatrix cov = new SymmetricMatrix(5); - double[] chisq = new double[2]; - int[] ndf = new int[2]; - Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>(); - Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>(); - HelicalTrackFit fit = new HelicalTrackFit(helix_parameters, cov, chisq, ndf,smap,msmap); - return fit;
+ SymmetricMatrix cov = new SymmetricMatrix(5); + double[] chisq = new double[2]; + int[] ndf = new int[2]; + Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>(); + Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>(); + HelicalTrackFit fit = new HelicalTrackFit(helix_parameters, cov, chisq, ndf,smap,msmap); + return fit;
}
- public Map<String,Double> calculateTrackHitResidual(HelicalTrackHit hth,HelicalTrackFit track, boolean includeMS ) {
+ public static Map<String,Double> calculateTrackHitResidual(HelicalTrackHit hth,HelicalTrackFit track, boolean includeMS ) {
Map<String,Double> residuals = new HashMap<String,Double>(); Map<HelicalTrackHit, MultipleScatter> msmap = track.ScatterMap();
diff -u -r1.1 -r1.2 --- StripMPAlignmentInput.java 19 Nov 2012 21:51:07 -0000 1.1 +++ StripMPAlignmentInput.java 27 Nov 2012 02:17:29 -0000 1.2 @@ -4,31 +4,24 @@
import hep.aida.ref.plotter.PlotterRegion; import hep.physics.matrix.BasicMatrix; import hep.physics.matrix.MatrixOp;
-import hep.physics.vec.*; -import java.io.FileWriter; -import java.io.IOException; -import java.io.PrintWriter;
+import hep.physics.vec.BasicHep3Vector; +import hep.physics.vec.Hep3Matrix; +import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp;
import java.util.ArrayList; import java.util.Arrays; import java.util.List; import java.util.Map;
-import java.util.logging.Level; -import java.util.logging.Logger; -import org.lcsim.constants.Constants;
import org.lcsim.detector.tracker.silicon.SiSensor; import org.lcsim.event.RawTrackerHit; import org.lcsim.event.Track; import org.lcsim.event.TrackerHit; import org.lcsim.fit.helicaltrack.*;
-import org.lcsim.hps.event.HPSTransformations;
import org.lcsim.hps.monitoring.AIDAFrame; import org.lcsim.hps.recon.tracking.SvtUtils; import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.hps.recon.tracking.TrackerHitUtils; -import org.lcsim.hps.users.mgraham.alignment.RunAlignment;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.util.aida.AIDA;
/** * Class to calculate and print the residuals and derivatives
@@ -52,7 +45,6 @@
private double[] _resid = new double[3]; private double[] _error = new double[3];
- private final int _nTrackParameters = 5; //the five track parameters
private AIDAFrame plotterFrame;
@@ -188,7 +180,7 @@
*/ //Find interception with the plane the hit belongs to i.e. the predicted hit position
- Hep3Vector p = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()));
+ Hep3Vector p = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(this._bfield.z()));
double pathLengthToInterception = HelixUtils.PathToXPlane(_trk, p.x(),0,0).get(0); //Find the unit vector of the track direction TrackDirection trkdir = HelixUtils.CalculateTrackDirection(_trk, pathLengthToInterception);
@@ -381,15 +373,13 @@
String side = SvtUtils.getInstance().isTopLayer((SiSensor)((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement()) ? "top" : "bottom"; double bfield = Math.abs(this._bfield.z());
- if(_DEBUG) { - System.out.printf("%s: Finding interception point for residual calculation (B=%s)\n",this.getClass().getSimpleName(),this._bfield.toString()); - double xint_wrong = trackUtil.calculateHelixInterceptXPlane(_trk, strip); - double s_wrong = HelixUtils.PathToXPlane(_trk, xint_wrong, 0, 0).get(0); - Hep3Vector trkpos_wrong= HelixUtils.PointOnHelix(_trk, s_wrong); - System.out.printf("%s: using wrong method xint_wrong=%.3f s_wrong=%.3f -> trkpos_wrong=%s\n",this.getClass().getSimpleName(),xint_wrong,s_wrong,trkpos_wrong.toString()); - }
+// if(_DEBUG) { +// System.out.printf("%s: Finding interception point for residual calculation (B=%s)\n",this.getClass().getSimpleName(),this._bfield.toString()); +// Hep3Vector trkpos_wrong= TrackUtils.helixPositionAtPlaneApprox(_trk, strip); +// System.out.printf("%s: using wrong method xint_wrong=%.3f s_wrong=%.3f -> trkpos_wrong=%s\n",this.getClass().getSimpleName(),trkpos_wrong.toString()); +// }
//Find interception with plane that the strips belongs to
- Hep3Vector trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield);
+ Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bfield);
if(_DEBUG) { System.out.printf("%s: found interception point at %s \n",this.getClass().getSimpleName(),trkpos.toString());
@@ -400,7 +390,7 @@
System.out.printf("%s: failed to get interception point (%s) \n",this.getClass().getSimpleName(),trkpos.toString()); System.out.printf("%s: track params\n%s\n",this.getClass().getSimpleName(),_trk.toString()); System.out.printf("%s: track pT=%.3f chi2=[%.3f][%.3f] \n",this.getClass().getSimpleName(),_trk.pT(bfield),_trk.chisq()[0],_trk.chisq()[1]);
- trkpos = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, bfield,true);
+ trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bfield);
System.exit(1); }
@@ -445,7 +435,7 @@
if (_DEBUG) { System.out.printf("%s: CalculateResidual Result ----\n",this.getClass().getSimpleName());
- System.out.printf("%s: Strip Origin: %s\n",this.getClass().getSimpleName(),corigin.toString());
+ System.out.printf("%s: Strip layer %d Origin: %s\n",this.getClass().getSimpleName(),strip.layer(),corigin.toString());
System.out.printf("%s: Position on the track (tracking coordinates) at the strip: %s\n",this.getClass().getSimpleName(),trkpos.toString()); System.out.printf("%s: vdiffTrk %s\n",this.getClass().getSimpleName(),vdiffTrk.toString()); System.out.printf("%s: vdiff %s\n",this.getClass().getSimpleName(),vdiff.toString());
@@ -805,12 +795,12 @@
private void CalculateGlobalDerivativesWRONG(HelicalTrackStrip strip) { //Find interception with plane that the strips belongs to
- Hep3Vector x_vec = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()));
+ Hep3Vector x_vec = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(this._bfield.z()));
if(Double.isNaN(x_vec.x())) { System.out.printf("%s: error this trkpos is wrong %s\n",this.getClass().getSimpleName(),x_vec.toString()); System.out.printf("%s: origin %s trk \n%s\n",this.getClass().getSimpleName(),strip.origin(),_trk.toString());
- x_vec = trackUtil.calculateIterativeHelixInterceptXPlane(_trk, strip, Math.abs(this._bfield.z()),true);
+ x_vec = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(this._bfield.z()));
System.exit(1); } double slope = _trk.slope();
diff -u -r1.1 -r1.2 --- MPAlignmentInputCalculator.java 19 Nov 2012 21:51:07 -0000 1.1 +++ MPAlignmentInputCalculator.java 27 Nov 2012 02:17:29 -0000 1.2 @@ -41,7 +41,6 @@
protected AlignmentUtils _alignUtils; protected AlignmentUtils.OldAlignmentUtils _oldAlignUtils; protected AlignmentUtils.NumDerivatives _numDerivatives;
- protected TrackUtils trackUtil;
protected TrackerHitUtils trackerHitUtil; protected boolean hideFrame = false; protected AIDA aida = AIDA.defaultInstance();
@@ -61,7 +60,6 @@
_alignUtils = new AlignmentUtils(_DEBUG); _oldAlignUtils = new AlignmentUtils(_DEBUG).new OldAlignmentUtils(); _numDerivatives = new AlignmentUtils(_DEBUG).new NumDerivatives();
- trackUtil = new TrackUtils();
trackerHitUtil = new TrackerHitUtils(_DEBUG); _bfield = new BasicHep3Vector(0., 0., 1.); _type = type;
diff -u -r1.4 -r1.5 --- WTrack.java 7 Nov 2012 20:54:13 -0000 1.4 +++ WTrack.java 27 Nov 2012 02:17:29 -0000 1.5 @@ -7,6 +7,8 @@
import hep.physics.vec.BasicHep3Vector; import hep.physics.vec.Hep3Vector; import hep.physics.vec.VecOp;
+import java.util.ArrayList; +import java.util.List;
import org.lcsim.constants.Constants; import org.lcsim.fit.helicaltrack.HelicalTrackFit;
@@ -21,11 +23,13 @@
public enum PARAM{TEST;} private double[] _parameters = new double[7]; public HelicalTrackFit _htf = null;
- double _bfield; - double _bfield_constant; -
+ private double _bfield; + private double _bfield_constant;
int _q = 0;
+ static int max_iterations_intercept = 10; + static double epsilon_intercept = 1e-4; +
public WTrack(double [] params, double bfield, int q) { _bfield = bfield; _q = q;
@@ -36,8 +40,10 @@
public WTrack(WTrack trk) { _bfield = trk._bfield; _q = trk.getCharge();
- _bfield_constant = -1*Constants.fieldConversion*_bfield*_q;
+ _bfield_constant = trk._bfield_constant;
_parameters = trk.getParameters();
+ _htf = trk._htf; + _debug = trk._debug;
} public WTrack(HelicalTrackFit track, double bfield) {
@@ -113,13 +119,16 @@
return ( new BasicHep3Vector(_parameters[4],_parameters[5],_parameters[6])); }
+ public String paramsToString() { + String str = ""; + for(int i=0;i<7;++i) str += _parameters[i] + ", "; + return str; + }
public String toString() {
- String str = "WTrack params ["; - for(int i=0;i<7;++i) str += _parameters[i] + ", "; - str += "]";
+ String str = "WTrack params [" + paramsToString() + "]";
if(this._htf!=null) {
- str += "\n WTrack corresponding HelicalTrackFit params:\n";
+ str += "\n with corresponding HelicalTrackFit:\n";
str += this._htf.toString(); } return str;
@@ -127,6 +136,212 @@
+ + + + 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)); + } + + private 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 + */ + boolean debug = true; + double a = track.getBfieldConstant(); + Hep3Vector p0 = track.getP0(); + Hep3Vector x0 = track.getX0(); + double p = p0.magnitude(); + double rho = a / p; + 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 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(" xp " + xp.toString()); + System.out.println(" eta " + eta.toString()); + System.out.println(" h " + h.toString()); + System.exit(1); + } + double root1 = (-B + Math.sqrt(t)) /(2*A); + double root2 = (-B - Math.sqrt(t)) /(2*A); + + // choose the smallest positive solution + // if both negative choose the smallest negative ??? + //if(root1==0 || root2==0) root=0; + double 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(" I should never get here! (root1=" + root1 + " root2=" + root2+")"); +// System.exit(1); +// } + if((debug)) { + System.out.println(" getPathLengthToPlaneApprox "); + System.out.println(" " + track.toString()); + System.out.println(" xp " + xp.toString()); + System.out.println(" eta " + eta.toString()); + System.out.println(" h " + h.toString()); + System.out.println(" p " + p + " rho " + rho + " t " + t + " A " + A + " B " + B + " C " + C); + System.out.println(" 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; + } + + + private 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 + */ + boolean debug = true; + // Find track parameters at that path length + Hep3Vector p = getMomentumOnHelix(track, s, h); + Hep3Vector x = getPointOnHelix(track,s, h); + + Hep3Vector p_tmp = getMomentumOnHelix(track, s); + Hep3Vector x_tmp = getPointOnHelix(track, s); + + if((debug)) { + System.out.println(" point on helix at s"); + System.out.println(" p " + p.toString() + " p_tmp " + p_tmp.toString()); + System.out.println(" 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(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 + * xp: point on the plane + * eta: unit vector of the plane + * h: unit vector of magnetic field + */ + boolean debug = true; + + + 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_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); + + // 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: updated trak params: [%s]\n",this.getClass().getSimpleName(),iteration,trk.paramsToString()); + + //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()); + + + } + + 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()); + + Hep3Vector point = trk.getX0(); + return point; + + } +
diff -u -r1.1 -r1.2 --- ModuleMPAlignmentInput.java 19 Nov 2012 21:51:07 -0000 1.1 +++ ModuleMPAlignmentInput.java 27 Nov 2012 02:17:29 -0000 1.2 @@ -18,6 +18,7 @@
import org.lcsim.fit.helicaltrack.*; import org.lcsim.hps.monitoring.AIDAFrame; import org.lcsim.hps.recon.tracking.SvtUtils;
+import org.lcsim.hps.recon.tracking.TrackUtils;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate; import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -335,7 +336,7 @@
private void CalculateResidual(HelicalTrackHit hit) { if(_DEBUG) System.out.printf("%s: --- Calculate Residual for HelicalTrackhit ---\n",this.getClass().getSimpleName());
- Map<String,Double> res_map = this.trackUtil.calculateTrackHitResidual(hit, _trk, _includeMS);
+ Map<String,Double> res_map = TrackUtils.calculateTrackHitResidual(hit, _trk, _includeMS);
_resid[0] = 0.; _error[0] = 0.320/Math.sqrt(12); _resid[1] = res_map.get("resy");
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