lcsim/src/org/lcsim/fit/helicaltrack
diff -u -r1.9 -r1.10
--- HelicalTrackFitter.java 20 Aug 2007 20:49:35 -0000 1.9
+++ HelicalTrackFitter.java 31 Aug 2007 22:32:31 -0000 1.10
@@ -4,7 +4,7 @@
*
* Created on March 25, 2006, 6:11 PM
*
- * $Id: HelicalTrackFitter.java,v 1.9 2007/08/20 20:49:35 lstevens Exp $
+ * $Id: HelicalTrackFitter.java,v 1.10 2007/08/31 22:32:31 lstevens Exp $
*/
import hep.physics.matrix.SymmetricMatrix;
@@ -18,64 +18,79 @@
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.*;
+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.
- *
+ *
* 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();
+ private ZSegmentFitter _zfitter = new ZSegmentFitter();
private HelicalTrackFit _fit;
private double[] _circleParameters = new double[3];
- private ArrayList<TrackerHit> hit = new ArrayList<TrackerHit>();
+ 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();
/**
* Creates a new instance of HelicalTrackFitter.
*/
- public 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) {
+ */
+ 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] = Math.sqrt(hits.get(i).getCovMatrix()[5]);
- }
- return fitting(hits, wrphi, wz);
- }
+ 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);
+ }
/**
* 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
@@ -85,28 +100,30 @@
* @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) {
-
- hit.clear();
-
- for(int i=0; i<np; i++) {
+ 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);
-
- hit.add(bth);
+ _hits.add(bth);
}
- order(hit);
-
+ _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]);
+ for(int i=0; i<np; ++i)
+ {
+ wrphi[i] = 1/(drphi[i]*drphi[i]);
}
- return fitting(hit, wrphi, dz);
- }
+ return fitting(_hitsord, wrphi, dz);
+ }
/**
* Fit a helix to set of TrackerHits.
* @param hits list of TrackerHits
@@ -114,24 +131,42 @@
* @param wz array of weights for line fitter
* @return true for successful fit
*/
- public boolean fitting(List<TrackerHit> hits, double[] wrphi, double[] wz) {
+ 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);
+ }
+ }
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];
- for(int i=0; i<np; i++) {
+ 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;
}
-
- boolean success = _cfitter.fit(x, y, wrphi, np);
- if(!success) return false;
+// if(success) _aida.cloud1D("successful circle fit").fill(1);
CircleFit cfit = _cfitter.getfit();
double radius = 1./cfit.curvature();
-// System.out.println("radius: "+radius);
double phi0 = cfit.phi();
double dca = -cfit.dca();
double xc = (radius-dca)*sin(phi0);
@@ -139,61 +174,84 @@
_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[np];
- double[] zfit = new double[np];
- double[] wrphibh = 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;
- for(int i=0; i<np; 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++;
+ 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++;
+ }
}
-
- success = _lfitter.fit(sfit, zfit, wz, nz);
- if(!success) return false;
- SlopeInterceptLineFit lfit = _lfitter.getFit();
-
- double[] pars = new double[5];
-
- SymmetricMatrix cov = new SymmetricMatrix(5);
+ double[] pars = new double[5];
double[] chisq = new double[2];
int[] ndf = new int[2];
chisq[0] = cfit.chisq();
- chisq[1] = lfit.chisquared();
-
- ndf[1] = lfit.ndf();
- ndf[0] = np - 3;
+ 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();
-
- // 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.;
@@ -203,25 +261,90 @@
//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);
+// _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(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(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();
+ }
_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() {
+ 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() {
+ public double[] getCircleParameters()
+ {
return _circleParameters;
- }
- double s(double radius, double xcenter, double ycenter, double x1, double y1, double x2, double y2) {
+ }
+ 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);
@@ -232,27 +355,70 @@
* Reorder TrackerHits from closest to furthest from origin in z.
* @param hits list of TrackerHits
*/
- public void order(List<TrackerHit> hits) {
+ private void order(List<TrackerHit> hitslist)
+ {
Comparator<TrackerHit> z = new Comp();
- Collections.sort(hits, z);
- double[] hit1 = hits.get(0).getPosition();
- double[] hit2 = hits.get(hits.size()-1).getPosition();
+ 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(hits);
+ 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};
+
+ return modInfo;
+ }
}
class Comp implements Comparator<TrackerHit>
{
- public int compare(TrackerHit hit1, TrackerHit hit2) {
+ 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)
+ else if (z1 > z2)
return 1;
else
return 0;
