lcd/hep/lcd/recon/tracking/fitters/sld
diff -u -r1.2 -r1.3
--- WeightMatrix.java 26 Mar 2002 18:28:26 -0000 1.2
+++ WeightMatrix.java 10 Aug 2005 20:31:14 -0000 1.3
@@ -31,6 +31,11 @@
SIGij = SUMk[(Li-Lk+dLk/sqrt(3))*THETk*(Lj-Lk+dLk/sqrt(3)*THETk]
(I follow here prescription from MOURS Benoit, who created fitting routines for SLD)
+
+ However, 1/sqrt(3) represent displacement in the layer filled with scattering material
+ uniformely. In our case, we have very thin dense scater layer followd by much larger
+ gap of air. If we ignore scatter in the air (as we do), when no 1/sqrt(3) coeff. should
+ be used, displacement will be THETk * (Lj-Lk) for last layer also.
During track fitting we will calculate weight matrix for each track only once, assuming that it
does not change considerably due to track parameters change during the fit.
@@ -75,6 +80,7 @@
private double[] lradii;
private boolean[] active;
private int npnts=0;
+ private int nlrstot = 0;
private boolean useBP = false;
private boolean useTPC = false;
private boolean useCCD = false;
@@ -125,7 +131,7 @@
private void doWeightMatrix()
{
- df.setMaximumFractionDigits(5);
+ df.setMaximumFractionDigits(7);
if (dbg) System.out.println("Calculating weight matrix");
PropertySet tpc = (PropertySet) d.get("MultiLayer_TRACK_Barrel");
useTPC = false;
@@ -152,12 +158,12 @@
trpmat.putParams(tpar);
double[] lengths = trpmat.getLengths();
lradii = trpmat.getRadii();
+ nlrstot = lradii.length;
active = trpmat.getActive();
npnts=0;
for(int i=0; i<trpmat.getNLayers(); i++)
if(active[i]) npnts++;
- useBP = true;
itpcr = 0;
PropertySet vtx = (PropertySet) d.get("MultiLayer_VERTX_Barrel");
if(vtx != null)
@@ -205,6 +211,7 @@
* REMEMBER - all resolutions, residuals and so on are in cm .
*/
double pim = 0.1396; // assume pion mass for calculating beta);
+// double pim = 1.; // assume proton mass for calculating beta);
double oosq3 = 0.57735;
double merr[][] = new double[2*npnts][2*npnts];
for(int i=0; i<2*npnts; i++)
@@ -224,7 +231,9 @@
double coef = 0.0136/(beta*ptot);
double thetk[] = new double[npnts];
df.setMaximumFractionDigits(5);
- if(useBP) itpcr++;
+// if(useBP) itpcr++;
+ itpcr++;
+ if(!useBP) active[0]=false;
for(int i=0; i<npnts; i++)
{
swmr.swim(lradii[i]);
@@ -234,10 +243,13 @@
if(sinxs > 0.09) xanc[i] = 1./(sinxs * sinxs);
if(xanc[i] < 1.)
System.out.println("Cros. angle corr. for point "+i+" : "+xanc[i]);
- if(xanc[i] > 100.) xanc[i] = 100.;
- double eln = lradlen[i];
+ if(debug) System.out.println("Cros. angle corr for pnt "+i+" : "+df.format(xanc[i]));
+ if(xanc[i] > 100.) xanc[i] = 100.;
+ double eln = lradlen[i]+0.000033*(pathlen[i+1]-pathlen[i]);
+ if(debug) System.out.println("radiation length of layer "+i+" : "+df.format(eln));
double lcor = 1. + 0.038 * Math.log(eln);
- thetk[i] = coef * lcor * Math.sqrt(eln);
+ thetk[i] = 1.2 * coef * lcor * Math.sqrt(eln);
+ if(debug) System.out.println("Scatter angle in layer "+i+" : "+df.format(thetk[i]));
}
double grp = (double) ngrp;
@@ -256,29 +268,41 @@
double msersq = 0.;
if(j > 0)
{
- for(int k = 0; k<j; k++)
+ for(int k = 0; k<j+1; k++)
{
- double dLkt = (oosq3+0.5) * Math.abs((pathlen[k+1]-pathlen[k]));
- if(k==i) dLkt=0.;
+ double dLkt = 0.;
+ if(k < nlrstot-1) dLkt= Math.abs((lradii[k+1]-lradii[k]));
+ if(k==nlrstot-1) dLkt = Math.abs((lradii[k]-lradii[k-1]));
+// if(k==i) dLkt=0.;
if(active[k])
{
double delpi = 0.;
double delpj = 0.;
- if(i > k+1) delpi = pathlen[i]-pathlen[k+1];
- if(j > k+1) delpj = pathlen[j]-pathlen[k+1];
- msersq+= thetk[k]*thetk[k]*(delpi+dLkt)*(delpj+dLkt);
+ double ersqad = 0.;
+// if(i > k+1) delpi = pathlen[i]-pathlen[k+1];
+// if(j > k+1) delpj = pathlen[j]-pathlen[k+1];
+ if((k < i+1) && (k < nlrstot-1)) delpi = lradii[i]-lradii[k];
+ if((k < j+1) && (k < nlrstot-1)) delpj = lradii[j]-lradii[k];
+ ersqad = thetk[k]*thetk[k]*(delpi)*(delpj);
+ if(((k > i-1) || (k > j-1)) && (i != j)) ersqad = 0.;
+ msersq+= ersqad;
}
}
}
- msersq = msersq * xanc[j];
+// msersq = msersq * xanc[j];
if(j==i) merr[i][j] = me*me + msersq ;
- if(i != j) merr[i][j] = msersq ;
- if(i != j) merr[j][i] = merr[i][j];
+// if(i != j) merr[i][j] = msersq ;
+ if(i != j)
+ {
+ merr[i][j] = 0.5 * msersq;
+ merr[j][i] = 0.5 * msersq;
+ }
+// if(i != j) merr[j][i] = merr[i][j];
}
}
for(int i = npnts; i < 2*npnts; i++)
{
- double me = 0.;
+ double me = 9999.;
for(int j = 0; j < i+1; j++)
{
int ie = i-npnts;
@@ -293,46 +317,77 @@
double msersq = 0.;
if(je > 1)
{
- for(int k = 1; k<je; k++)
+ for(int k = 1; k<je+1; k++)
{
- double dLkt = (oosq3+0.5) * (pathlen[k+1]-pathlen[k]);
- if(k==ie) dLkt = 0.;
+ double dLkt = 0.;
+ if(k<nlrstot-1) dLkt= (pathlen[k+1]-pathlen[k]);
+// if(k==ie) dLkt = 0.;
double delpi = 0.;
double delpj = 0.;
- if(ie > k+1) delpi = pathlen[ie]-pathlen[k+1];
- if(je > k+1) delpj = pathlen[je]-pathlen[k+1];
- if(active[k])
- msersq+= thetk[k]*thetk[k]*(delpi+dLkt)*(delpj+dLkt);
+ double ersqad = 0.;
+ if(ie > k+1) delpi = pathlen[ie]-pathlen[k+1];
+ if(je > k+1) delpj = pathlen[je]-pathlen[k+1];
+ ersqad = thetk[k]*thetk[k]*(delpi)*(delpj);
+ if(((k > ie-1) || (k > je-1)) && (i != j)) ersqad = 0.;
+ if(active[k])
+ msersq+= ersqad;
}
}
msersq*=oneovcls;
if(i==j) merr[i][j] = me*me + msersq;
- if(i != j) merr[i][j] = msersq ;
- if(i != j) merr[j][i] = merr[i][j];
+// if(i != j) merr[i][j] = msersq ;
+// if(i != j) merr[j][i] = merr[i][j];
+ if(i != j)
+ {
+ merr[i][j]=0.5 * msersq;
+ merr[j][i]=0.5 * msersq;
+ }
}
}
if(useBP)
{
- merr[0][0] = 999.;
- merr[npnts][npnts] = 999.;
+ merr[0][0] = 0.000001;
+ merr[npnts][npnts] = 0.000001;
+ }
+ else
+ {
+ merr[0][0] = 99999.;
+ merr[npnts][npnts] = 99999.;
}
- /*
+ if(debug)
+ {
for(int i=0; i<npnts; i++)
{
- double errinxy = merr[i][i] * 10000.;
- double errinz = merr[i+npnts][i+npnts] * 10000.;
+ double errinxy = Math.sqrt(merr[i][i]);
+ double errinz = Math.sqrt(merr[i+npnts][i+npnts]);
System.out.println("layer "+i+" RPHI err = "+df.format(errinxy)+" Z err = "+df.format(errinz));
}
- */
+ }
Matrix ers = new Matrix(merr);
wm = ers.inverse();
- if(debug) System.out.println("Weight matrix calculated");
- for(int i=0; i<npnts; i++)
+/* for(int i=0; i<2*npnts; i++)
+ {
+ for(int j=0; j<2*npnts; j++)
+ {
+ if(i==j)
+ {
+ if(merr[i][j] > 0.0000000001) merr[i][j]=1./merr[i][j];
+ else merr[i][j]=0.;
+ }
+ else merr[i][j]=0.;
+ }
+ }
+ wm = new Matrix(merr); */
+ if(debug)
{
- double wtxy = wm.get(i,i);
- double wtz = wm.get(i+npnts,i+npnts);
-// System.out.println("layer "+i+" RPHI wt = "+df.format(wtxy)+" Z wt = "+df.format(wtz));
- }
+ System.out.println("Weight matrix calculated");
+ for(int i=0; i<npnts; i++)
+ {
+ double wtxy = wm.get(i,i);
+ double wtz = wm.get(i+npnts,i+npnts);
+ System.out.println("layer "+i+" RPHI wt = "+df.format(wtxy)+" Z wt = "+df.format(wtz));
+ }
+ }
}
