lcd/hep/lcd/recon/tracking/fitters/sld
diff -u -r1.4 -r1.5
--- SLDFitTrack.java 21 Mar 2005 17:31:05 -0000 1.4
+++ SLDFitTrack.java 10 Aug 2005 20:30:23 -0000 1.5
@@ -24,6 +24,7 @@
public boolean useVXD = true;
public boolean useTPC = true;
+ public int lastTrToTr = 3;
private final int MaxLayerHits = 1000;
private final int MaxLayersNmbr = 160;
public double GoodChi2 = 10.;
@@ -83,12 +84,12 @@
private DecimalFormat df = new DecimalFormat();
private boolean printRes = false;
-// Histogram h1;
-// Histogram h2;
-// Histogram h3;
-// Histogram h4;
-// Histogram h5;
-// Histogram h6;
+/* Histogram h1;
+ Histogram h2;
+ Histogram h3;
+ Histogram h4;
+ Histogram h5;
+ Histogram h6; */
double[] initpar = new double[6];
double[] finalpar = new double[6];
private boolean dohist = true;
@@ -98,8 +99,7 @@
public SLDFitTrack()
{
- /*
- if(dohist)
+/* if(dohist)
{
String folder = "/SLDFitTrack";
HistogramFolder.setDefaultFolder(folder);
@@ -109,8 +109,7 @@
h4 = new Histogram("change in z0");
h5 = new Histogram("change in tanlam");
h6 = new Histogram("change in chi2");
- }
- */
+ } */
}
public void fit(int it,VertexDetector v,Tracker t)
@@ -158,6 +157,10 @@
tchg = (tkpar.omega > 0.)? 1 : -1;
Field_Strength field = (Field_Strength) det.get("Field_Strength");
BField = field.getField();
+ df.setMaximumFractionDigits(5);
+ if(debug) System.out.println("Initial trk par - d0: "+df.format(tkpar.d0)+" phi0: "+
+ df.format(tkpar.phi0)+" omega: "+df.format(tkpar.omega)+" z0: "+df.format(tkpar.z0)+
+ " tanlam: "+df.format(tkpar.s));
double r_cyl = 100.;
double z_pl = 500.;
swmr.setBField(BField);
@@ -231,6 +234,12 @@
else return veryLarge;
}
+ public double getUnfitChi2()
+ {
+ if(!badTrack) return initpar[5];
+ else return veryLarge;
+ }
+
public double getNDF()
{
return NDoF;
@@ -291,7 +300,7 @@
| Routine takes care about grouping layers for large detector, so that about 12 measurement
| points will be generated for large tracker.
| It also gives estimate for layers active radii. This radii will be corrected during hit
- | list filling, when real hitt coordinate will be used. For now we will be using geometry file
+ | list filling, when real hit coordinate will be used. For now we will be using geometry file
| inner layer radius + 1/2 of layer thickness as an estimation for hit radius.
|_______________________________________________________________________________________________
*/
@@ -443,7 +452,7 @@
private void getTrkInters(TrkParams tp)
{
- double[] fr0 = new double[3];
+/* double[] fr0 = new double[3];
double[] fp0 = new double[3];
@@ -455,14 +464,19 @@
fp0[1] = fpt * Math.sin(tp.phi0);
fp0[2] = fpt * tp.s;
int chrg = (tp.omega > 0) ? 1 : -1 ;
- swmr.swim(fp0,fr0,chrg);
+ swmr.swim(fp0,fr0,chrg); */
+ swmr.swim(tp);
for(int i=1; i<nvlrs+ntlrs+1; i++)
{
double rl = lrsr[i];
- swmr.swim(rl);
+ double zm = 500.;
+ swmr.swim(rl,zm);
double[] xx = swmr.getIntersect();
trkz[i] = xx[2];
trkphi[i] = Math.atan2(xx[1],xx[0]);
+ df.setMaximumFractionDigits(5);
+ if(debug) System.out.println("Lr "+i+" R: "+df.format(rl)+" trk x: "+df.format(xx[0])+
+ " trk y: "+df.format(xx[1]));
}
}
@@ -525,8 +539,8 @@
double[] bestch2 = new double[MaxLayersNmbr];
- double[] fr0 = new double[3];
- double[] fp0 = new double[3];
+/* double[] fr0 = new double[3];
+ double[] fp0 = new double[3]; */
double[] bestmch= new double[MaxLayersNmbr];
double fpt = BField * 0.003 /Math.abs(tkpar.omega);
@@ -546,18 +560,19 @@
// now swim track through correct radii
- fr0[0] = tkpar.d0 * Math.cos(tkpar.phi0+halfpi);
+/* fr0[0] = tkpar.d0 * Math.cos(tkpar.phi0+halfpi);
fr0[1] = tkpar.d0 * Math.sin(tkpar.phi0+halfpi);
fr0[2] = tkpar.z0;
fp0[0] = fpt * Math.cos(tkpar.phi0);
fp0[1] = fpt * Math.sin(tkpar.phi0);
fp0[2] = fpt * tkpar.s;
- int chrg = (tkpar.omega > 0) ? 1 : -1 ;
- swmr.swim(fp0,fr0,chrg);
+ int chrg = (tkpar.omega > 0) ? 1 : -1 ; */
+ swmr.swim(tkpar);
for(int i=1; i<nvlrs+ntlrs*ngrp+1; i++)
{
double rl = slrsr[i];
- swmr.swim(rl);
+ double zm = 500.;
+ swmr.swim(rl,zm);
double[] xx = swmr.getIntersect();
sltz[i] = xx[2];
sltphi[i] = Math.atan2(xx[1],xx[0]);
@@ -744,8 +759,11 @@
if(nlr>127) nlr = nlr-128;
someHit hit = new someHit();
hit.putPoint(point);
- layerhits[nlr+1][0]=hit;
- nlhits[nlr+1]=1;
+ if(!vh.isEndcap())
+ {
+ layerhits[nlr+1][0]=hit;
+ nlhits[nlr+1]=1;
+ }
}
int nth = tpc.getNPoints(it);
for(int ip=1; ip<nth+1; ip++)
@@ -755,7 +773,7 @@
TrackerHit th = tpc.getTrkHit(it,ip);
int nlr = th.getLayer();
nlr+=6;
- if(hra > maxvr)
+ if((hra > maxvr) && !th.isEndcap())
{
someHit hit = new someHit();
hit.putPoint(point);
@@ -805,16 +823,27 @@
if(!active[i]){ res[i]=0; res[i+nps] = 0; }
if((i != 0) && active[i])
{
- res[i] = (trkphi[i]-hitsphi[i]) * lrsr[i];
+ double dph = trkphi[i]-hitsphi[i];
+ if(Math.abs(dph) > Math.PI)
+ {
+ if(dph > 0.) dph-= (2.* Math.PI);
+ if(dph < 0.) dph+= (2. * Math.PI);
+ }
+ res[i] = dph * lrsr[i];
res[i+nps] = trkz[i] - hitsz[i];
}
}
if(printRes)
{
- // printRes = false;
- for(int i=10; i< nps; i++)
- System.out.println("point "+i+" trkphi= "+df.format(trkphi[i])+
- " hitphi= "+df.format(hitsphi[i])+" res= "+df.format(res[i]));
+ printRes = false;
+ for(int i=1; i< nps; i++)
+ {
+ double hitx = lrsr[i]*Math.cos(hitsphi[i]);
+ double hity = lrsr[i]*Math.sin(hitsphi[i]);
+ System.out.println("pnt "+i+" trkphi= "+df.format(trkphi[i])+
+ " hitphi= "+df.format(hitsphi[i])+" res= "+df.format(res[i])+
+ "res Z: "+df.format(res[i+nps])+" htx: "+df.format(hitx)+" hty: "+df.format(hity));
+ }
}
}
@@ -895,18 +924,23 @@
// wp.print(8,6);
Matrix y = new Matrix(Ym,nrs);
// y.print(8,6);
- Matrix dp = wp.solve(y);
- if(debug) System.out.println("Parameter increments: ");
- if(debug) dp.print(8,6);
+ LUDecomposition lud = new LUDecomposition(wp);
+ if(lud.isNonsingular())
+ {
+ Matrix dp = wp.solve(y);
+ if(debug) System.out.println("Parameter increments: ");
+ if(debug) dp.print(8,6);
// System.out.println("Error matrix: ");
// erma.print(9,7);
- if(debug && success)
- {
+ if(debug && success)
+ {
for(int i=0; i<5; i++)
System.out.println("Error in parameter "+i+" is "+Math.sqrt(erma.get(i,i)));
- }
- for(int i=0; i<5; i++)
+ }
+ for(int i=0; i<5; i++)
dpar[i] = dp.get(i,0);
+ }
+ else badTrack = true;
}
private void calculateChi2()
@@ -940,14 +974,15 @@
{
success = false;
badTrack = false;
- for(int nit = 0; nit<2; nit++)
+ for(int nit = 0; nit<4; nit++)
{
- if(nit == 0) delpar = 0.000001;
- else delpar = 0.0000001;
+ if(debug) printRes = true;
+ if(nit == 0) delpar = 0.00001;
+ else delpar = 0.000001;
getTrkInters(tkpar);
getResiduals(residuals);
calculateChi2();
- if(nit == 0) initpar[5]=chi2;
+ if(nit == 0) initpar[5]=chi2;
if(chi2 > veryLarge-1.)
{
if(debug)
@@ -957,10 +992,13 @@
if(debug) System.out.println("itteratin "+nit+" chi2= "+chi2);
getDerivatives();
getSolution();
+ if(badTrack) { chi2 = veryLarge; return; }
+ boolean nochange = true;
for(int i=0; i<5; i++)
{
double maxdev = 0.1;
- if((i==0) || (i==3)) maxdev = 1.0;
+ if((i==0) || (i==3)) maxdev = 1.0;
+ if(Math.abs(dpar[i]) > delpar) nochange = false;
if(Math.abs(dpar[i]) > maxdev)
{
badTrack = true;
@@ -975,6 +1013,7 @@
tkpar.omega+=dpar[2];
tkpar.z0+=dpar[3];
tkpar.s+=dpar[4];
+ if(nochange) break;
}
getTrkInters(tkpar);
getResiduals(residuals);
@@ -985,8 +1024,7 @@
finalpar[2]=tkpar.omega;
finalpar[3]=tkpar.z0;
finalpar[4]=tkpar.s;
- /*
- if(dohist)
+ /* if(dohist)
{
h1.fill(finalpar[0]-initpar[0]);
h2.fill(finalpar[1]-initpar[1]);
@@ -994,8 +1032,7 @@
h4.fill(finalpar[3]-initpar[3]);
h5.fill(finalpar[4]-initpar[4]);
if(NDoF > 1.) h6.fill((finalpar[5]-initpar[5])/NDoF);
- }
- */
+ } */
if(debug)
{
System.out.println("final chi2= "+df.format(chi2));
