lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.4 -r1.5
--- HelicalTrackFit.java 23 Oct 2006 19:42:36 -0000 1.4
+++ HelicalTrackFit.java 11 Dec 2007 20:35:04 -0000 1.5
@@ -3,57 +3,67 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFit.java,v 1.4 2006/10/23 19:42:36 tonyj Exp $
+ * $Id: HelicalTrackFit.java,v 1.5 2007/12/11 20:35:04 partridge Exp $
*/
package org.lcsim.fit.helicaltrack;
import hep.physics.matrix.SymmetricMatrix;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.lcsim.fit.circle.CircleFit;
+import org.lcsim.fit.line.SlopeInterceptLineFit;
+import org.lcsim.fit.zsegment.ZSegmentFit;
+
/**
* Represents the result of a fit to a helical track.
* @author Norman Graf
*/
-public class HelicalTrackFit
-{
+public class HelicalTrackFit {
// fit independently to a circle in x-y and a line in s-z
// first is circle, second is line
private double[] _chisq = new double[2];
private int[] _ndf = new int[2];
private double[] _parameters;
private SymmetricMatrix _covmatrix;
- private double[] _circleParameters = new double[3];
+ private Map<HelicalTrackHit, Double> _smap;
/** Creates a new instance of HelicalTrackFit */
- public HelicalTrackFit(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf)
- {
+ public HelicalTrackFit(double[] pars, SymmetricMatrix cov, double[] chisq, int[] ndf) {
+ this(pars, cov, chisq, ndf, new HashMap<HelicalTrackHit, Double>());
+ }
+
+ public HelicalTrackFit(double[]pars, SymmetricMatrix cov, double[] chisq, int[] ndf,
+ Map<HelicalTrackHit, Double> smap) {
_parameters = pars;
_covmatrix = cov;
_chisq = chisq;
_ndf = ndf;
+ _smap = smap;
}
- public double[] parameters()
- {
+ public double[] parameters() {
return _parameters;
}
- public SymmetricMatrix covariance()
- {
+ public SymmetricMatrix covariance() {
return _covmatrix;
}
- public double[] chisq()
- {
+ public double[] chisq() {
return _chisq;
}
- public int[] ndf()
- {
+ public int[] ndf() {
return _ndf;
}
- public String toString()
- {
+ public Map<HelicalTrackHit, Double> pathmap() {
+ return _smap;
+ }
+
+ public String toString() {
StringBuffer sb = new StringBuffer("HelicalTrackFit: \n");
sb.append("d0= "+_parameters[0]+"\n");
@@ -65,14 +75,6 @@
return sb.toString();
}
- public void setCircleParameters(double[] pars)
- {
- System.arraycopy(pars,0,_circleParameters,0,3);
- }
- public double[] circleParameters()
- {
- return _circleParameters;
- }
}
lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.11 -r1.12
--- HelicalTrackFitter.java 6 Nov 2007 18:31:16 -0000 1.11
+++ HelicalTrackFitter.java 11 Dec 2007 20:35:04 -0000 1.12
@@ -4,34 +4,30 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.11 2007/11/06 18:31:16 tknelson Exp $
+ * $Id: HelicalTrackFitter.java,v 1.12 2007/12/11 20:35:04 partridge Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
-import static java.lang.Math.atan2;
-import static java.lang.Math.abs;
-import static java.lang.Math.PI;
-import static java.lang.Math.cos;
-import static java.lang.Math.sin;
+import hep.physics.vec.BasicHep3Vector;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+// import org.lcsim.event.TrackerHit;
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.List;
-import java.util.Comparator;
-import java.util.Arrays;
import org.lcsim.fit.zsegment.ZSegmentFit;
import org.lcsim.fit.zsegment.ZSegmentFitter;
-import hep.aida.*;
-import org.lcsim.util.aida.AIDA;
+
/**
* 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.
+ * A straight-line fit is then performed on the s-z coordinates. If there are two few
+ * 3D hits to perform the s-z fit, the ZSegmentFitter is used to find the s-z parameters.
*
* 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,
@@ -39,77 +35,28 @@
*
* The resulting track parameters follow the "L3 Convention" adopted by org.lcsim.
*
+ * Modified November 2007 by Richard Partridge
* 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();
private ZSegmentFitter _zfitter = new ZSegmentFitter();
private HelicalTrackFit _fit;
- private double[] _circleParameters = new double[3];
-
- private ArrayList<TrackerHit> _hits = new ArrayList<TrackerHit>();
- private ArrayList<TrackerHit> _hitsord = new ArrayList<TrackerHit>();
- private ArrayList<TrackerHit> _3D = new ArrayList<TrackerHit>();
- private ArrayList<TrackerHit> _2D = new ArrayList<TrackerHit>();
-
- AIDA _aida = AIDA.defaultInstance();
-
- 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;
- }
- }
-
/**
* Creates a new instance of HelicalTrackFitter.
*/
- public HelicalTrackFitter()
- {
- }
- /**
- * 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];
- _hitsord.clear();
- _hitsord = ordhits(hits);
-
- for(int i=0; i<np; ++i)
- {
- //(1/(dxdx + dydy))
- wrphi[i] = 1/(_hitsord.get(i).getCovMatrix()[0]+_hitsord.get(i).getCovMatrix()[2]);
- //dz
- wz[i] = Math.sqrt(_hitsord.get(i).getCovMatrix()[5]);
- }
- return fitting(_hitsord, wrphi, wz);
+ public HelicalTrackFitter() {
}
+
/**
- * 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.
+ * Original constructor for HelicalTrackFitter kept for backwards compatibility.
+ * Not recommened for new code, may be deprecated in the future.
* @param x array of x coordinates
* @param y array of y coordinates
* @param z array of z coordinates
@@ -118,335 +65,229 @@
* @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);
- //flags that hit is 2D
- if (dz[i]<0) bth.setType(1);
- _hits.add(bth);
- }
- _hitsord.clear();
- _hitsord = ordhits(_hits);
- _hits.clear();
-
- double[] wrphi = new double[np];
-
- for(int i=0; i<np; ++i)
- {
- wrphi[i] = 1/(drphi[i]*drphi[i]);
+ public boolean 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++) {
+ if (dz[i] > 0.) {
+ hitcol.add(new HelicalTrack3DHit(x[i], y[i], z[i], drphi[i], dz[i]));
+ } else {
+ double zmin = z[i] - Math.abs(dz[i]);
+ double zmax = z[i] + Math.abs(dz[i]);
+ hitcol.add(new HelicalTrack2DHit(x[i], y[i], z[i], drphi[i], zmin, zmax));
+ }
}
- return fitting(_hitsord, wrphi, dz);
+ Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
+ return fit(hitcol, msmap);
}
- /**
- * 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)
- {
- SymmetricMatrix cov = new SymmetricMatrix(5);
- _3D.clear();
- _2D.clear();
- for(TrackerHit hit: hits)
- {
- if(hit.getType()==0) _3D.add(hit);
- if(hit.getType()==1)
- {
- //added this line for problem with cheater, don't want more
- //than 5 outer tracker barrel hits per track
-// if(_2D.size()==5) continue;
- _2D.add(hit);
- }
+
+ public boolean fit(List<HelicalTrackHit> hitcol) {
+ Map<HelicalTrackHit, MultipleScatter> msmap = new HashMap<HelicalTrackHit, MultipleScatter>();
+ return fit(hitcol, msmap);
}
- int np = hits.size();
- int n3D = _3D.size();
- int n2D = _2D.size();
- double[] x = new double[np];
- double[] y = new double[np];
- double[] z = new double[np];
+
+ public boolean fit(List<HelicalTrackHit> hitcol, Map<HelicalTrackHit, MultipleScatter> msmap) {
- 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);
- if(!success)
- {
-// _aida.cloud1D("failed circle fit").fill(1);
- return false;
+ // Create lists for the various types of hits
+ List<HelicalTrackHit> circle_hits = new ArrayList<HelicalTrackHit>();
+ List<HelicalTrackHit> pixel_hits = new ArrayList<HelicalTrackHit>();
+ List<HelicalTrackHit> strip_hits = new ArrayList<HelicalTrackHit>();
+ // 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);
+ // Pixel hits
+ if (hit instanceof HelicalTrack3DHit) pixel_hits.add(hit);
+ // Strip hits
+ if (hit instanceof HelicalTrack2DHit) strip_hits.add(hit);
}
-// if(success) _aida.cloud1D("successful circle fit").fill(1);
-
- CircleFit cfit = _cfitter.getfit();
- double radius = 1./cfit.curvature();
- double phi0 = cfit.phi();
- double dca = -cfit.dca();
- double xc = (radius-dca)*sin(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
- //need to add 1 to array size if there is 1 3D hit
- int add3D=0;
- if(n3D==1) add3D=1;
- double[] s = new double[np];
- double[] sfit = new double[n3D];
- double[] zfit = new double[n3D];
- double[] v_wz = new double[n3D];
- double[] bs = new double[n2D+add3D];
- double[] zmin = new double[n2D+add3D];
- double[] zmax = new double[n2D+add3D];
- double x0 = dca*Math.sin(phi0);
- double y0 = dca*Math.cos(phi0);
- int nz = 0;
- int nb = 0;
- for(int i=0; i<np; i++)
- {
- if (i==0)
- {
- s[0] = s(radius, xc, yc, x[0], y[0], x0, y0);
-// System.out.println("s: "+s[0]);
- }
- else
- {
- s[i] = s(radius, xc, yc, x[i], y[i], x[i-1], y[i-1]) + s[i-1];
-// System.out.println(" next s: "+s[i]+" last s: "+s[i-1]);
- }
- if (hits.get(i).getType()==0)
- {
- sfit[nz] = s[i];
- zfit[nz] = z[i];
- v_wz[nz] = wz[i];
- nz++;
- }
- if (hits.get(i).getType()==1)
- {
- //added line for problem with cheater
-// if (nb==5) continue;
- bs[nb] = s(radius, xc, yc, x[i], y[i], x0, y0);
- zmin[nb] = moduleInfo(z[i])[0];
- zmax[nb] = moduleInfo(z[i])[1];
- nb++;
- }
- }
- double[] pars = new double[5];
+ // Create the arrays that will hold the fit output
double[] chisq = new double[2];
- int[] ndf = new int[2];
-
- chisq[0] = cfit.chisq();
-
- ndf[0] = np - 3;
+ int[] ndof = new int[2];
+ double[] par = new double[5];
+ SymmetricMatrix cov = new SymmetricMatrix(5);
- // 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
- pars[2] = cfit.curvature();
+ // Setup for doing the circle fit
+ int nc = circle_hits.size();
+ if (nc < 3) return false;
+ double[] x = new double[nc];
+ double[] y = new double[nc];
+ double[] wrphi = new double[nc];
+
+ // Store the hit coordinates and weights in arrays for the circle fitter
+ for(int i=0; i<nc; i++) {
+ HelicalTrackHit hit = circle_hits.get(i);
+ x[i] = hit.x();
+ y[i] = hit.y();
+ double drphi_ms = 0.;
+ if (msmap.containsKey(hit)) drphi_ms = msmap.get(hit).drphi();
+ 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
+ boolean success = _cfitter.fit(x, y, wrphi, nc);
+ if(!success) return false;
- if(n3D >= 2)
- {
- success = _lfitter.fit(sfit, zfit, v_wz, nz);
- if(!success)
- {
- _aida.cloud1D("failed line fit").fill(1);
- return false;
+ // Store the circle fit output
+ CircleFit cfit = _cfitter.getfit();
+ chisq[0] = cfit.chisq();
+ ndof[0] = nc - 3;
+ par[0] = -cfit.dca();
+ par[1] = cfit.phi();
+ par[2] = cfit.curvature();
+ cov.setElement(0, 0, cfit.cov()[0]);
+ cov.setElement(0, 1, cfit.cov()[1]);
+ cov.setElement(1, 1, cfit.cov()[2]);
+ 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);
+
+ // Check if we have enough pixel hits to do a straight-line fit of s(z)
+ int npix = pixel_hits.size();
+ if (npix > 1) {
+
+ // Setup for the line fit
+ double[] s = new double[npix];
+ double[] z = new double[npix];
+ double[] dz = new double[npix];
+
+ // Store the coordinate and errors for the line fit
+ for(int i=0; i<npix; i++) {
+ HelicalTrack3DHit hit = (HelicalTrack3DHit) pixel_hits.get(i);
+ z[i] = hit.z();
+ double dz_ms = 0.;
+ if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
+ 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]);
}
- if(success) _aida.cloud1D("successful line fit").fill(1);
- SlopeInterceptLineFit lfit = _lfitter.getFit();
+ // 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
+ SlopeInterceptLineFit lfit = _lfitter.getFit();
chisq[1] = lfit.chisquared();
+ ndof[1] = npix - 2;
+ par[3] = lfit.intercept();
+ par[4] = lfit.slope();
+ cov.setElement(3, 3, Math.pow(lfit.interceptUncertainty(),2));
+ cov.setElement(3, 4, lfit.covariance());
+ cov.setElement(4, 4, Math.pow(lfit.slopeUncertainty(),2));
- ndf[1] = lfit.ndf();
-
- // z0
- pars[3] = lfit.intercept();
- // tan lambda
- pars[4] = lfit.slope();
+ } else {
- //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());
+ // Not enough pixel hits for a line fit, try a ZSegment fit
- if(n2D >= 1)
- {
- int range = n2D;
- //had this line to prevent more than 5 outer tracker barrel hits per track
- //but put in 2 lines earlier in code that does same thing; for cheater problem
-// if (n2D >= 5) range = 5;
- for(int i=0; i<range; i++)
- {
- //z = z0+s*tanlambda
- double zpred = lfit.intercept() + lfit.slope()*bs[i];
- double zerr = 3.*Math.sqrt(Math.abs(cov.e(3,3)+2*bs[i]*cov.e(4,3)+cov.e(4,4)));
- if (!(zmin[i] <= (zpred+zerr)) || !((zpred-zerr) <= zmax[i]))
- {
-// _aida.cloud1D("failed 2D zcheck").fill(1);
-// _aida.cloud1D("#hit that failed").fill(i+1);
-// _aida.cloud1D("# barrel hits for failed 2D check").fill(_bhits.size());
- double mindiff = Math.abs(zmax[i]+zerr-zpred);
- double maxdiff = Math.abs(zmin[i]-(zpred+zerr));
- double diff=mindiff;
- if(maxdiff<mindiff) diff=maxdiff;
- //tells how far outside zmin or zmax the predicted z lies
- _aida.cloud1D("z difference").fill(diff);
-
- return false;
- }
- }
+ // If we have one pixel hit, turn it into a pseudo strip hit
+ if (npix == 1) {
+ HelicalTrack3DHit hit = (HelicalTrack3DHit) pixel_hits.get(0);
+ double dz_ms = 0.;
+ if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
+ double dz = Math.sqrt(Math.pow(hit.dz(),2)+Math.pow(dz_ms,2));
+ double zmin = hit.z() - dz;
+ double zmax = hit.z() + dz;
+ strip_hits.add(new HelicalTrack2DHit(hit.x(), hit.y(), hit.z(), hit.drphi(), zmin, zmax));
}
- }
- if(n3D <=1 && n2D >=2)
- {
- if(n3D==1)
- {
- double length = 3.*v_wz[0];
- zmin[n2D] = _3D.get(0).getPosition()[2]-length;
- zmax[n2D] = _3D.get(0).getPosition()[2]+length;
- //3D segment gets put at end of array
- bs[n2D] = sfit[0];
+
+ // 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];
- success = _zfitter.fit(bs, zmin, zmax);
- if(!success)
- {
- // _aida.cloud1D("failed zseg w/1 vtx hit").fill(1);
- return false;
+ // Store the coordinates and errors for the ZSegment fit
+ for(int i=0; i<npix; i++) {
+ HelicalTrack2DHit hit = (HelicalTrack2DHit) strip_hits.get(i);
+ s[i] = smap.get(hit);
+ zmin[i] = hit.zmin();
+ zmax[i] = hit.zmax();
}
-// if(success) _aida.cloud1D("successful zseg fit w/1 vtx hit");
- }
- if(n3D==0 && n2D >=3)
- {
- success = _zfitter.fit(bs, zmin, zmax);
- if(!success)
- {
-// _aida.cloud1D("failed zseg").fill(1);
- return false;
- }
-// if(success) _aida.cloud1D("successful zseg fit");
+
+ // Call the ZSegment fitter and check for success
+ success = _zfitter.fit(s, zmin, zmax);
+ if(!success) return false;
+
+ // Store the ZSegment fit output
+ ZSegmentFit zfit = _zfitter.getFit();
+ chisq[1] = 0.;
+ ndof[1] = 0;
+ par[3] = zfit.getCentroid()[0];
+ par[4] = zfit.getCentroid()[1];
+ cov.setElement(3, 3, zfit.getCovariance().e(0, 0));
+ cov.setElement(3, 4, zfit.getCovariance().e(0,1));
+ cov.setElement(4, 4, zfit.getCovariance().e(1, 1));
}
- ZSegmentFit _zfit = _zfitter.getFit();
- //chisq doesn't make sense for zfitter, set to zero
- chisq[1] = 0.;
- //ndf doesn't make sens for zfitter, set to zero
- ndf[1] = 0;
- // z0
- pars[3] = _zfit.getCentroid()[0];
- // tan lambda
- pars[4] = _zfit.getCentroid()[1];
-
- cov = _zfit.getCovariance();
}
- _fit = new HelicalTrackFit(pars, cov, chisq, ndf);
- _fit.setCircleParameters(_circleParameters);
+ // Create the HelicalTrackFit for this helix and exit
+ _fit = new HelicalTrackFit(par, cov, chisq, ndof, smap);
return true;
}
+
/**
* Get the results of the most recent fit.
* @return HelicalTrackFit corresponding to most recent fit
*/
- public HelicalTrackFit getFit()
- {
+ 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;
- }
- /**
- * Reorder TrackerHits from closest to furthest from origin in z.
- * @param hits list of TrackerHits
+ * Find the arc lengths for a list of hits
*/
- private void order(List<TrackerHit> hitslist)
- {
- Comparator<TrackerHit> z = new Comp();
- Collections.sort(hitslist, z);
- double[] hit1 = hitslist.get(0).getPosition();
- double[] hit2 = hitslist.get(hitslist.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(hitslist);
- }
- }
- private ArrayList<TrackerHit> ordhits(List<TrackerHit> hits)
- {
- ArrayList<TrackerHit> hitlist = new ArrayList<TrackerHit>();
- ArrayList<TrackerHit> vhits = new ArrayList<TrackerHit>();
- ArrayList<TrackerHit> bhits = new ArrayList<TrackerHit>();
- vhits.clear();
- bhits.clear();
- hitlist.clear();
- for(TrackerHit hit: hits)
- {
- if(hit.getType()==0) vhits.add(hit);
- if(hit.getType()==1) bhits.add(hit);
- }
- if(!vhits.isEmpty())
- {
- order(vhits);
- hitlist.addAll(vhits);
- }
- if(!bhits.isEmpty())
- { //need to order according to s
- order(bhits); //should I put this line in?
- hitlist.addAll(bhits);
- }
- return hitlist;
- }
- public double[] moduleInfo(double z)
- {
- double segs = 100., minz=0, maxz=0;
- //first positive z segment has zmin=0
- if (z>=0)
- {
- minz = ((int)(z/segs))*segs;
- maxz = minz+segs;
- }
- //neg numbers round in positive direction. ex. -2.2 becomes -2
- else if (z<0)
- {
- maxz = ((int)(z/segs))*segs;
- minz = maxz-segs;
- if (maxz==minz)
- {
- maxz += segs;
- }
- }
- double [ ] modInfo = {minz, maxz};
+ public Map<HelicalTrackHit, Double> getPathLengths(double[] par, List<HelicalTrackHit> hits) {
- return modInfo;
+ // 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.
+ int nhits = hits.size();
+ Collections.sort(hits);
+ 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.;
+ for(int i=0; i<nhits; 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;
+ }
+ return smap;
}
-}
+}
\ No newline at end of file