lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.10 -r1.11
--- HelicalTrackFitter.java 31 Aug 2007 22:32:31 -0000 1.10
+++ HelicalTrackFitter.java 6 Nov 2007 18:31:16 -0000 1.11
@@ -4,7 +4,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.10 2007/08/31 22:32:31 lstevens Exp $
+ * $Id: HelicalTrackFitter.java,v 1.11 2007/11/06 18:31:16 tknelson Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
@@ -18,8 +18,8 @@
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 java.util.ArrayList;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseTrackerHit;
import java.util.Collections;
import java.util.List;
@@ -32,19 +32,19 @@
/**
* 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();
@@ -60,37 +60,55 @@
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.
+ * 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];
+ 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);
- }
+ _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);
+ }
/**
* Convert hits in Cartesian coordinates to BaseTrackerHits for helix fitting.
- * Order hits from closest to furthest from origin in z. Calculate weights
+ * 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
@@ -100,8 +118,8 @@
* @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)
- {
+ 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();
@@ -115,15 +133,15 @@
_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]);
+ wrphi[i] = 1/(drphi[i]*drphi[i]);
}
return fitting(_hitsord, wrphi, dz);
- }
+ }
/**
* Fit a helix to set of TrackerHits.
* @param hits list of TrackerHits
@@ -132,14 +150,16 @@
* @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){
+ 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
+ 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);
@@ -157,9 +177,10 @@
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){
+ }
+ boolean success = _cfitter.fit(x, y, wrphi, np);
+ if(!success)
+ {
// _aida.cloud1D("failed circle fit").fill(1);
return false;
}
@@ -193,141 +214,153 @@
int nb = 0;
for(int i=0; i<np; i++)
{
- if (i==0) {
+ if (i==0)
+ {
s[0] = s(radius, xc, yc, x[0], y[0], x0, y0);
// System.out.println("s: "+s[0]);
}
- else {
+ 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 (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++;
- }
+ 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];
+ double[] pars = new double[5];
double[] chisq = new double[2];
int[] ndf = new int[2];
chisq[0] = cfit.chisq();
ndf[0] = np - 3;
-
+
// 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();
-
- if(n3D >= 2){
- success = _lfitter.fit(sfit, zfit, v_wz, nz);
- if(!success){
- _aida.cloud1D("failed line fit").fill(1);
- return false;
- }
- if(success) _aida.cloud1D("successful line fit").fill(1);
- SlopeInterceptLineFit lfit = _lfitter.getFit();
-
- chisq[1] = lfit.chisquared();
-
- ndf[1] = lfit.ndf();
-
- // z0
- pars[3] = lfit.intercept();
- // tan lambda
- pars[4] = lfit.slope();
-
- //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());
-
- 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(n3D >= 2)
+ {
+ success = _lfitter.fit(sfit, zfit, v_wz, nz);
+ if(!success)
+ {
+ _aida.cloud1D("failed line fit").fill(1);
+ return false;
+ }
+ if(success) _aida.cloud1D("successful line fit").fill(1);
+
+ SlopeInterceptLineFit lfit = _lfitter.getFit();
+
+ chisq[1] = lfit.chisquared();
+
+ ndf[1] = lfit.ndf();
+
+ // z0
+ pars[3] = lfit.intercept();
+ // tan lambda
+ pars[4] = lfit.slope();
+
+ //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());
+
+ 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);
+ 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);
+ 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;
+ return false;
+ }
}
}
}
- }
- 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];
-
- success = _zfitter.fit(bs, zmin, zmax);
- if(!success){
- // _aida.cloud1D("failed zseg w/1 vtx hit").fill(1);
- return false;
- }
+ 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];
+
+ success = _zfitter.fit(bs, zmin, zmax);
+ if(!success)
+ {
+ // _aida.cloud1D("failed zseg w/1 vtx hit").fill(1);
+ return false;
+ }
// 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(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");
- }
- 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();
- }
+ }
+ 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);
-
+
return true;
}
- /**
+ /**
* Get the results of the most recent fit.
* @return HelicalTrackFit corresponding to most recent fit
*/
@@ -335,14 +368,14 @@
{
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) );
@@ -363,10 +396,11 @@
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);
- }
- }
+ if(dist1 > dist2)
+ {
+ Collections.reverse(hitslist);
+ }
+ }
private ArrayList<TrackerHit> ordhits(List<TrackerHit> hits)
{
ArrayList<TrackerHit> hitlist = new ArrayList<TrackerHit>();
@@ -375,52 +409,44 @@
vhits.clear();
bhits.clear();
hitlist.clear();
- for(TrackerHit hit: hits){
- if(hit.getType()==0) vhits.add(hit);
- if(hit.getType()==1) bhits.add(hit);
+ for(TrackerHit hit: hits)
+ {
+ if(hit.getType()==0) vhits.add(hit);
+ if(hit.getType()==1) bhits.add(hit);
}
- if(!vhits.isEmpty()){
+ if(!vhits.isEmpty())
+ {
order(vhits);
hitlist.addAll(vhits);
- }
- if(!bhits.isEmpty()) { //need to order according to s
+ }
+ 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};
-
- return modInfo;
- }
-}
-class Comp implements Comparator<TrackerHit>
-{
- public int compare(TrackerHit hit1, TrackerHit hit2)
+ public double[] moduleInfo(double z)
{
- double z1 = hit1.getPosition()[2];
- double z2 = hit2.getPosition()[2];
+ 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};
- if (z1 < z2)
- return -1;
- else if (z1 > z2)
- return 1;
- else
- return 0;
+ return modInfo;
}
}