4 modified files
lcsim/src/org/lcsim/contrib/NickSinev/event/util
diff -u -r1.1 -r1.2
--- SmearTrackerHits.java 2 Aug 2006 01:29:46 -0000 1.1
+++ SmearTrackerHits.java 9 Aug 2006 20:28:35 -0000 1.2
@@ -23,10 +23,18 @@
/**
* @author sinev U of Oregon; SLAC x2970; [log in to unmask]
- * @version $Id: SmearTrackerHits.java,v 1.1 2006/08/02 01:29:46 sinev Exp $
+ * @version $Id: SmearTrackerHits.java,v 1.2 2006/08/09 20:28:35 sinev Exp $
*/
public class SmearTrackerHits extends Driver
{
+ private double bvtxClSpanR = 0.1;
+ private double bvtxClSpanCS = 0.1;
+ private double evtxClSpanZ = 0.1;
+ private double evtxClSpanDS = 0.1;
+ private double btrkClSpanR = 0.5;
+ private double btrkClSpanCS = 0.5;
+ private double etrkClSpanZ = 0.5;
+ private double etrkClSpanDS = 0.5;
private Vector svbhits = null ;
private Vector svehits = null;
private Vector trbhits = null;
@@ -73,6 +81,14 @@
public void setTrkECRPResoloution(double res) { trkecRPhism = res;}
public void setTrkEcRResoloution(double res) { trkecRsm = res;}
public void setDoSmear(boolean dosm) { doSmear=dosm; }
+ public void setVtxBarrHitClustRSpan(double rsp) { bvtxClSpanR = rsp; }
+ public void setVtxBarrHitClustTransSpan(double tsp) { bvtxClSpanCS = tsp; }
+ public void setVtxECHitClustRZpan(double zsp) { evtxClSpanZ = zsp; }
+ public void setVtxECHitClustTransSpan(double tsp) { evtxClSpanDS = tsp; }
+ public void setTrkBarrHitClustRSpan(double rsp) { btrkClSpanR = rsp; }
+ public void setTrkBarrHitClustTransSpan(double tsp) { btrkClSpanCS = tsp; }
+ public void setTrkECHitClustRZpan(double zsp) { etrkClSpanZ = zsp; }
+ public void setTrkECHitClustTransSpan(double tsp) { etrkClSpanDS = tsp; }
public void process(EventHeader event)
{
@@ -110,12 +126,15 @@
double rcl = Math.sqrt(clcc[0]*clcc[0]+clcc[1]*clcc[1]);
double rht = Math.sqrt(htc[0]*htc[0]+htc[1]*htc[1]);
double dr = rcl-rht;
- if(Math.abs(dr)<0.1)
+ if(Math.abs(dr)<bvtxClSpanR)
{
- double dd = Math.sqrt((clcc[0]-htc[0])*(clcc[0]-htc[0])+
- (clcc[1]-htc[1])*(clcc[1]-htc[1])+
- (clcc[2]-htc[2])*(clcc[2]-htc[2]));
- if(dd<0.1)
+ double htaX = htc[0]*rcl/rht;
+ double htaY = htc[1]*rcl/rht;
+ double htaZ = htc[2]*rcl/rht;
+ double dd = Math.sqrt((clcc[0]-htaX)*(clcc[0]-htaX)+
+ (clcc[1]-htaY)*(clcc[1]-htaY)+
+ (clcc[2]-htaZ)*(clcc[2]-htaZ));
+ if(dd<bvtxClSpanCS)
{
newcl=false;
cof.addHit(th);
@@ -185,12 +204,15 @@
double zcl = clcc[2];
double zht = htc[2];
double dz = zcl-zht;
- if(Math.abs(dz)<0.1)
+ double azcl = Math.abs(zcl);
+ double azht = Math.abs(zht);
+ if(Math.abs(dz)<evtxClSpanZ)
{
- double dd = Math.sqrt((clcc[0]-htc[0])*(clcc[0]-htc[0])+
- (clcc[1]-htc[1])*(clcc[1]-htc[1])+
- (clcc[2]-htc[2])*(clcc[2]-htc[2]));
- if(dd<0.1)
+ double htaX = htc[0]*azcl/azht;
+ double htaY = htc[1]*azcl/azht;
+ double dd = Math.sqrt((clcc[0]-htaX)*(clcc[0]-htaX)+
+ (clcc[1]-htaY)*(clcc[1]-htaY));
+ if(dd<evtxClSpanDS)
{
newcl=false;
cof.addHit(th);
@@ -260,12 +282,14 @@
double rcl = Math.sqrt(clcc[0]*clcc[0]+clcc[1]*clcc[1]);
double rht = Math.sqrt(htc[0]*htc[0]+htc[1]*htc[1]);
double dr = rcl-rht;
- if(Math.abs(dr)<0.3)
+ if(Math.abs(dr)<btrkClSpanR)
{
- double dd = Math.sqrt((clcc[0]-htc[0])*(clcc[0]-htc[0])+
- (clcc[1]-htc[1])*(clcc[1]-htc[1])+
- (clcc[2]-htc[2])*(clcc[2]-htc[2]));
- if(dd<0.3)
+ double htaX = htc[0]*rcl/rht;
+ double htaY = htc[1]*rcl/rht;
+ double htaZ = htc[2]*rcl/rht;
+ double dd = Math.sqrt((clcc[0]-htaX)*(clcc[0]-htaX)+
+ (clcc[1]-htaY)*(clcc[1]-htaY)+(clcc[2]-htaZ)*(clcc[2]-htaZ));
+ if(dd<btrkClSpanCS)
{
newcl=false;
cof.addHit(th);
@@ -323,12 +347,15 @@
double zcl = clcc[2];
double zht = htc[2];
double dz = zcl-zht;
- if(Math.abs(dz)<0.3)
+ double azcl = Math.abs(zcl);
+ double azht = Math.abs(zht);
+ if(Math.abs(dz)<etrkClSpanZ)
{
- double dd = Math.sqrt((clcc[0]-htc[0])*(clcc[0]-htc[0])+
- (clcc[1]-htc[1])*(clcc[1]-htc[1])+
- (clcc[2]-htc[2])*(clcc[2]-htc[2]));
- if(dd<0.3)
+ double htaX = htc[0]*azcl/azht;
+ double htaY = htc[1]*azcl/azht;
+ double dd = Math.sqrt((clcc[0]-htaX)*(clcc[0]-htaX)+
+ (clcc[1]-htaY)*(clcc[1]-htaY));
+ if(dd<etrkClSpanDS)
{
newcl=false;
cof.addHit(th);
lcsim/src/org/lcsim/contrib/NickSinev/tracking/util
diff -u -r1.1 -r1.2
--- WMTrackPropagator.java 2 Aug 2006 01:29:46 -0000 1.1
+++ WMTrackPropagator.java 9 Aug 2006 20:28:35 -0000 1.2
@@ -27,7 +27,7 @@
/**
* @author Nick Sinev U of Oregon, [log in to unmask] SLAC x2970
- * @version $Id: WMTrackPropagator.java,v 1.1 2006/08/02 01:29:46 sinev Exp $
+ * @version $Id: WMTrackPropagator.java,v 1.2 2006/08/09 20:28:35 sinev Exp $
*/
public class WMTrackPropagator
@@ -116,6 +116,7 @@
MaterialManager mama=null;
double airSpecX0 = 0.;
double airDens = 0.;
+ double totelos = 0.;
private DecimalFormat df = new DecimalFormat();
private double mPi = 0.139;
@@ -177,10 +178,12 @@
public Matrix getProjectedErrorMatrix() { return projErrorMatrix; }
public int getNActiveLrs() { return naclt; }
+
+ public double getTotalEnergyLoss() { return totelos; }
public void init(Detector detector)
{
- df.setMaximumFractionDigits(1);
+ df.setMaximumFractionDigits(5);
if((detector != null) &&(!initialized))
{
det = detector;
@@ -584,6 +587,7 @@
void propagate( double[] tpar)
{
+ totelos = 0.;
for(int i=0; i<5; i++) trpar[i]=tpar[i];
double ttl = trpar[4];
trackPt = Bfield*Constants.fieldConversion/Math.abs(trpar[2]);
@@ -601,6 +605,8 @@
if(!Double.isNaN(len)) maxlen = Math.abs(len);
else maxlen = swm.getDistanceToZ(ttls*maxz);
if(maxlen < 0.) maxlen = -maxlen;
+ if(maxlen > 6.*Math.PI*trkCurvRad) maxlen = 6.*Math.PI*trkCurvRad;
+ if(maxlen > 3.* maxz) maxlen = 3.* maxz;
// System.out.println("Max trk length: "+df.format(maxlen));
hitbytr = new boolean[nSubLrs];
lentoh = new double[nSubLrs];
@@ -632,7 +638,8 @@
{
boolean hit = true;
len = swm.getDistanceToZ(ttls*minZ[sl]);
- if((len > Math.PI*trkCurvRad) && (trkCurvRad > 0.5*maxR[sl])) hit=false;
+ if((len > Math.abs(ttl)*Math.PI*trkCurvRad) && (trkCurvRad > 0.5*maxR[sl])) hit=false;
+ if(6.*Math.PI*trkCurvRad*Math.abs(ttl) < Math.abs(minZ[sl])) hit=false;
if(Math.abs(len) > maxlen) hit=false;
if(hit)
{
@@ -647,11 +654,11 @@
entryCyl[sl]=false;
nhlrs++;
}
- else
+ else if(trkCurvRad > 0.5*minR[sl])
{
len = swm.getDistanceToRadius(minR[sl]);
pnt = swm.getPointAtDistance(len);
- if((Math.abs(pnt.z())<maxZ[sl])&&(Math.abs(pnt.z())>minZ[sl]))
+ if((Math.abs(pnt.z())<maxZ[sl])&&(Math.abs(pnt.z())>minZ[sl])&&(len<maxlen))
{
// System.out.println("Hit inner radius of disk at Z= "+pnt.z()+" and radius "+minR[sl]);
hitbytr[sl]=true;
@@ -792,6 +799,7 @@
double eloss = 0.;
if(doelos) eloss = 0.0017*avDens[i]*(oslol[j]-oslel[j]); //as we do not have de/dx for particular mat,
// we estimate it as 1.7 MeV/g/cm**2
+ totelos+=eloss;
double mpath = oslol[j]-oslel[j];
if(rllen<0.) System.out.println("Negative path! Output length "+oslol[j]+" while inp "+oslel[j]);
if(oslel[j] < 1.) System.out.println("Zero track length to layer "+j);
@@ -799,7 +807,7 @@
trktotrlp+=rllen;
double outscang = scc*Math.sqrt(accx)*(1.+0.038*Math.log(accx));
double clsca = Math.sqrt(outscang*outscang-acscasq);
- clsca = scc * Math.sqrt(rllen)*(1.+0.038*Math.log(rllen));
+// clsca = scc * Math.sqrt(rllen)*(1.+0.038*Math.log(rllen)); //+++++++++++++++++++++++++
acscasq+=clsca*clsca;
// scang[i]=scc*Math.sqrt(rllen)*(1.+0.038*Math.log(rllen));
indinhl[nlrwa]=j;
@@ -885,7 +893,7 @@
naclt=0;
for(int i=0; i<nlrwa; i++)
{
- int k = indinhl[i]; //k==rank of layer (order from start of tr) from all hit lrs
+ int k = indinhl[i]; //k==rank of layer (order from start of tr) in all hit lrs
if(k!=-1)
{
activeli[k]=-1;
@@ -994,7 +1002,8 @@
double accer=0.;
for(int k=0; k<posi; k++)
{
- double armlen = chordl[indinhl[posi]]-lentoesp[k];
+// double armlen = chordl[indinhl[posi]]-lentoesp[k];
+ double armlen = lentoesp[posi]-lentoesp[k];
accer+=armlen*armlen*scinesp[k]*scinesp[k];
if(accer > 10000.) countab = false;
}
@@ -1033,7 +1042,7 @@
double levj = lenj-lentoesp[k];
double dRi = hitR[posi]-hitR[k];
double dRj = hitR[posj]-hitR[k];
- double dPt = scinesp[k] * trkfulmom*trsinl;
+ double dPt = scinesp[k] * trkfulmom * trsinl;
double dRc = trkCurvRad*dPt/trackPt;
if(dopri)
{
@@ -1047,7 +1056,10 @@
double dispfrci = Math.abs(0.5 * dRi*dRi * dRc / (trkCurvRad * trkCurvRad));
double dispfrcj = Math.abs(0.5 * dRj*dRj * dRc / (trkCurvRad * trkCurvRad));
// mel+= scinesp[k]*scinesp[k]*levj*levi;
- mel+= (scinesp[k]*scinesp[k]*levj*levi + dispfrci*dispfrcj);
+ mel+= (scinesp[k]*scinesp[k]*levj*levi + dispfrci*dispfrcj);
+ if(scinesp[k]*scinesp[k]*levj*levi > 100000.)
+ System.out.println("from layer "+k+" to lrs "+posi+" "+posj+" sc. ang "+df.format(scinesp[k])+
+ " lever i "+df.format(levi)+" lever j "+df.format(levj));
}
if(j==i)
{
@@ -1057,8 +1069,8 @@
{
if(inccor)
{
- errorMatrix.set(j,i, 0.5*mel);
- errorMatrix.set(i,j, 0.5*mel);
+ errorMatrix.set(j,i, mel); //+++++++++++++++++++++++++++++++++++++++++++
+ errorMatrix.set(i,j, mel); //+++++++++++++++++++++++++++++++++++++++++++
}
if(!inccor)
{
@@ -1103,18 +1115,22 @@
{
int ii = actLrInd[i];
int lrini = subLrInd[ii];
- double norm = 1./Math.sqrt(actLrPnts[i].x()*actLrPnts[i].x()+actLrPnts[i].y()*actLrPnts[i].y());
- Hep3Vector dir1i = new BasicHep3Vector(norm*actLrPnts[i].y(),-norm*actLrPnts[i].x(),0.);
+// double norm = 1./Math.sqrt(actLrPnts[i].x()*actLrPnts[i].x()+actLrPnts[i].y()*actLrPnts[i].y());
+ Hep3Vector dir1inn = new BasicHep3Vector(actLrPnts[i].y(),-actLrPnts[i].x(),0.);
+ Hep3Vector dir1i = VecOp.unit(dir1inn);
Hep3Vector dir2i = new BasicHep3Vector(0.,0.,1.);
// Hep3Vector trkchdir = new BasicHep3Vector(actLrPnts[i].x(), actLrPnts[i].y(), actLrPnts[i].z());
Hep3Vector trkchdir = trkMomAt[i];
+ Hep3Vector trkdir = VecOp.unit(trkchdir);
// if(i>0) trkchdir = VecOp.sub(actLrPnts[i], actLrPnts[i-1]);
sigi = 1.;
if(!isCylinder[lrini])
{
- dir2i = new BasicHep3Vector(norm*actLrPnts[i].x(),norm*actLrPnts[i].y(),0.);
+ Hep3Vector dir2inn = new BasicHep3Vector(actLrPnts[i].x(),actLrPnts[i].y(),0.);
+ dir2i = VecOp.unit(dir2inn);
sigi = -1.;
}
+ /*
Hep3Vector dirsc1 = VecOp.cross(dir2i,trkchdir);
Hep3Vector dirsc2 = VecOp.cross(dirsc1,trkchdir);
double cosp1i = Math.abs(VecOp.dot(dir1i,dirsc1)/
@@ -1126,27 +1142,22 @@
if(cosp1i < 0.1) cosp1i = 0.1;
if(cosp2i < 0.1) cosp2i = 0.1;
double cor1i = 1./cosp1i;
- double cor2i = 1./cosp2i;
-/* if(isCylinder[lrini])
- {
- if(Math.abs(cor2i-Math.abs(trkfulmom/trackPt))>0.01)
- {
- System.out.println("Lr: "+i+" ii "+ii+" lrini "+lrini+" cor2i= "+cor2i+" while P/Pt= "+trkfulmom/trackPt);
- System.out.println("dir2i: "+dir2i.x()+" "+dir2i.y()+" "+dir2i.z());
- System.out.println("dirsc2: "+dirsc2.x()+" "+dirsc2.y()+" "+dirsc2.z());
- System.out.println("trkchdir "+trkchdir.x()+" "+trkchdir.y()+" "+trkchdir.z());
- System.out.println("trpar: "+trpar[0]+" "+trpar[1]+" "+trpar[2]+" "+trpar[3]+" "+trpar[4]);
- System.out.println("Total number of act. lrs traversed by track: "+naclt);
- System.out.println("Total number of sublrs hit by tr: "+nhlrs);
- System.out.println("ind in SL rank length to hit");
- for(int jj=0; jj<nhlrs; jj++)
- System.out.println(hlrn[jj]+" "+rank[jj]+" "+tltoh[jj]);
- }
- }*/
-// if(cor1i > 5.) System.out.println("Too large projection coeff. cor1i: "+df.format(cor1i));
-// if(cor2i > 5.) System.out.println("Too large projection coeff. cor2i: "+df.format(cor2i));
- corco1[i]=cor1i;
- corco2[i]=cor2i;
+ double cor2i = 1./cosp2i; */
+ double cosd1 = VecOp.dot(trkdir, dir1i);
+ if(Math.abs(cosd1) < 0.99) corco1[i]=1./Math.sqrt(1.-cosd1*cosd1);
+ else corco1[i] = 10.;
+ double cosd2 = VecOp.dot(trkdir,dir2i);
+ if(Math.abs(cosd2) < 0.99)
+ {
+ corco2[i]=1./Math.sqrt(1.-cosd2*cosd2);
+ }
+ else
+ {
+ System.out.println(" Track dir: "+df.format(trkdir.x())+" "+df.format(trkdir.y())+" "+
+ df.format(trkdir.z()));
+ System.out.println(" proj dir: "+df.format(dir2i.x())+" "+df.format(dir2i.y())+" "+df.format(dir2i.z()));
+ corco2[i] = 10.;
+ }
if(i>0)
{
if(corco1[i]<corco1[i-1])
@@ -1170,11 +1181,6 @@
{
int jj = actLrInd[j];
int lrinj = subLrInd[jj];
-// norm = 1./Math.sqrt(actLrPnts[j].x()*actLrPnts[j].x()+actLrPnts[j].y()*actLrPnts[j].y());
-// Hep3Vector dir1j = new BasicHep3Vector(norm*actLrPnts[j].y(),-norm*actLrPnts[j].x(),0.);
-// Hep3Vector dir2j = new BasicHep3Vector(0.,0.,1.);
-// trkchdir = new BasicHep3Vector(actLrPnts[j].x(), actLrPnts[j].y(), actLrPnts[j].z());
-// if(j>0) trkchdir = VecOp.sub(actLrPnts[j], actLrPnts[j-1]);
sigj = 1.;
if(!isCylinder[lrinj])
{
@@ -1182,19 +1188,13 @@
}
double cor1j = corco1[j];
double cor2j = corco2[j];
-// double mel = errorMatrix.get(i,j)*cor1i*cor1j;
- double mel = errorMatrix.get(i,j);
-// double mel = errorMatrix.get(i,j);
-// if(i!=j) projErrorMatrix.set(i,j, 0.);
+ double mel = errorMatrix.get(i,j)*cor1i*cor1j; //++++++++++++++++++++++++++++++++
projErrorMatrix.set(i,j, mel);
-// if(j>i) projErrorMatrix.set(i,j, 0.);
mel = errorMatrix.get(i,j)*cor2i*cor2j;
-// mel = errorMatrix.get(i,j);
sigcor= sigi*sigj;
if(sigcor<0) sigcor*= ttls;
if(i==j) sigcor = 1.;
projErrorMatrix.set(i+naclt, j+naclt,sigcor*mel);
-// if(j>i) projErrorMatrix.set(i+naclt, j+naclt,0.);
}
}
boolean badm = false;
@@ -1218,11 +1218,22 @@
}
if(badm)
{
- projErrorMatrix.print(8,6);
+ errorMatrix.print(8,4);
+ }
+/* if(corco2[naclt-1]>9.9)
+ {
System.out.println("coefficients were: ");
for(int i=0; i<naclt; i++)
- System.out.print(" "+df.format(corco1[i]));
- }
+ System.out.print(" "+df.format(corco1[i]));
+ System.out.println(" ");
+ for(int i=0; i<naclt; i++)
+ System.out.print(" "+df.format(corco2[i]));
+ System.out.println(" ");
+ System.out.println(" for track with parameters: ");
+ for(int i=0; i<5; i++)
+ System.out.print(df.format(trpar[i])+" ");
+ System.out.println(" ");
+ } */
}
public double getTrkRadLenAtR(double R)
lcsim/src/org/lcsim/contrib/NickSinev/tracking/wmfitter
diff -u -r1.1 -r1.2
--- SLDFitTrack.java 2 Aug 2006 01:29:47 -0000 1.1
+++ SLDFitTrack.java 9 Aug 2006 20:28:36 -0000 1.2
@@ -55,12 +55,6 @@
private Matrix ermx;
private double[][] wma;
private double[][] ermxa;
- /*
- note: number of residuals and derivatives equals to 2*(npoints+1), where npoints is number
- of points on track. Residuals and derivatives in rphi are the first, in z - next. The
- very first residual is always = 0., because first point corresponds to beam pipe, for which
- there is no measurement, but it is included in weight matrix calculations.
- */
private Detector det;
private double veryLarge = 10000000000.;
private double[] dpar=null;
@@ -86,7 +80,7 @@
private boolean printRes = false;
double[] initpar = new double[6];
double[] finalpar = new double[6];
- WMTrackPropagator trkPr=null;
+// WMTrackPropagator trkPr=null;
private boolean onlyOmegaFit = false;
private boolean oripar = true;
@@ -210,8 +204,8 @@
if(debug) System.out.println("Setting track in Weight Matrix\n");
Wm.setTrack((Track) track);
if(debug) System.out.println("Extracting track propagator\n");
- trkPr = Wm.getPropagator();
- if(trkPr == null) {System.out.println("Propagator is null! can't fit\n"); badTrack=true; }
+// trkPr = Wm.getPropagator();
+// if(trkPr == null) {System.out.println("Propagator is null! can't fit\n"); badTrack=true; }
}
public void assignHits()
@@ -233,7 +227,6 @@
for(int j=0; j<2*naclrs; j++)
{
wma[i][j]=mtr.get(i,j);
-// if(i!=j) wma[i][j]=0.;
if(Double.toString(wma[i][j]).equals("NaN")) badTrack=true;
if((i==j) && (wma[i][j] < 0.)) System.out.println("Negative diagonal element in weight matrix: "+i+","+j);
}
@@ -256,12 +249,9 @@
{
double[] hpo = hit.getPosition();
double rhit = Math.sqrt(hpo[0]*hpo[0]+hpo[1]*hpo[1]);
-// if(rhit < 80.) hasHit[i]=false;
}
if(hasHit[i]) nlrswh++;
}
-// if(debug)
-// System.out.println("Active mask is set for "+nlrswh+" layers");
NDoF = nlrswh * 2 - nFreePar;
if(planarFit) NDoF=nlrswh-nFreePar;
}
@@ -301,42 +291,27 @@
private void getTrkInters()
{
Hep3Vector pos;
- trkPr.setTrack(tkpar);
- Hep3Vector[] pnts = trkPr.getTrkActLrPnts();
+ Wm.setTrackParams(tkpar);
+ Wm.calculateTrkInt();
+ Hep3Vector[] pnts = Wm.getTrkPnts();
int ntpnts = pnts.length;
for(int i=0; i<naclrs; i++)
{
-// if(oripar && !hasHit[i]) System.out.println("Layer "+i+" does not have hit assigned");
- if((i<ntpnts) && (hasHit[i]))
+ if(i<ntpnts)
{
-// hasHit[i]=true;
- isCylinder[i]=trkPr.isActiveCylinder(i);
+ isCylinder[i]=Wm.isCylinder()[i];
if(isCylinder[i])
{
- pos = trkPr.getTrkPntAtR(rH[i]);
+ pos = pnts[i];
phiX[i] = Math.atan2(pos.y(),pos.x());
rX[i] = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
if(Math.abs(rH[i]-rX[i]) > 0.001) System.out.println("wrong inters. radius: "+rX[i]+
"while hit is at radius "+rH[i]);
zX[i] = pos.z();
-/* if((Math.abs(zH[i]-zX[i]) > 100.)&& oripar)
- {
- System.out.println("Cylinder: wrong inters. z: "+zX[i]+
- "while hit is at z "+zH[i]+" for layer "+i+" itter "+nitt);
- System.out.println("params: "+tkpar[0]+" "+tkpar[1]+" "+tkpar[2]+" "+tkpar[3]+" "+tkpar[4]);
- pos = trkPr.getTrkPntAtZ(zX[i]);
- double ratz = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
- System.out.println("Radius at correct z is"+ratz+" while radius at wrong Z "+rX[i]+" and hit r "+rH[i]);
- for(int k=0; k<ntpnts; k++)
- {
- double rtrpnt = Math.sqrt(pnts[k].x()*pnts[k].x()+pnts[k].y()*pnts[k].y());
- System.out.println("Z "+df.format(pnts[k].z())+" r "+df.format(rtrpnt));
- }
- } */
}
if(!isCylinder[i])
{
- pos = trkPr.getTrkPntAtZ(zH[i]);
+ pos = pnts[i];
phiX[i] = Math.atan2(pos.y(),pos.x());
rX[i] = Math.sqrt(pos.x()*pos.x()+pos.y()*pos.y());
zX[i] = pos.z();
@@ -351,7 +326,7 @@
private void getHitPositions()
{
- trkPr.setTrack(tkpar);
+ Wm.setTrackParams(tkpar);
TrackerHit hit;
double[] hpo;
for(int i=0; i<naclrs; i++)
@@ -405,60 +380,6 @@
}
}
- private void findSweetSpot()
- {
- double[] scandp = new double[5];
- scandp[0]=0.001;
- scandp[1]=0.0001;
- scandp[2]=0.04*Math.abs(freePar[2]);
- scandp[3]=0.001;
- scandp[4]=0.0001;
- int keepi = 0;
- int nsteps = 0;
- boolean newmin = true;
- double stchi2 = calculateChi2For(residuals);
- double newchi2 = 10000.;
- while(newmin && (nsteps<100))
- {
- nsteps++;
- newmin = false;
- setTrackPars();
- getTrkInters();
- getResiduals(residuals);
- for(int i=0; i<nFreePar; i++)
- {
- freePar[i]+=scandp[i];
- setTrackPars();
- getTrkInters();
- getResiduals(residuals);
- newchi2=calculateChi2For(residuals);
- if(Math.abs(newchi2)<stchi2)
- {
- stchi2=newchi2;
- keepi = i+1;
- newmin = true;
- }
- freePar[i]-=scandp[i];
- freePar[i]-=scandp[i];
- setTrackPars();
- getTrkInters();
- getResiduals(residuals);
- newchi2=calculateChi2For(residuals);
- if(Math.abs(newchi2)<stchi2)
- {
- stchi2=newchi2;
- keepi = -i-1;
- newmin = true;
- }
- freePar[i]+=scandp[i];
- }
- if(newmin)
- {
- if(keepi > 0) freePar[keepi-1]+=scandp[keepi-1];
- if(keepi < 0) freePar[-keepi-1]-=scandp[-keepi-1];
- }
- }
- }
private void getDerivatives()
{
@@ -530,10 +451,11 @@
{
for(int m=0; m<nFreePar; m++)
{
+ ym=0.;
for(int n=0; n<nFreePar; n++)
{
double me = 0.;
- ym = 0.;
+// ym = 0.;
for(int i=0; i<nps; i++)
{
int lin = 0;
@@ -662,7 +584,7 @@
if(i < naclrs) lin = i;
if(i > naclrs-1) lin = i - naclrs;
if(i > (wma[0].length - 1)) break;
- for(int j=0; j<i+1; j++)
+ for(int j=0; j<nps; j++)
{
int ljn = 0;
if(j<naclrs) ljn = j;
lcsim/src/org/lcsim/contrib/NickSinev/tracking/wmfitter
diff -u -r1.1 -r1.2
--- WeightMatrix.java 2 Aug 2006 01:29:47 -0000 1.1
+++ WeightMatrix.java 9 Aug 2006 20:28:36 -0000 1.2
@@ -52,15 +52,17 @@
private boolean ignoreMS = false;
private DecimalFormat df = new DecimalFormat();
private Detector d;
- private Track tr;
+// private Track tr;
private Matrix wm;
private Matrix ermat;
+ private Matrix ermats;
private Matrix erm;
private Hep3Vector[] trkpnts;
private int[] hitlrs;
private HelixSwimmer swmr = null;
private int tchg;
private int naclrs=0;
+ private int nlrswh = 0;
private double Bfield = 5.;
private WMTrackPropagator mser = new WMTrackPropagator();
private CMTransform cmtr = new CMTransform();
@@ -101,7 +103,7 @@
public void setIncludeCorrelations(boolean incc) { inccor = incc; initialise(); }
public Hep3Vector[] getTrkPnts() {return trkPnts; }
public WMTrackPropagator getPropagator() { return mser; }
- public Matrix getErrorMatrix() { return ermat; }
+ public Matrix getErrorMatrix() { return ermats; }
public void setTrack(Track tr)
{
@@ -110,36 +112,31 @@
else
{
if(!initialized) initialise();
-// System.out.println("Track has "+thits.size()+" hits assigned");
mser.setTrack(track);
erm = mser.getProjectedErrorMatrix();
trkpnts = mser.getTrkActLrPnts();
hitlrs = mser.getTrkHitActLrs();
if(hitlrs == null) System.out.println("Track did not hit any layers !");
naclrs = mser.getNActiveLrs();
-// int m=erm.getRowDimension();
-// int n=erm.getColumnDimension();
-// System.out.println("WMTrackPropogator found naclrs "+naclrs+" with track inters: ");
-// for(int i=0; i<naclrs; i++)
-// {
-// System.out.println(" "+df.format(trkpnts[i].x())+" "+df.format(trkpnts[i].y())+" "+
-// df.format(trkpnts[i].z()));
-// }
}
}
+ public void setTrackParams(double[] prms)
+ {
+ mser.setTrack(prms);
+ }
+
public void assignHits()
{
wm = null;
- int nlrswh = 0;
- int mrank=2*naclrs;
+ nlrswh = 0;
if(naclrs < 3) return;
double[] lrsxra = new double[naclrs];
double minra = 10000.;
double maxra = 0.;
-// System.out.println("Setting WM with number of active lrs "+naclrs);
+ int mrank=2*naclrs;
Matrix msm = new Matrix(mrank,mrank);
- double averme = erm.get(mrank/4,mrank/4);
+// double averme = erm.get(mrank/4,mrank/4);
List<TrackerHit> thits = track.getTrackerHits();
for(int i=0; i<mrank; i++)
for(int j=0; j<mrank; j++) msm.set(i,j,0.);
@@ -162,7 +159,7 @@
// System.out.println("Searching hits close to radius "+df.format(rtr));
// System.out.println("Helix crossing at X: "+df.format(tx)+" Y: "+df.format(ty)+" Z: "+df.format(tz)+" radius: "+df.format(rtr));
// System.out.println("Expected position errors sq. from MS: "+df.format(exper1s)+" "+df.format(exper2s));
- double cldists = 10000.;
+ double cldists = 1000.;
TrackerHit clht=null;
active[i]=false;
closeHits[i] = null;
@@ -298,11 +295,11 @@
double norres = dphx*rt/Math.sqrt(c1);
// double norres = dphx*rt;
// System.out.println("Norm. residual for layer "+hitlrs[i]+" "+df.format(norres)+" cldists "+df.format(cldists));
- if(Math.abs(norres) < 5.)
+ if(Math.abs(norres) < 30.)
aida.cloud1D("WeightMatrix/ residual XY in lr "+hitlrs[i]).fill(norres);
norres = (tzr-hpo[2])/Math.sqrt(c2);
// norres = (tzr-hpo[2]);
- if(Math.abs(norres) < 5.)
+ if(Math.abs(norres) < 30.)
aida.cloud1D("WeightMatrix/ residual Z in lr "+hitlrs[i]).fill(norres);
}
}
@@ -340,11 +337,11 @@
double norres = dphx*rt/Math.sqrt(c1);
// double norres = dphx*rt;
// System.out.println("Norm. residual for layer "+hitlrs[i]+" "+df.format(norres)+" cldists "+df.format(cldists));
- if(Math.abs(norres) < 5.)
+ if(Math.abs(norres) < 30.)
aida.cloud1D("WeightMatrix/ residual XY in lr "+hitlrs[i]).fill(norres);
norres = (rh-rt)/Math.sqrt(c2);
// norres = (rh-rt);
- if(Math.abs(norres) < 5.)
+ if(Math.abs(norres) < 30.)
aida.cloud1D("WeightMatrix/ residual R in lr "+hitlrs[i]).fill(norres);
}
}
@@ -403,12 +400,50 @@
if(i != j) { ermat.set(i,j,0.); ermat.set(j,i,0.); }
}
}
- wm = ermat.inverse();
-// wm = msm.inverse();
+ nlrswh = 0;
+ for(int i=0; i<naclrs; i++)
+ if(active[i]) nlrswh++;
+ int mranks = 2*nlrswh;
+ ermats = new Matrix(mranks,mranks);
+ int in=0;
+ int jn=0;
+ for(int i=0; i<naclrs; i++)
+ {
+ jn=0;
+ for(int j=0; j<naclrs; j++)
+ {
+ if(active[i]&&active[j])
+ {
+ ermats.set(in,jn,ermat.get(i,j));
+ ermats.set(in+nlrswh,jn,ermat.get(i+naclrs,j));
+ ermats.set(in,jn+nlrswh,ermat.get(i,j+naclrs));
+ ermats.set(in+nlrswh,jn+nlrswh,ermat.get(i+naclrs,j+naclrs));
+ }
+ if(active[j]) jn++;
+ }
+ if(active[i]) in++;
+ }
+
+ wm = ermats.inverse();
+ in = 0;
+ for(int i=0; i<naclrs; i++)
+ {
+ if(active[i])
+ {
+ closeHits[in]=closeHits[i];
+ Cylinder[in]=Cylinder[i];
+ hitlrs[in]=hitlrs[i];
+ trkpnts[in]=trkPnts[i];
+ trkPnts[in]=trkPnts[i];
+ active[in]=true;
+ in++;
+ }
+ }
+ for(int i=nlrswh; i<MAXLRS; i++) active[i]=false;
boolean negel = false;
- for(int i=0; i<mrank; i++)
+ for(int i=0; i<mranks; i++)
{
- for(int j=0; j<mrank; j++)
+ for(int j=0; j<mranks; j++)
{
double mel = wm.get(i,j);
if(Double.toString(mel).equals("NaN"))
@@ -437,18 +472,15 @@
}
public void calculateTrkInt()
{
- for(int i=0; i< naclrs; i++)
+ for(int i=0; i< nlrswh; i++)
{
- double tx = trkpnts[i].x();
- double ty = trkpnts[i].y();
- double tz = trkpnts[i].z();
TrackerHit clht = closeHits[i];
double[] hpo = {0.,0.,0.};
- if(active[i]&&(clht!=null)) hpo=clht.getPosition();
+ if(clht!=null) hpo=clht.getPosition();
if(Cylinder[i])
{
- double rref=Math.sqrt(tx*tx+ty*ty);
- if(active[i]&&(clht!=null))
+ double rref=0.;
+ if(clht!=null)
{
rref=Math.sqrt(hpo[0]*hpo[0]+hpo[1]*hpo[1]);
}
@@ -456,14 +488,13 @@
double txr = ctrp.x();
double tyr = ctrp.y();
double tzr = ctrp.z();
- double rtin = Math.sqrt(txr*txr+tyr*tyr);
Hep3Vector trkpnt = new BasicHep3Vector(txr,tyr,tzr);
trkPnts[i]=trkpnt;
}
else
{
- double zref=tz;
- if(active[i]&&(clht!=null)) zref=hpo[2];
+ double zref=0.;
+ if(clht!=null) zref=hpo[2];
Hep3Vector ctrp = mser.getTrkPntAtZ(zref);
double txz = ctrp.x();
double tyz = ctrp.y();
@@ -473,10 +504,12 @@
}
}
}
- public int getNActiveLrs() { return naclrs; }
+ public int getNActiveLrs() { return nlrswh; }
public Matrix getWeightMatrix() { return wm; }
public boolean[] getActive() { return active; }
public TrackerHit[] getHits() { return closeHits; }
public boolean[] isCylinder() { return Cylinder; }
+
+ public int[] getHitLrs() { return hitlrs; }
}
\ No newline at end of file
CVSspam 0.2.8
