lcsim/src/org/lcsim/fit/helicaltrack
diff -N HelixUtils.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixUtils.java 4 Feb 2008 20:50:01 -0000 1.1
@@ -0,0 +1,183 @@
+/*
+ * HelixUtils.java
+ *
+ * Created on February 1, 2008, 5:48 PM
+ *
+ */
+
+package org.lcsim.fit.helicaltrack;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.lcsim.fit.circle.CircleFit;
+
+/**
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class HelixUtils {
+ private CircleFit _cfit;
+ private HelicalTrackFit _helix;
+ private double _dca;
+ private double _phi0;
+ private double _omega;
+ private double _RC;
+ private double _z0;
+ private double _slope;
+ private double _xr;
+ private double _yr;
+ private double _xc;
+ private double _yc;
+ private double _x0;
+ private double _y0;
+ private double _sth;
+ private double _cth;
+
+ /** Creates a new instance of HelixUtils */
+ public HelixUtils() {
+ }
+
+ public double PathLength(CircleFit cfit, HelicalTrackHit hit1, HelicalTrackHit hit2) {
+
+ // Cache the circle parameters
+ CircleCache(cfit);
+
+ // Return the path length between hit1 and hit2
+ return PathCalc(hit1.x(), hit1.y(), hit2.x(), hit2.y());
+ }
+
+ public double PathLength(CircleFit cfit, HelicalTrackHit hit) {
+ // Cache the circle parameters
+ CircleCache(cfit);
+
+ // Return the path length from the DCA
+ return PathCalc(_x0, _y0, hit.x(), hit.y());
+ }
+
+ public double PathToZPlane(HelicalTrackFit helix, double z) {
+
+ // Cache the helix parameters
+ HelixCache(helix);
+
+ // Find the z distace between the DCA and the z plane and calculate the x-y path length
+ double zdist = z - _z0;
+ double safeslope = _slope;
+ if (Math.abs(_slope) < 1e-6) safeslope = 1e-6 * Math.signum(_slope);
+ return zdist / safeslope;
+ }
+
+ public List<Double> PathToCylinder(HelicalTrackFit helix, double r, double smax) {
+
+ // Cache the helix parameters
+ HelixCache(helix);
+
+ // Create a list to hold the path lengths
+ List<Double> pathlist = new ArrayList<Double>();
+
+ // Find the angle from the point of closest approach to the first crossing
+ double cdphi = 1 - (r * r - _dca * _dca) / (2. * _RC * (_RC - _dca));
+
+ if (Math.abs(cdphi) < 1.) {
+ double dphi = Math.acos(cdphi);
+ Double s = dphi * Math.abs(_RC);
+ if (s < smax) pathlist.add(s);
+ s += 2. * (Math.PI - dphi);
+ if (s < smax) pathlist.add(s);
+ s += 2. * dphi;
+ while (s < smax && pathlist.size() < 10) {
+ pathlist.add(s);
+ s += 2. * (Math.PI - dphi);
+ if (s < smax) pathlist.add(s);
+ s += 2. * dphi;
+ }
+ }
+ return pathlist;
+ }
+
+ public Hep3Vector Direction(HelicalTrackFit helix, double s) {
+ double phi = _phi0 - s / _RC;
+ double ux = Math.cos(phi) * _sth;
+ double uy = Math.sin(phi) * _sth;
+ double uz = _cth;
+ return new BasicHep3Vector(ux, uy, uz);
+ }
+
+ public Hep3Vector PointOnHelix(HelicalTrackFit helix, double s) {
+
+ // Cache the helix parameters
+ HelixCache(helix);
+
+ // Find the azimuthal direction at this path length
+ double phi = _phi0 - s / _RC;
+
+ // Calculate the position on the helix at this path length
+ double x = _xc - _RC * Math.sin(phi);
+ double y = _yc + _RC * Math.cos(phi);
+ double z = _z0 + s * _slope;
+
+ // Return the position as a Hep3Vector
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ private double PathCalc(double x1, double y1, double x2, double y2) {
+
+ // Find the angle between these points measured wrt the circle center
+ double phi1 = Math.atan2(y1 - _yc, x1 - _xc);
+ double phi2 = Math.atan2(y2 - _yc, x2 - _xc);
+ double dphi = phi2 - phi1;
+
+ // Make sure dphi is in the valid range (-pi, pi)
+ if (dphi > Math.PI) dphi -= 2. * Math.PI;
+ if (dphi < -Math.PI) dphi += 2. * Math.PI;
+
+ // Return the arc length
+ double s = - _RC * dphi;
+ return s;
+ }
+
+ private void CircleCache(CircleFit cfit) {
+
+ // Cache some calculations for this circle if we haven't already saved them
+ if (_cfit != cfit) {
+ _cfit = cfit;
+ _helix = null;
+ _dca = _cfit.dca();
+ _phi0 = _cfit.phi();
+ _RC = 1.0 / _cfit.curvature();
+ _xr = _cfit.xref();
+ _yr = _cfit.yref();
+ _xc = _xr + (_RC - _dca) * Math.sin(_phi0);
+ _yc = _yr - (_RC - _dca) * Math.cos(_phi0);
+ _x0 = _xr - _dca * Math.sin(_phi0);
+ _y0 = _yr + _dca * Math.cos(_phi0);
+ }
+ return;
+ }
+
+ private void HelixCache(HelicalTrackFit helix) {
+ if (_helix != helix) {
+ _cfit = null;
+ _helix = helix;
+ double[] param = _helix.parameters();
+ _dca = _helix.dca();
+ _phi0 = _helix.phi0();
+ _RC = 1.0 / _helix.curvature();
+ _z0 = _helix.z0();
+ _slope = _helix.slope();
+ _xr = 0.;
+ _yr = 0.;
+ _xc = _xr + (_RC - _dca) * Math.sin(_phi0);
+ _yc = _yr - (_RC - _dca) * Math.cos(_phi0);
+ _x0 = _xr - _dca * Math.sin(_phi0);
+ _y0 = _yr + _dca * Math.cos(_phi0);
+ _sth = _helix.sth();
+ _cth = _helix.cth();
+ }
+ return;
+ }
+}
\ No newline at end of file
lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.5 -r1.6
--- HelicalTrackFit.java 11 Dec 2007 20:35:04 -0000 1.5
+++ HelicalTrackFit.java 4 Feb 2008 20:50:01 -0000 1.6
@@ -3,7 +3,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFit.java,v 1.5 2007/12/11 20:35:04 partridge Exp $
+ * $Id: HelicalTrackFit.java,v 1.6 2008/02/04 20:50:01 partridge Exp $
*/
package org.lcsim.fit.helicaltrack;
@@ -47,6 +47,27 @@
public double[] parameters() {
return _parameters;
}
+
+ public double dca() {
+ return _parameters[0];
+ }
+
+ public double phi0() {
+ return _parameters[1];
+ }
+
+ public double curvature() {
+ return _parameters[2];
+ }
+
+ public double z0() {
+ return _parameters[3];
+ }
+
+ public double slope() {
+ return _parameters[4];
+ }
+
public SymmetricMatrix covariance() {
return _covmatrix;
}
@@ -63,6 +84,22 @@
return _smap;
}
+ public double cth() {
+ return _parameters[4] / Math.sqrt(1 + Math.pow(_parameters[4], 2));
+ }
+
+ public double sth() {
+ return 1. / Math.sqrt(1 + Math.pow(_parameters[4], 2));
+ }
+
+ public double pT(double bfield) {
+ return 0.3 * bfield / Math.abs(_parameters[2]);
+ }
+
+ public double p(double bfield) {
+ return pT(bfield) / sth();
+ }
+
public String toString() {
StringBuffer sb = new StringBuffer("HelicalTrackFit: \n");
@@ -74,7 +111,4 @@
return sb.toString();
}
-
-
-
-}
+}
\ No newline at end of file
lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.13 -r1.14
--- HelicalTrackFitter.java 11 Dec 2007 20:51:25 -0000 1.13
+++ HelicalTrackFitter.java 4 Feb 2008 20:50:01 -0000 1.14
@@ -4,7 +4,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.13 2007/12/11 20:51:25 partridge Exp $
+ * $Id: HelicalTrackFitter.java,v 1.14 2008/02/04 20:50:01 partridge Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
@@ -45,6 +45,7 @@
private CircleFitter _cfitter = new CircleFitter();
private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
private ZSegmentFitter _zfitter = new ZSegmentFitter();
+ private HelixUtils _hutil = new HelixUtils();
private HelicalTrackFit _fit;
@@ -79,12 +80,12 @@
Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
return fit(hitcol, msmap);
}
-
- public boolean fit(List<HelicalTrackHit> hitcol) {
- Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
- return fit(hitcol, msmap);
- }
-
+
+ public boolean fit(List<HelicalTrackHit> hitcol) {
+ Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
+ return fit(hitcol, msmap);
+ }
+
public boolean fit(List<HelicalTrackHit> hitcol, Map<HelicalTrackHit, MultipleScatter> msmap) {
// Create lists for the various types of hits
@@ -120,9 +121,8 @@
x[i] = hit.x();
y[i] = hit.y();
double drphi_ms = 0.;
- if (msmap.containsKey(hit)) drphi_ms = msmap.get(hit).bendMS();
+ if (msmap.containsKey(hit)) drphi_ms = msmap.get(hit).bendMS();
wrphi[i] = 1. / (Math.pow(hit.drphi(),2) + Math.pow(drphi_ms,2));
- System.out.println(" Circle fit - x: "+x[i]+" y "+y[i]+" w "+wrphi[i]);
}
// Call the circle fitter and check for success
@@ -142,9 +142,8 @@
cov.setElement(0, 2, cfit.cov()[3]);
cov.setElement(1, 2, cfit.cov()[4]);
cov.setElement(2, 2, cfit.cov()[5]);
- System.out.println(cfit.toString());
// Calculate the arc lengths from the DCA to each hit
- Map<HelicalTrackHit, Double> smap = getPathLengths(par, hitcol);
+ Map<HelicalTrackHit, Double> smap = getPathLengths(cfit, hitcol);
// Check if we have enough pixel hits to do a straight-line fit of s(z)
int npix = pixel_hits.size();
@@ -163,13 +162,12 @@
if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).straightMS();
dz[i] = Math.sqrt(Math.pow(hit.dz(),2)+Math.pow(dz_ms,2));
s[i] = smap.get(hit);
- System.out.println("Line fit - z: "+z[i]+" s "+s[i]+" dz "+dz[i]);
}
// Call the line fitter and check for success
success = _lfitter.fit(s, z, dz, npix);
if(!success) return false;
-
+
// TODO: see if the strip hits are consistent with the track hits
// Store the line fit output
@@ -200,12 +198,12 @@
// Setup for the ZSegment fit
int nstrip = strip_hits.size();
if (nstrip > 1) {
- double[] s = new double[npix];
- double[] zmin = new double[npix];
- double[] zmax = new double[npix];
+ double[] s = new double[nstrip];
+ double[] zmin = new double[nstrip];
+ double[] zmax = new double[nstrip];
// Store the coordinates and errors for the ZSegment fit
- for(int i=0; i<npix; i++) {
+ for(int i=0; i<nstrip; i++) {
HelicalTrack2DHit hit = (HelicalTrack2DHit) strip_hits.get(i);
s[i] = smap.get(hit);
zmin[i] = hit.zmin();
@@ -244,7 +242,7 @@
/**
* Find the arc lengths for a list of hits
*/
- public Map<HelicalTrackHit, Double> getPathLengths(double[] par, List<HelicalTrackHit> hits) {
+ public Map<HelicalTrackHit, Double> getPathLengths(CircleFit cfit, List<HelicalTrackHit> hits) {
// Sort the hits to be monotonic in z so that we can follow a curling track
// It is assumed that the first hit on a track is closer to the origin than the last hit
@@ -255,38 +253,25 @@
double zfirst = hits.get(0).z();
double zlast = hits.get(nhits-1).z();
if (Math.abs(zfirst) > Math.abs(zlast)) Collections.reverse(hits);
-
+
// Create a map to store the arc lengths
Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>();
- // Find the circle radius and the cartesian coordinates of the circle origin and
- // point of closest approach
- double dca = par[0];
- double phi0 = par[1];
- double radius = 1. / par[2];
- double xc = (radius - dca) * Math.sin(phi0);
- double yc = -(radius - dca) * Math.cos(phi0);
- double x0 = dca * Math.sin(phi0);
- double y0 = dca * Math.cos(phi0);
-
- // Measure the arc length starting from the point of closest approach
- double last_phi = Math.atan2(y0 - yc, x0 - xc);
- double last_s = 0.;
+ // Loop over the hits and store the arc lengths
+ double s = 0.;
for(int i=0; i<nhits; i++) {
+ if (i == 0) {
+
+ // For the first hit, measure from the DCA
+ s += _hutil.PathLength(cfit, hits.get(i));
+ } else {
+
+ // For subsequent hits, add in the arc length from the previous hit
+ s+= _hutil.PathLength(cfit, hits.get(i-1), hits.get(i));
+ }
- // Find the arc lengths between a hit and the previous hit
- // Note: positive charged tracks rotate clockwise (i.e., dphi < 0)
- HelicalTrackHit hit = hits.get(i);
- double phi = Math.atan2(hit.y() - yc, hit.x() - xc);
- double dphi = phi - last_phi;
- if (dphi > Math.PI) dphi -= 2. * Math.PI;
- if (dphi < -Math.PI) dphi += 2. * Math.PI;
- double s = last_s - radius * dphi;
-
- // Save the arc length for this hit and update the previous hit info
- smap.put(hit, s);
- last_phi = phi;
- last_s = s;
+ // Save the arc length for this hit
+ smap.put(hits.get(i), s);
}
return smap;
}