lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.31 -r1.32
--- HelicalTrackFitter.java 19 Jan 2011 20:44:05 -0000 1.31
+++ HelicalTrackFitter.java 1 Feb 2011 22:54:24 -0000 1.32
@@ -4,7 +4,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.31 2011/01/19 20:44:05 mgraham Exp $
+ * $Id: HelicalTrackFitter.java,v 1.32 2011/02/01 22:54:24 partridge Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
@@ -43,8 +43,8 @@
* @author Norman Graf
* @version 2.0 (R. Partridge)
*/
-
public class HelicalTrackFitter {
+
private CircleFitter _cfitter = new CircleFitter();
private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
private ZSegmentFitter _zfitter = new ZSegmentFitter();
@@ -58,6 +58,7 @@
* Status of the HelicalTrackFit.
*/
public enum FitStatus {
+
/**
* Successful Fit.
*/
@@ -77,7 +78,8 @@
/**
* ZSegmentFit failed.
*/
- ZSegmentFitFailed};
+ ZSegmentFitFailed
+ };
/**
* Creates a new instance of HelicalTrackFitter.
@@ -101,7 +103,7 @@
*/
public FitStatus fit(double[] x, double[] y, double[] z, double[] drphi, double[] dz, int np) {
List<HelicalTrackHit> hitcol = new ArrayList<HelicalTrackHit>();
- for(int i=0; i<np; i++) {
+ for (int i = 0; i < np; i++) {
Hep3Vector pos = new BasicHep3Vector(x[i], y[i], z[i]);
if (dz[i] > 0.) {
HelicalTrackHit hit = new HelicalTrack3DHit(pos, MakeCov(pos, drphi[i], 0., dz[i]),
@@ -143,6 +145,23 @@
*/
public FitStatus fit(List<HelicalTrackHit> hitcol, Map<HelicalTrackHit, MultipleScatter> msmap, HelicalTrackFit oldhelix) {
+ // Check that we have at least 3 hits
+ int nhit = hitcol.size();
+
+ // 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
+ // It is also assumed that the track won't curl through an angle > 180 degrees between
+ // neighboring points. This might be a problem for curlers with small dip angle.
+ Collections.sort(hitcol);
+
+ // See if the first hit is closer to the origin than the last hit
+ // If not, reverse the order of the hits
+ double zfirst = hitcol.get(0).z();
+ double zlast = hitcol.get(nhit - 1).z();
+ if (Math.abs(zfirst) > Math.abs(zlast)) {
+ Collections.reverse(hitcol);
+ }
+
// Initialize the various fitter outputs
_cfit = null;
_lfit = null;
@@ -156,30 +175,21 @@
// Sort the hits into the appropriate lists
for (HelicalTrackHit hit : hitcol) {
- // Hits to be used in the circle fit
- if (hit.drphi() > 0) circle_hits.add(hit);
+ // Hits to be used in the circle fit
+ if (hit.drphi() > 0) circle_hits.add(hit);
// Pixel hits
- if (hit instanceof HelicalTrack3DHit) pixel_hits.add(hit);
+ if (hit instanceof HelicalTrack3DHit) pixel_hits.add(hit);
// Strip hits
- if (hit instanceof HelicalTrack2DHit) strip_hits.add(hit);
+ if (hit instanceof HelicalTrack2DHit) strip_hits.add(hit);
// Cross hits
- if (hit instanceof HelicalTrackCross) pixel_hits.add(hit);
- }
+ if (hit instanceof HelicalTrackCross) pixel_hits.add(hit);
+ }
// Check to make sure we have at least 3 circle hits before proceeding
int nc = circle_hits.size();
- if (nc < 3) return FitStatus.CircleFitFailed;
-
- // 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
- // It is also assumed that the track won't curl through an angle > 180 degrees between
- // neighboring points. This might be a problem for curlers with small dip angle.
- Collections.sort(circle_hits);
- // See if the first hit is closer to the origin than the last hit
- // If not, reverse the order of the hits
- double zfirst = circle_hits.get(0).z();
- double zlast = circle_hits.get(nc-1).z();
- if (Math.abs(zfirst) > Math.abs(zlast)) Collections.reverse(circle_hits);
+ if (nc < 3) {
+ return FitStatus.CircleFitFailed;
+ }
// Create the objects that will hold the fit output
double[] chisq = new double[2];
@@ -193,7 +203,7 @@
double[] wrphi = new double[nc];
// Store the hit coordinates and weights in arrays for the circle fitter
- for(int i=0; i<nc; i++) {
+ for (int i = 0; i < nc; i++) {
HelicalTrackHit hit = circle_hits.get(i);
// Store the hit position
x[i] = hit.x();
@@ -201,18 +211,43 @@
// Find the weight (= 1/uncertainty^2) for this hit
// First get the multiple scattering uncertainty
double drphi_ms = 0.;
- if (msmap.containsKey(hit)) drphi_ms = msmap.get(hit).drphi();
+ if (msmap.containsKey(hit)) {
+ drphi_ms = msmap.get(hit).drphi();
+ }
// Get the hit resolution and combine uncertainties in quadrature
double drphi_res = hit.drphi();
- wrphi[i] = 1. / (drphi_res*drphi_res + drphi_ms*drphi_ms);
+ wrphi[i] = 1. / (drphi_res * drphi_res + drphi_ms * drphi_ms);
}
// Call the circle fitter and check for success
boolean success = _cfitter.fit(x, y, wrphi, nc);
- if(!success) return FitStatus.CircleFitFailed;
+ if (!success) {
+ return FitStatus.CircleFitFailed;
+ }
+
+ // Get the results of the fit
+ _cfit = _cfitter.getfit();
+
+ // Calculate the arc lengths from the DCA to each hit and check for backwards hits
+ Map<HelicalTrackHit, Double> smap = getPathLengths(hitcol);
// If we are going around the circle in the wrong direction, fix it
- _cfit = CircleFix(_cfitter.getfit(), circle_hits);
+ if (smap.get(circle_hits.get(nc-1)) < 0.) {
+
+ // Change the circle fit parameters to reverse direction
+ _cfit = CircleFix(_cfitter.getfit());
+
+ // Flip the signs of the path lengths
+ for (HelicalTrackHit hit : hitcol) {
+ double oldpath = smap.get(hit);
+ smap.put(hit, -oldpath);
+ }
+ }
+
+ // Check that things make sense
+ for (HelicalTrackHit hit : smap.keySet()) {
+ if (smap.get(hit) < 0.) return FitStatus.InconsistentSeed;
+ }
// Save the chi^2 and dof for the circle fit
chisq[0] = _cfit.chisq();
@@ -238,12 +273,6 @@
cov.setElement(HelicalTrackFit.phi0Index, HelicalTrackFit.dcaIndex, -1. * _cfit.cov()[4]); // fix d0 sign convention
cov.setElement(HelicalTrackFit.dcaIndex, HelicalTrackFit.dcaIndex, _cfit.cov()[5]);
- // Calculate the arc lengths from the DCA to each hit and check for backwards hits
- Map<HelicalTrackHit, Double> smap = getPathLengths(hitcol);
- for (HelicalTrackHit hit : smap.keySet()) {
- if (smap.get(hit) < 0.) return FitStatus.InconsistentSeed;
- }
-
// Check if we have enough pixel hits to do a straight-line fit of s vs z
int npix = pixel_hits.size();
if (npix > 1) {
@@ -254,7 +283,7 @@
double[] dz = new double[npix];
// Store the coordinates and errors for the line fit
- for(int i=0; i<npix; i++) {
+ for (int i = 0; i < npix; i++) {
HelicalTrackHit hit = pixel_hits.get(i);
z[i] = hit.z();
dz[i] = HitUtils.zres(hit, msmap, oldhelix);
@@ -263,7 +292,9 @@
// Call the line fitter and check for success
success = _lfitter.fit(s, z, dz, npix);
- if(!success) return FitStatus.LineFitFailed;
+ if (!success) {
+ return FitStatus.LineFitFailed;
+ }
// Save the line fit, chi^2, and DOF
_lfit = _lfitter.getFit();
@@ -275,9 +306,9 @@
par[HelicalTrackFit.slopeIndex] = _lfit.slope();
// Save the line fit covariance matrix elements
- cov.setElement(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index, Math.pow(_lfit.interceptUncertainty(),2));
+ cov.setElement(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index, Math.pow(_lfit.interceptUncertainty(), 2));
cov.setElement(HelicalTrackFit.z0Index, HelicalTrackFit.slopeIndex, _lfit.covariance());
- cov.setElement(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex, Math.pow(_lfit.slopeUncertainty(),2));
+ cov.setElement(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex, Math.pow(_lfit.slopeUncertainty(), 2));
} else {
@@ -293,7 +324,9 @@
// We should always have enough hits for a ZSegment fit
int nstrip = strip_hits.size();
- if (nstrip < 2) throw new RuntimeException("Too few hits for a ZSegment fit");
+ if (nstrip < 2) {
+ throw new RuntimeException("Too few hits for a ZSegment fit");
+ }
// Setup for the ZSegment fit
double[] s = new double[nstrip];
@@ -301,19 +334,23 @@
double[] zmax = new double[nstrip];
// Store the path lengths and z limits for the ZSegment fit
- for(int i=0; i<nstrip; i++) {
+ for (int i = 0; i < nstrip; i++) {
HelicalTrack2DHit hit = (HelicalTrack2DHit) strip_hits.get(i);
- s[i] = smap.get(hit);
- // Get the multiple scattering uncertainty and adjust z limits accordingly
+ s[i] = smap.get(hit);
+ // Get the multiple scattering uncertainty and adjust z limits accordingly
double dz = 0.;
- if (msmap.containsKey(hit)) dz = msmap.get(hit).dz();
- zmin[i] = hit.zmin() - _tolerance * dz;
- zmax[i] = hit.zmax() + _tolerance * dz;
+ if (msmap.containsKey(hit)) {
+ dz = msmap.get(hit).dz();
}
+ zmin[i] = hit.zmin() - _tolerance * dz;
+ zmax[i] = hit.zmax() + _tolerance * dz;
+ }
// Call the ZSegment fitter and check for success
success = _zfitter.fit(s, zmin, zmax);
- if(!success) return FitStatus.ZSegmentFitFailed;
+ if (!success) {
+ return FitStatus.ZSegmentFitFailed;
+ }
// Save the ZSegment fit, chi^2, and DOF
_zfit = _zfitter.getFit();
@@ -328,13 +365,15 @@
cov.setElement(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index,
_zfit.getCovariance().e(0, 0));
cov.setElement(HelicalTrackFit.z0Index, HelicalTrackFit.slopeIndex,
- _zfit.getCovariance().e(0,1));
+ _zfit.getCovariance().e(0, 1));
cov.setElement(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex,
_zfit.getCovariance().e(1, 1));
}
// Include a chisq penalty in the s-z fit if we miss any axial strips
- if (strip_hits.size() > 0) chisq[1] += MissedStripPenalty(strip_hits, par, cov, smap, msmap);
+ if (strip_hits.size() > 0) {
+ chisq[1] += MissedStripPenalty(strip_hits, par, cov, smap, msmap);
+ }
// Create the HelicalTrackFit for this helix
_fit = new HelicalTrackFit(par, cov, chisq, ndof, smap, msmap);
@@ -421,10 +460,10 @@
// Create a new covariance matrix and set the non-zero elements
SymmetricMatrix cov = new SymmetricMatrix(3);
- cov.setElement(0, 0, (y*y * drphi*drphi + x*x * dr*dr) / r2);
- cov.setElement(0, 1, x * y * (dr*dr - drphi*drphi) / r2);
- cov.setElement(1, 1, (x*x * drphi*drphi + y*y * dr*dr) / r2);
- cov.setElement(2, 2, dz*dz);
+ cov.setElement(0, 0, (y * y * drphi * drphi + x * x * dr * dr) / r2);
+ cov.setElement(0, 1, x * y * (dr * dr - drphi * drphi) / r2);
+ cov.setElement(1, 1, (x * x * drphi * drphi + y * y * dr * dr) / r2);
+ cov.setElement(2, 2, dz * dz);
return cov;
}
@@ -440,18 +479,14 @@
* @param hitlist list of hits used for the circle fit
* @return fixed circle fit
*/
- private CircleFit CircleFix(CircleFit oldfit, List<HelicalTrackHit> hitlist) {
-
- // Check if we are going around the circle in the right direction by getting
- // the path length to the last hit and checking that it is positive
- int nhits = hitlist.size();
- double s = HelixUtils.PathLength(oldfit, hitlist.get(nhits - 1));
- if (s > 0.) return oldfit;
+ private CircleFit CircleFix(CircleFit oldfit) {
// Reverse the direction by changing the sign of dca, curv, and adding pi to phi0
double dca = -oldfit.dca();
double phi0 = oldfit.phi() + Math.PI;
- if (phi0 > 2. * Math.PI) phi0 -= 2. * Math.PI;
+ if (phi0 > 2. * Math.PI) {
+ phi0 -= 2. * Math.PI;
+ }
double curv = -oldfit.curvature();
// Also fix the affected covariance matrix elements
@@ -473,18 +508,37 @@
// Create a map to store the arc lengths
Map<HelicalTrackHit, Double> smap = new HashMap<HelicalTrackHit, Double>();
- // Loop over the hits and store the arc lengths
- double s = 0.;
- for(int i=0; i<hits.size(); i++) {
- // For the first hit, measure from the DCA
- if (i == 0) {
- s = HelixUtils.PathLength(_cfit, hits.get(i));
+ // Initialize looper tracking and iterate over ordered list of hits
+ double slast = 0.;
+ int ilast = -1;
+ double s;
+ for (int i = 0; i < hits.size(); i++) {
+
+ // Retrieve the next hit ordered by z coordinate and check hit type
+ HelicalTrackHit hit = hits.get(i);
+ if (hit instanceof HelicalTrack2DHit) {
+
+ // Axial hit - measure from the DCA (can't handle loopers)
+ s = HelixUtils.PathLength(_cfit, hit);
+
} else {
- // For subsequent hits, add in the arc length from the previous hit
- s += HelixUtils.PathLength(_cfit, hits.get(i-1), hits.get(i));
+
+ if (ilast < 0) {
+ // For the first 3D hit, measure from the DCA
+ s = HelixUtils.PathLength(_cfit, hit);
+
+ } else {
+ // For subsequent hits, add in the arc length from the previous 3D hit
+ s = slast + HelixUtils.PathLength(_cfit, hits.get(ilast), hit);
+ }
+
+ // Update info on the last 3D hit
+ ilast = i;
+ slast = s;
}
+
// Save the arc length for this hit
- smap.put(hits.get(i), s);
+ smap.put(hit, s);
}
return smap;
}
@@ -515,20 +569,26 @@
// Find the predicted z coordinate
double zhelix = z0 + s * slope;
// Find the uncertainty^2 in the z coordinate
- double dzsq = cov_z0z0 + 2. * s * cov_z0sl + s*s * cov_slsl;
+ double dzsq = cov_z0z0 + 2. * s * cov_z0sl + s * s * cov_slsl;
// Find the multiple scattering error
double dz_ms = 0.;
- if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
+ if (msmap.containsKey(hit)) {
+ dz_ms = msmap.get(hit).dz();
+ }
// Add the multiple scattering uncertainty
- dzsq += dz_ms*dz_ms;
+ dzsq += dz_ms * dz_ms;
// Find the limits in z for the strip
double zmin = ((HelicalTrack2DHit) hit).zmin();
double zmax = ((HelicalTrack2DHit) hit).zmax();
// Calculate a chisq penalty if the predicted z is not on the strip
- if (zhelix < zmin) chisq += (zhelix-zmin)*(zhelix-zmin) / dzsq;
- if (zhelix > zmax) chisq += (zhelix-zmax)*(zhelix-zmax) / dzsq;
+ if (zhelix < zmin) {
+ chisq += (zhelix - zmin) * (zhelix - zmin) / dzsq;
}
+ if (zhelix > zmax) {
+ chisq += (zhelix - zmax) * (zhelix - zmax) / dzsq;
}
+ }
+ }
return chisq;
}
-}
\ No newline at end of file
+}