lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.5 -r1.6
--- HelicalTrackFitter.java 23 Oct 2006 19:42:36 -0000 1.5
+++ HelicalTrackFitter.java 1 Aug 2007 20:13:53 -0000 1.6
@@ -1,12 +1,12 @@
+package org.lcsim.fit.helicaltrack;
/*
* HelicalTrackFitter.java
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.5 2006/10/23 19:42:36 tonyj Exp $
+ * $Id: HelicalTrackFitter.java,v 1.6 2007/08/01 20:13:53 lstevens Exp $
*/
-package org.lcsim.fit.helicaltrack;
import hep.physics.matrix.SymmetricMatrix;
import static java.lang.Math.atan2;
import static java.lang.Math.abs;
@@ -15,13 +15,30 @@
import static java.lang.Math.sin;
import org.lcsim.fit.circle.CircleFit;
import org.lcsim.fit.circle.CircleFitter;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.fit.line.SlopeInterceptLineFit;
import org.lcsim.fit.line.SlopeInterceptLineFitter;
+import java.util.ArrayList;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseTrackerHit;
+import java.util.Collections;
+import java.util.*;
/**
- *
+ * Fit a helix to a set of space points. First, a circle is fit to the x-y coordinates.
+ * A straight-line fit is then performed on the s-z coordinates.
+ *
+ * The r-phi and z coordinate measurements are assumed to be uncorrelated. A block
+ * diagonal covariance matrix is formed from the results of the circle and s-z fits,
+ * ignoring any correlations between these fits.
+ *
+ * The resulting track parameters follow the "L3 Convention" adopted by org.lcsim.
+ *
+ * Modified July 2007 by Lori Stevens
+ * Modified August 2006 by Richard Partridge
* @author Norman Graf
*/
-public class HelicalTrackFitter
+
+public class HelicalTrackFitter
{
private CircleFitter _cfitter = new CircleFitter();
private SlopeInterceptLineFitter _lfitter = new SlopeInterceptLineFitter();
@@ -30,69 +47,135 @@
private double[] _circleParameters = new double[3];
- /** Creates a new instance of HelicalTrackFitter */
+ private ArrayList<TrackerHit> hit = new ArrayList<TrackerHit>();
+ /**
+ * Creates a new instance of HelicalTrackFitter.
+ */
public HelicalTrackFitter()
{
- }
-
- public boolean fit(double[] x, double[] y, double[] z, double[] drphi, double[] dz, int np)
- {
- // fit a circle in the x-y plane
+ }
+ /**
+ * Reorder list of TrackerHits from closest to furthest from origin in z.
+ * Calculate weights for circle and line fits.
+ * @param hits list of TrackerHits
+ * @return True for successful fit
+ */
+ public boolean fit(List<TrackerHit> hits)
+ {
+ int np = hits.size();
+ double[] wrphi = new double[np];
+ double[] wz = new double[np];
+
+ order(hits);
+
+ for(int i=0; i<np; ++i){
+ //(1/(dxdx + dydy))
+ wrphi[i] = 1/(hits.get(i).getCovMatrix()[0]+hits.get(i).getCovMatrix()[2]);
+ //(1/dzdz)
+ wz[i] = 1/hits.get(i).getCovMatrix()[5];
+ }
+ return fitting(hits, wrphi, wz);
+ }
+ /**
+ * Convert hits in Cartesian coordinates to BaseTrackerHits for helix fitting.
+ * Order hits from closest to furthest from origin in z. Calculate weights
+ * for circle and line fits.
+ * @param x array of x coordinates
+ * @param y array of y coordinates
+ * @param z array of z coordinates
+ * @param drphi error in r-phi hit position
+ * @param dz error in z coordinate. A negative value indicates a barrel hit.
+ * @param np number of points
+ * @return true for successful fit
+ */
+ public boolean fit(double[] x, double[] y, double[] z, double[] drphi, double[] dz, int np)
+ {
+ for(int i=0; i<np; i++)
+ {
+ BaseTrackerHit bth = new BaseTrackerHit();
+ double[] pos = {x[i], y[i], z[i]};
+ bth.setPosition(pos);
+ if (dz[i]>0) bth.setType(0);
+ if (dz[i]<0) bth.setType(1);
+
+ hit.add(bth);
+ }
+ order(hit);
+
double[] wrphi = new double[np];
+ double[] wz = new double[np];
+
for(int i=0; i<np; ++i)
{
- wrphi[i] = 1/(drphi[i]*drphi[i]);
+ wrphi[i] = 1/(drphi[i]*drphi[i]);
+ wz[i] = 1/(dz[i]*dz[i]);
+ }
+ return fitting(hit, wrphi, wz);
+ }
+ /**
+ * Fit a helix to set of TrackerHits.
+ * @param hits list of TrackerHits
+ * @param wrphi array of weights for circle fitter
+ * @param wz array of weights for line fitter
+ * @return true for successful fit
+ */
+ public boolean fitting(List<TrackerHit> hits, double[] wrphi, double[] wz)
+ {
+ int np = hits.size();
+ double[] x = new double[np];
+ double[] y = new double[np];
+ double[] z = new double[np];
+
+ for(int i=0; i<np; i++)
+ {
+ x[i] = hits.get(i).getPosition()[0];
+ y[i] = hits.get(i).getPosition()[1];
+ z[i] = hits.get(i).getPosition()[2];
}
- boolean success = _cfitter.fit(x, y, wrphi, np);
+
+ boolean success = _cfitter.fit(x, y, wrphi, np);
if(!success) return false;
CircleFit cfit = _cfitter.getfit();
double radius = 1./cfit.curvature();
+// System.out.println("radius: "+radius);
double phi0 = cfit.phi();
- double dca = cfit.dca();
+ double dca = -cfit.dca();
double xc = (radius-dca)*sin(phi0);
- double yc = (radius-dca)*cos(phi0);
+ double yc = -(radius-dca)*cos(phi0);
_circleParameters[0] = xc;
_circleParameters[1] = yc;
_circleParameters[2] = radius;
// fit a line in the s-z plane
// assumes points are in increasing order
- //TODO check whether hit has z information
- // TODO double-check path length calculation
double[] s = new double[np];
- double[] wz = new double[np];
- for(int i=0; i<np; ++i)
+ double[] sfit = new double[np];
+ double[] zfit = new double[np];
+ double[] wrphibh = new double[np];
+ int nz = 0;
+ for(int i=0; i<np; i++)
{
- s[i] = s(radius, xc, yc, x[i], y[i], 0., 0.);
- wz[i] = 1/(dz[i]*dz[i]);
- }
+ if (i==0) {
+ s[0] = s(radius, xc, yc, x[0], y[0], 0., 0.);
+ }
+ else {
+ s[i] = s(radius, xc, yc, x[i], y[i], x[i-1], y[i-1]) + s[i-1];
+ }
+ //type 0 = 3D vertex hit, can only do line fit for vertex hits
+ if (hits.get(i).getType()==0)
+ {
+ sfit[nz] = s[i];
+ zfit[nz] = z[i];
+ nz++;
+ }
+ }
- success = _lfitter.fit(s, z, dz, np);
+ success = _lfitter.fit(sfit, zfit, wz, nz);
if(!success) return false;
SlopeInterceptLineFit lfit = _lfitter.getFit();
- // TODO convert fit results to track parameters and covariance matrix.
- // now have the fits, need to assemble the track parameters...
+ double[] pars = new double[5];
- // see org.lcsim.event.Track.Parameters
- // The track parameters for LCD are defined as follows
- // <table>
- // <tr><th>Index</th><th>Meaning</th></tr>
- // <tr><td> 0 </td><td> d0 = XY impact parameter </td><tr>
- // <tr><td> 1 </td><td> phi0 </td><tr> </td><tr>
- // <tr><td> 2 </td><td> omega = 1/curv.radius (negative for negative tracks) </td><tr>
- // <tr><td> 3 </td><td> z0 = z of track origin (z impact parameter) </td><tr>
- // <tr><td> 4 </td><td> s = tan lambda </td><tr>
- // </table>
-
- // tanlambda = line fit slope
- // z0 = line fit intercept
- //TODO should shift s-z line fit to centroid of x-y fit and introduce covariance terms to be correct.
-
-
-// System.out.println("circle fit: "+cfit);
-// System.out.println("line fit: "+lfit);;
- double[] pars = new double[5];
SymmetricMatrix cov = new SymmetricMatrix(5);
double[] chisq = new double[2];
int[] ndf = new int[2];
@@ -101,11 +184,10 @@
chisq[1] = lfit.chisquared();
ndf[1] = lfit.ndf();
- // TODO fix this so CircleFit returns ndf
- ndf[0] = lfit.ndf();
+ ndf[0] = np - 3;
- // d0, xy impact parameter.
- pars[0] = cfit.dca();
+ // d0, xy impact parameter. Note: - sign due to opposite sign convention in CircleFitter
+ pars[0] = -cfit.dca();
// phi0
pars[1] = cfit.phi();
// omega signed curvature
@@ -114,34 +196,80 @@
pars[3] = lfit.intercept();
// tan lambda
pars[4] = lfit.slope();
-
- // TODO fill in covariance matrix
-// for(int i=0; i<5; ++i)
-// {
-// System.out.println("pars["+i+"]= "+pars[i]);
-// }
+
+ // fill in covariance matrix
+ cov.setElement(0,0, cfit.cov()[0]);
+ cov.setElement(0,1,-cfit.cov()[1]); // Need - sign due to sign convention in CircleFitter
+ cov.setElement(1,1, cfit.cov()[2]);
+ cov.setElement(0,2,-cfit.cov()[3]); // Need - sign due to sign convention in CircleFitter
+ cov.setElement(1,2, cfit.cov()[4]);
+ cov.setElement(2,2, cfit.cov()[5]);
+ //cov[0][3] = 0.;
+ //cov[1][3] = 0.;
+ //cov[2][3] = 0.;
+ cov.setElement(3,3, lfit.interceptUncertainty());
+ //cov[0][4] = 0.;
+ //cov[1][4] = 0.;
+ //cov[2][4] = 0.;
+ cov.setElement(3,4, lfit.covariance());
+ cov.setElement(4,4, lfit.slopeUncertainty());
_fit = new HelicalTrackFit(pars, cov, chisq, ndf);
_fit.setCircleParameters(_circleParameters);
return true;
}
-
+ /**
+ * Get the results of the most recent fit.
+ * @return HelicalTrackFit corresponding to most recent fit
+ */
public HelicalTrackFit getFit()
{
return _fit;
}
-
+ /**
+ * Get the circle paramaters (xc, yc, R) from the most recent fit.
+ * @return Circle parameters (xc, yc, R)
+ */
public double[] getCircleParameters()
{
return _circleParameters;
- }
-
+ }
double s(double radius, double xcenter, double ycenter, double x1, double y1, double x2, double y2)
{
double phi1 = atan2( (y1-ycenter), (x1-xcenter) );
double phi2 = atan2( (y2-ycenter), (x2-xcenter) );
double dphi = abs(phi1-phi2);
if(dphi>PI) dphi = 2.*PI-dphi;
- return abs(radius*dphi);
+ return abs(radius)*dphi;
}
-
-}
\ No newline at end of file
+ /**
+ * Reorder TrackerHits from closest to furthest from origin in z.
+ * @param hits list of TrackerHits
+ */
+ public void order(List<TrackerHit> hits)
+ {
+ Comparator<TrackerHit> z = new Comp();
+ Collections.sort(hits, z);
+ double[] hit1 = hits.get(0).getPosition();
+ double[] hit2 = hits.get(hits.size()-1).getPosition();
+ double dist1 = Math.sqrt(hit1[0]*hit1[0]+hit1[1]*hit1[1]+hit1[2]*hit1[2]);
+ double dist2 = Math.sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]+hit2[2]*hit2[2]);
+ if(dist1 > dist2){
+ Collections.reverse(hits);
+ }
+ }
+}
+class Comp implements Comparator<TrackerHit>
+{
+ public int compare(TrackerHit hit1, TrackerHit hit2)
+ {
+ double z1 = hit1.getPosition()[2];
+ double z2 = hit2.getPosition()[2];
+
+ if (z1 < z2)
+ return -1;
+ else if (z1 > z2)
+ return 1;
+ else
+ return 0;
+ }
+}
