1 added + 5 modified, total 6 files
lcsim/src/org/lcsim/contrib/NickSinev/event/base
diff -N BaseSimTrackerHit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BaseSimTrackerHit.java 8 Sep 2006 00:51:06 -0000 1.1
@@ -0,0 +1,59 @@
+package org.lcsim.contrib.NickSinev.event.base;
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.event.*;
+import org.lcsim.event.base.*;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.Subdetector;
+
+public class BaseSimTrackerHit implements SimTrackerHit
+{
+ double[] point = new double[3];
+ MCParticle mcparticle;
+ double time=0.;
+ double dedx=0.;
+ int cellID = 0;
+ int layer = 0;
+ double[] momentum;
+ double pathlength = 0.;
+ IDDecoder iddecoder;
+ Subdetector subdetector;
+
+
+
+ public BaseSimTrackerHit(double[] pos,MCParticle mcp)
+ {
+ mcparticle = mcp;
+ for(int i=0; i<3; i++) point[i]=pos[i];
+ }
+
+
+ public void setPoint(double[] pnt)
+ {
+ for(int i=0; i<3; i++)
+ point[i]=pnt[i];
+ }
+
+ public void setMCParticle(MCParticle mcp) { mcparticle=mcp; }
+
+ public int getLayer() {return layer; }
+
+ public double[] getPoint() { return point; }
+
+ public double getTime() { return time; }
+
+ public double getdEdx() { return dedx; }
+
+ public MCParticle getMCParticle() { return mcparticle; }
+
+ public int getCellID() { return cellID; }
+
+ public IDDecoder getIDDecoder() { return iddecoder; }
+
+ public Subdetector getSubdetector() { return subdetector; }
+
+ public double getPathLength() { return pathlength; }
+
+ public double[] getMomentum() { return momentum; }
+
+}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/NickSinev/event/base
diff -u -r1.1 -r1.2
--- SmearedTrackerHit.java 2 Aug 2006 01:29:46 -0000 1.1
+++ SmearedTrackerHit.java 8 Sep 2006 00:51:06 -0000 1.2
@@ -36,8 +36,8 @@
double phi = Math.atan2(_pos[1],_pos[0]);
double R = Math.sqrt(_pos[0]*_pos[0]+_pos[1]*_pos[1]);
double Z = _pos[2];
- double phis = phi+RPsm/R;
double Rs = R+Rsm;
+ double phis = phi+RPsm/Rs;
double Zs = Z+Zsm;
_pos[0]=Rs*Math.cos(phis);
_pos[1]=Rs*Math.sin(phis);
lcsim/src/org/lcsim/contrib/NickSinev/event/util
diff -u -r1.2 -r1.3
--- SmearTrackerHits.java 9 Aug 2006 20:28:35 -0000 1.2
+++ SmearTrackerHits.java 8 Sep 2006 00:51:10 -0000 1.3
@@ -23,7 +23,7 @@
/**
* @author sinev U of Oregon; SLAC x2970; [log in to unmask]
- * @version $Id: SmearTrackerHits.java,v 1.2 2006/08/09 20:28:35 sinev Exp $
+ * @version $Id: SmearTrackerHits.java,v 1.3 2006/09/08 00:51:10 sinev Exp $
*/
public class SmearTrackerHits extends Driver
{
@@ -61,25 +61,28 @@
CMTransform cmt = new CMTransform();
private AIDA aida = AIDA.defaultInstance();
private boolean doSmear=true;
- private boolean doHist = false;
+ private boolean doHist = false;
+ private static SmearTrackerHits instance;
public SmearTrackerHits()
{
+ instance = this;
}
- /**
- * process method called every time new event has been read
- * @param event <code>EventHeader</code> object containing list of all lists in the event
- */
+ public static SmearTrackerHits getInstance()
+ {
+ if(instance == null) instance = new SmearTrackerHits();
+ return instance;
+ }
- public void setVtxBarrRPResoloution(double res) { vtxbarRPhism = res;}
- public void setVtxBarrZResoloution(double res) { vtxbarZsm = res;}
- public void setVtxECRPResoloution(double res) { vtxecRPhism = res;}
- public void setVtxEcRResoloution(double res) { vtxecRsm = res;}
- public void setTrkBarrRPResoloution(double res) { trkbarRPhism = res;}
- public void setTrkBarrZResoloution(double res) { trkbarZsm = res;}
- public void setTrkECRPResoloution(double res) { trkecRPhism = res;}
- public void setTrkEcRResoloution(double res) { trkecRsm = res;}
+ public void setVtxBarrRPResolution(double res) { vtxbarRPhism = res;}
+ public void setVtxBarrZResolution(double res) { vtxbarZsm = res;}
+ public void setVtxECRPResolution(double res) { vtxecRPhism = res;}
+ public void setVtxECRResolution(double res) { vtxecRsm = res;}
+ public void setTrkBarrRPResolution(double res) { trkbarRPhism = res;}
+ public void setTrkBarrZResolution(double res) { trkbarZsm = res;}
+ public void setTrkECRPResolution(double res) { trkecRPhism = res;}
+ public void setTrkECRResolution(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; }
@@ -89,7 +92,20 @@
public void setTrkBarrHitClustTransSpan(double tsp) { btrkClSpanCS = tsp; }
public void setTrkECHitClustRZpan(double zsp) { etrkClSpanZ = zsp; }
public void setTrkECHitClustTransSpan(double tsp) { etrkClSpanDS = tsp; }
+ public double getVtxBarrRPResolution() { return vtxbarRPhism; }
+ public double getVtxBarrZResolution() { return vtxbarZsm;}
+ public double getVtxECRPResolution() { return vtxecRPhism;}
+ public double getVtxECRResolution() { return vtxecRsm;}
+ public double getTrkBarrRPResolution() { return trkbarRPhism;}
+ public double getTrkBarrZResolution() { return trkbarZsm;}
+ public double getTrkECRPResolution() { return trkecRPhism;}
+ public double getTrkECRResolution() { return trkecRsm;}
+ /**
+ * process method called every time new event has been read
+ * @param event <code>EventHeader</code> object containing list of all lists in the event
+ */
+
public void process(EventHeader event)
{
double[] oldcor = new double[3];
lcsim/src/org/lcsim/contrib/NickSinev/tracking/util
diff -u -r1.2 -r1.3
--- WMTrackPropagator.java 9 Aug 2006 20:28:35 -0000 1.2
+++ WMTrackPropagator.java 8 Sep 2006 00:51:12 -0000 1.3
@@ -27,7 +27,7 @@
/**
* @author Nick Sinev U of Oregon, [log in to unmask] SLAC x2970
- * @version $Id: WMTrackPropagator.java,v 1.2 2006/08/09 20:28:35 sinev Exp $
+ * @version $Id: WMTrackPropagator.java,v 1.3 2006/09/08 00:51:12 sinev Exp $
*/
public class WMTrackPropagator
@@ -78,8 +78,11 @@
double hcenterY=0.;
double trackPt = 0.;
double trsinl = 0.;
+ double trcosl = 0.;
protected Matrix errorMatrix;
protected Matrix projErrorMatrix;
+ protected Matrix corrMatrix;
+ protected Matrix errMatZ;
boolean initialized = false;
int nSubLrs=0;
int nhlrs = 0;
@@ -159,6 +162,8 @@
public int[] getTrkHitActLrs() {return actLrInd;}
+ public int[] getActLrsSubDet() {return subDetInd; }
+
public Hep3Vector[] getTrkActLrPnts() {return actLrPnts; }
public Hep3Vector[] getTrkMomAtLrs() {return trkMomAt; }
@@ -471,8 +476,8 @@
}
if(!iscyl)
{
- minR[slind]=rmin[i];
- maxR[slind]=rmax[i];
+ minR[slind]=rmin[i]; // blure rings edges
+ maxR[slind]=rmax[i]; // same
if(detailed)
{
minZ[slind]=zcur;
@@ -594,6 +599,7 @@
trkfulmom = trackPt * Math.sqrt(1. + ttl*ttl);
trkCurvRad = Math.abs(1./trpar[2]);
trsinl = ttl/Math.sqrt(1.+ttl*ttl);
+ trcosl = 1./Math.sqrt(1.+ttl*ttl);
// System.out.println("WMTrackPropagator is processing track with Pt "+trackPt+" full mom "+trkfulmom);
double maxlen = 100000.;
iq=0;
@@ -605,7 +611,7 @@
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 > 2.*Math.PI*trkCurvRad) maxlen = 2.*Math.PI*trkCurvRad;
if(maxlen > 3.* maxz) maxlen = 3.* maxz;
// System.out.println("Max trk length: "+df.format(maxlen));
hitbytr = new boolean[nSubLrs];
@@ -626,7 +632,7 @@
if(Math.abs(pnt.z()) < maxZ[sl])
{
// System.out.println("Sublayer "+sl+" is hit at Z "+df.format(pnt.z())+
-// " length along tr: "+df.format(len));
+// " length along tr: "+df.format(len)+" radius "+df.format(minR[sl]));
hitbytr[sl]=true;
lentoh[sl]=len;
entryCyl[sl]=true;
@@ -639,7 +645,7 @@
boolean hit = true;
len = swm.getDistanceToZ(ttls*minZ[sl]);
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(2.*Math.PI*trkCurvRad*Math.abs(ttl) < Math.abs(minZ[sl])) hit=false;
if(Math.abs(len) > maxlen) hit=false;
if(hit)
{
@@ -919,7 +925,8 @@
if(activeli[i] != -1)
{
actLrPnts[activeli[i]] = new BasicHep3Vector(swc.x(),swc.y(),swc.z());
- actIsCyl[activeli[i]]=entcyl[i];
+// actIsCyl[activeli[i]]=entcyl[i];
+ actIsCyl[activeli[i]]=isCylinder[subLrInd[actLrInd[activeli[i]]]];
Trajectory trj = swm.getTrajectory();
Hep3Vector swmd = ((Helix) trj).getTangentAtDistance(oslel[i]);
double swmma = Math.sqrt(swmd.x()*swmd.x()+swmd.y()*swmd.y()+swmd.z()*swmd.z());
@@ -1012,6 +1019,8 @@
}
naclt = ncntactl;
errorMatrix = new Matrix(naclt,naclt);
+ corrMatrix = new Matrix(naclt,naclt);
+ errMatZ = new Matrix(naclt,naclt);
double[] lngths = new double[naclt];
// System.out.println("Total number of active layers for this track: "+naclt);
for(int i=0; i<naclt; i++)
@@ -1033,6 +1042,9 @@
// double lenj = chordl[indinhl[posj]];
double lenj = lentoesp[posj];
double mel = 0.;
+ double melcoxy = 0.;
+ double melcoz = 0.;
+ double melz = 0.;
for(int k=0; k<posj; k++)
{
int kk = indinhl[k];
@@ -1055,15 +1067,20 @@
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);
+ mel+= (scinesp[k]*scinesp[k]*levj*levi);
+ melz+= (scinesp[k]*scinesp[k]*levj*levi);
+// melcoz+= scinesp[k]*levi*scinesp[k]*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));
}
+ corrMatrix.set(i,j,melcoz);
+ corrMatrix.set(j,i,melcoz);
if(j==i)
{
errorMatrix.set(i,j, mel);
+ errMatZ.set(i,j, melz);
}
else
{
@@ -1071,27 +1088,33 @@
{
errorMatrix.set(j,i, mel); //+++++++++++++++++++++++++++++++++++++++++++
errorMatrix.set(i,j, mel); //+++++++++++++++++++++++++++++++++++++++++++
+ errMatZ.set(i,j, melz);
+ errMatZ.set(j,i, melz);
}
if(!inccor)
{
errorMatrix.set(j,i,0.);
errorMatrix.set(i,j,0.);
+ errMatZ.set(i,j,0.);
+ errMatZ.set(j,i,0.);
}
}
}
}
boolean badm = false;
- for(int i=0; i<naclt; i++)
+ for(int i=1; i<naclt; i++)
{
- for(int j=0; j<i+1; j++)
+ for(int j=0; j<i; j++)
{
- if(errorMatrix.get(i,j) > errorMatrix.get(i,i))
+ if((errorMatrix.get(i,j) - errorMatrix.get(i,i))>0.00000001)
{
System.out.println("Non-diag element: "+i+" "+j+" "+df.format(errorMatrix.get(i,j))+" larger then diag: "+df.format(errorMatrix.get(i,i)));
badm =true;
}
}
}
+// errorMatrix.print(8,6);
+// corrMatrix.print(8,6);
if(badm)
{
errorMatrix.print(8,6);
@@ -1107,7 +1130,11 @@
int nmeast = naclt*2;
projErrorMatrix = new Matrix(nmeast,nmeast);
double[] corco1 = new double[naclt];
- double[] corco2 = new double[naclt];
+ double corco2c = 1./trcosl;
+ double corco2d = 1./Math.abs(trsinl);
+ Hep3Vector[] dirm1 = new Hep3Vector[naclt];
+ Hep3Vector[] dirm2 = new Hep3Vector[naclt];
+ double[] rproj = new double[naclt];
double sigi=1.;
double sigj=1.;
double sigcor = 1.;
@@ -1115,95 +1142,91 @@
{
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 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());
+ dirm1[i]=dir1i;
Hep3Vector trkchdir = trkMomAt[i];
Hep3Vector trkdir = VecOp.unit(trkchdir);
-// if(i>0) trkchdir = VecOp.sub(actLrPnts[i], actLrPnts[i-1]);
+ Hep3Vector pltrdir = new BasicHep3Vector(trkchdir.x(),trkchdir.y(),0.);
+ Hep3Vector pltrk = VecOp.unit(pltrdir);
sigi = 1.;
+ Hep3Vector dir2inn = new BasicHep3Vector(actLrPnts[i].x(),actLrPnts[i].y(),0.);
+ Hep3Vector dir2i = VecOp.unit(dir2inn);
if(!isCylinder[lrini])
{
- 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)/
- Math.sqrt(VecOp.dot(dir1i,dir1i)*VecOp.dot(dirsc1,dirsc1)));
- double cosp2i = Math.abs(VecOp.dot(dir2i,dirsc2)/
- Math.sqrt(VecOp.dot(dir2i,dir2i)*VecOp.dot(dirsc2,dirsc2)));
- if(Double.toString(cosp1i).equals("NaN")) cosp1i = 0.1;
- if(Double.toString(cosp2i).equals("NaN")) cosp2i = 0.1;
- if(cosp1i < 0.1) cosp1i = 0.1;
- if(cosp2i < 0.1) cosp2i = 0.1;
- double cor1i = 1./cosp1i;
- double cor2i = 1./cosp2i; */
+ dirm2[i]=dir2i;
double cosd1 = VecOp.dot(trkdir, dir1i);
+ rproj[i]=1.;
+ if(!isCylinder[lrini]) rproj[i] = VecOp.dot(pltrk, dir2i);
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])
- corco1[i]=corco1[i-1];
- if(corco2[i]<corco2[i-1])
- corco2[i]=corco2[i-1];
- }
}
for(int i=0; i<naclt; i++)
{
+ double corfr = 1.0;
+ double corfr2 = 1.;
int ii = actLrInd[i];
int lrini = subLrInd[ii];
double cor1i = corco1[i];
- double cor2i = corco2[i];
+ double cor2i = corco2c;
sigi = 1.;
if(!isCylinder[lrini])
{
sigi = -1.;
+ cor1i=1.;
+ cor2i=corco2d;
+ corfr = 1./corco1[i];
+ corfr2 = rproj[i];
}
for(int j=0; j<naclt; j++)
{
int jj = actLrInd[j];
int lrinj = subLrInd[jj];
sigj = 1.;
+ double cor1j = corco1[j];
+ double cor2j = corco2c;
if(!isCylinder[lrinj])
{
sigj=-1.;
- }
- double cor1j = corco1[j];
- double cor2j = corco2[j];
+ cor1j = 1.;
+ cor2j = corco2d;
+ corfr = 1./corco1[j];
+ corfr2 = rproj[j];
+ }
+ if(sigi*sigj > 0.) { corfr = 1.; corfr2 = 1.; }
double mel = errorMatrix.get(i,j)*cor1i*cor1j; //++++++++++++++++++++++++++++++++
+ if(i != j) mel*=corfr;
projErrorMatrix.set(i,j, mel);
- mel = errorMatrix.get(i,j)*cor2i*cor2j;
+ mel = errMatZ.get(i,j)*cor2i*cor2j;
+ if(i!=j) mel*= corfr2;
sigcor= sigi*sigj;
if(sigcor<0) sigcor*= ttls;
if(i==j) sigcor = 1.;
projErrorMatrix.set(i+naclt, j+naclt,sigcor*mel);
+ double siginco = -1.*iq*ttls;
+/* if(sigi < 0) sigi = sigi*(-1.*ttls);
+ if(sigj < 0) sigj = sigj*(-1.*ttls);
+ sigi = sigi*siginco;
+ sigj = sigj*siginco; */
+ sigcor*=siginco;
+// sigcor = 1.; //++++++++++++++++++++++++++++++++++++++++++!!! Ignore sign !
+/* if(true)
+ {
+ mel = corrMatrix.get(i,j)*cor1i*cor2j*sigcor;
+ projErrorMatrix.set(i+naclt,j,mel);
+ projErrorMatrix.set(i,j+naclt,mel);
+ } */
}
}
boolean badm = false;
- for(int i=0; i<naclt; i++)
+/* for(int i=1; i<naclt; i++)
{
- for(int j=0; j<i+1; j++)
+ for(int j=0; j<i; j++)
{
if(Double.toString(projErrorMatrix.get(i,j)).equals("NaN")) badm=true;
- if(projErrorMatrix.get(i,j) > projErrorMatrix.get(i,i))
+ if((projErrorMatrix.get(i,j) - projErrorMatrix.get(i,i))>0.0000001)
{
System.out.println("Non-diag element: "+i+" "+j+" "+
df.format(projErrorMatrix.get(i,j))+" is larger then diag: "+df.format(projErrorMatrix.get(i,i)));
@@ -1215,25 +1238,11 @@
badm=true;
}
}
- }
+ } */
if(badm)
{
- errorMatrix.print(8,4);
+ projErrorMatrix.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.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)
@@ -1295,7 +1304,11 @@
{
mddr = Math.abs(ral-r);
icll=i;
- if(ral > r) icll = i-1;
+ if((ral-r)>0.000001) icll=i-1;
+ if(Math.abs(ral-r) <= 0.000001)
+ {
+ return new BasicHep3Vector(pos.x(), pos.y(), pos.z());
+ }
}
}
}
@@ -1336,7 +1349,6 @@
}
}
setSwimmerAtPoint(icll);
- setSwimmerAtPoint(icll);
Trajectory trj = swm.getTrajectory();
double leno = trj.getDistanceToZPlane(z);
Hep3Vector npos = trj.getPointAtDistance(leno);
lcsim/src/org/lcsim/contrib/NickSinev/tracking/wmfitter
diff -u -r1.2 -r1.3
--- SLDFitTrack.java 9 Aug 2006 20:28:36 -0000 1.2
+++ SLDFitTrack.java 8 Sep 2006 00:51:15 -0000 1.3
@@ -15,15 +15,15 @@
import hep.physics.vec.VecOp;
import hep.physics.vec.BasicHep3Vector;
import java.util.List;
-import org.lcsim.event.Track;
+import org.lcsim.event.*;
import org.lcsim.event.base.*;
-import org.lcsim.event.TrackerHit;
import org.lcsim.util.swim.*;
import org.lcsim.geometry.FieldMap;
import org.lcsim.material.MaterialManager;
import org.lcsim.constants.Constants;
import org.lcsim.geometry.Detector;
import org.lcsim.contrib.NickSinev.tracking.util.*;
+import org.lcsim.util.aida.AIDA;
public class SLDFitTrack
@@ -77,12 +77,15 @@
private double[] delpar = new double[5];;
private DecimalFormat df = new DecimalFormat();
+ private DecimalFormat dfr = new DecimalFormat();
private boolean printRes = false;
double[] initpar = new double[6];
double[] finalpar = new double[6];
// WMTrackPropagator trkPr=null;
private boolean onlyOmegaFit = false;
private boolean oripar = true;
+ private AIDA aida = AIDA.defaultInstance();
+ private boolean doHist = true;
@@ -91,10 +94,12 @@
{
IPConstrained = false;
planarFit=false;
+ dfr.setMaximumFractionDigits(1);
}
public void setConstrainedFit(boolean set) { IPConstrained = set; }
public void setPlanarFit(boolean set) { planarFit=set; }
+ public void setDoHist(boolean doh) {doHist=doh;}
public void setDetector(Detector d)
{
@@ -123,7 +128,24 @@
public void fit()
{
- assignHits();
+ assignHits();
+ fitHits();
+ }
+
+ public void simulateFit(BaseTrack tr, MCParticle mcp)
+ {
+ setTrack(tr);
+ simulateFit(mcp);
+ }
+
+ public void simulateFit(MCParticle mcp)
+ {
+ simulateHits(mcp);
+ fitHits();
+ }
+
+ public void fitHits()
+ {
if(NDoF < 3)
{
System.out.println("Too little hits on track!");
@@ -203,17 +225,25 @@
tchg = (tkpar[2] > 0.)? 1 : -1;
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; }
}
public void assignHits()
{
- Wm.assignHits();
+ Wm.assignHits();
+ assignHits_c();
+ }
+
+ public void simulateHits(MCParticle mcp)
+ {
+ Wm.simulateHits(mcp);
+ assignHits_c();
+ }
+
+ private void assignHits_c()
+ {
if(debug) System.out.println("Extracting weight matrix\n");
- Matrix mtr = Wm.getWeightMatrix();
- if(mtr == null)
+ wm = Wm.getWeightMatrix();
+ if(wm == null)
{
NDoF = 0;
return;
@@ -226,9 +256,9 @@
{
for(int j=0; j<2*naclrs; j++)
{
- wma[i][j]=mtr.get(i,j);
+ wma[i][j]=wm.get(i,j);
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);
+// if((i==j) && (wma[i][j] < 0.)) System.out.println("Negative diagonal element in weight matrix: "+i+","+j);
}
}
if(badTrack)
@@ -244,12 +274,6 @@
for(int i=0; i<naclrs; i++)
{
hasHit[i] = Wm.getActive()[i];
- TrackerHit hit=Wm.getHits()[i];
- if(hit != null)
- {
- double[] hpo = hit.getPosition();
- double rhit = Math.sqrt(hpo[0]*hpo[0]+hpo[1]*hpo[1]);
- }
if(hasHit[i]) nlrswh++;
}
NDoF = nlrswh * 2 - nFreePar;
@@ -318,6 +342,10 @@
if(Math.abs(zH[i]-zX[i]) > 0.001) System.out.println("Disk: wrong inters. z: "+zX[i]+
"while hit is at z "+zH[i]);
}
+ if((nitt==0)&& doHist)
+ {
+ System.out.println("Layer "+i+" int. rad "+df.format(rX[i])+" Z "+df.format(zX[i]));
+ }
}
if(i>=ntpnts) hasHit[i]=false;
}
@@ -326,7 +354,7 @@
private void getHitPositions()
{
- Wm.setTrackParams(tkpar);
+// Wm.setTrackParams(tkpar);
TrackerHit hit;
double[] hpo;
for(int i=0; i<naclrs; i++)
@@ -369,15 +397,19 @@
private void printResiduals(double[] res)
{
+ df.setMaximumFractionDigits(8);
for(int i=0; i<naclrs; i++)
{
+ double phiHC = phiH[i];
+ double phiXC = phiX[i];
double xH=rX[i]*Math.cos(phiH[i]);
double yH=rX[i]*Math.sin(phiH[i]);
double xX=rX[i]*Math.cos(phiX[i]);
double yX=rX[i]*Math.sin(phiX[i]);
- System.out.println("Lr "+i+" R "+df.format(rX[i])+" xH "+df.format(xH)+" yH "+df.format(yH)+
- " xX "+df.format(xX)+" yX "+df.format(yX)+" res "+df.format(res[i]));
+ System.out.println("Lr "+i+" R "+df.format(rX[i])+" phiH "+df.format(phiHC)+
+ " phiX "+df.format(phiXC)+" res "+df.format(res[i]));
}
+ df.setMaximumFractionDigits(6);
}
@@ -407,20 +439,25 @@
{
if(nFreePar==5)
{
- derivatives[j][1] = rX[j];
+// derivatives[j][1] = rX[j];
// derivatives[j][2] = -0.5 * rX[j]*rX[j];
+// derivatives[j][3] = 0.;
+// derivatives[j][4] = 0.;
}
- if((nFreePar == 3)&&!planarFit)
- {
- derivatives[j][0] = rX[j];
-// derivatives[j][1] = -0.5 * rX[j]*rX[j];
- }
- if((nFreePar == 3)&&planarFit)
+ }
+ else
+ {
+ if(planarFit)
{
- derivatives[j][1] = rX[j];
-// derivatives[j][2] = -0.5 * rX[j]*rX[j];
+ derivatives[j][1]=0.;
+ derivatives[j][2]=0.;
+ }
+ if((nFreePar == 3)&&!planarFit)
+ {
+ derivatives[j][0]=0.;
+ derivatives[j][1]=0.;
}
- }
+ }
}
freePar[i]-=delpar[i];
}
@@ -433,10 +470,24 @@
derivatives[j+naclrs][0]=0.;
}
}
+/* if(nitt==0)
+ {
+ System.out.println("N active layers: "+naclrs);
+ for(int i=0; i<2*naclrs; i++)
+ {
+ System.out.print(i+" ");
+ for(int j=0; j<nFreePar; j++)
+ {
+ System.out.print(df.format(derivatives[i][j])+" : ");
+ }
+ System.out.println(" ");
+ }
+ } */
}
private void getSolution()
{
+ df.setMaximumFractionDigits(5);
double ym=0.;
int wml = wma[0].length;
int nps = 2*naclrs;
@@ -478,7 +529,55 @@
Ym[m] = -ym;
}
wp = new Matrix(WPmn);
-// if(debug) wp.print(8,6);
+ if(doHist)
+ {
+ System.out.println("Itteration "+nitt+ " wp matrix: ");
+ wp.print(8,6);
+ }
+ if(doHist && (nitt==0))
+ {
+ for(int i=0; i<nFreePar; i++)
+ {
+ double pwt = Math.sqrt(wp.get(i,i));
+ aida.cloud1D("SLDFitTrack/ parameter "+i+" weight").fill(pwt);
+ if(i==2)
+ {
+ if(true)
+ {
+// System.out.println("Normal weights");
+// wp.print(8, 1);
+ System.out.println(
+ "n cyl? Xrad Z derivat diag wt diag cont full contr running sum");
+ double mes = 0.;
+ for(int ip=0; ip<nps; ip++)
+ {
+ double me = 0.;
+ int lin = 0;
+ if(ip < naclrs) lin = ip;
+ if(ip > naclrs-1) lin = ip - naclrs;
+ for(int j=0; j<nps; j++)
+ {
+ int ljn = 0;
+ if(j<naclrs) ljn = j;
+ if(j>naclrs-1) ljn = j-naclrs;
+ if(hasHit[lin] && hasHit[ljn])
+ {
+ me+=derivatives[ip][i]*wma[ip][j]*derivatives[j][i];
+ }
+ }
+ mes+=me;
+ double dime = derivatives[ip][i]*wma[ip][ip]*derivatives[ip][i];
+ System.out.println(ip+" "+isCylinder[lin]+" "+dfr.format(rX[lin])+" "+dfr.format(zX[lin])+
+ " "+df.format(derivatives[ip][i])+" "+dfr.format(wma[ip][ip])+
+ " "+dfr.format(dime)+" "+dfr.format(me)+" "+dfr.format(mes));
+
+ }
+
+
+ }
+ }
+ }
+ }
int nrs = nFreePar;
Matrix y = new Matrix(Ym,nrs);
// if(debug) y.print(8,6);
@@ -539,17 +638,18 @@
}
}
}
- if((chi2<0.) && (nlrswh==naclrs))
+ if(doHist && (nlrswh==naclrs) && (nitt==0))
{
- System.out.println("Negative chi2 ! Weight matrix is: ");
- for(int i=0; i<nps; i++)
- {
- for(int j=0; j<nps; j++)
- System.out.print(df.format(wma[i][j])+" ");
- System.out.println(" ");
- }
+ if(chi2 < 0.) System.out.println("Negative chi2 ! Weight matrix is: ");
+ else System.out.println("Weight matrix is: ");
+ Matrix prmatr = wm.getMatrix(0,naclrs-1,0,naclrs-1);
+ prmatr.print(8,6);
System.out.println("Error matrix was: ");
- ermx.print(8,6);
+ prmatr = ermx.getMatrix(0,naclrs-1,0, naclrs-1);
+ prmatr.print(8,6);
+// Matrix invpr = prmatr.inverse();
+// System.out.println(" and if we inverse it once more: ");
+// invpr.print(8,6);
}
if((chi2 > veryLarge-1.) ||(Double.toString(chi2).equals("NaN")))
{
@@ -608,23 +708,24 @@
success = false;
badTrack = false;
df.setMaximumFractionDigits(6);
- delpar[0] = 0.00002;
- delpar[1] = 0.000002;
+ delpar[0] = 0.0002;
+ delpar[1] = 0.0002;
delpar[2] = 0.000001;
delpar[3] = 0.00002;
- delpar[4] = 0.000002;
+ delpar[4] = 0.00002;
if((nFreePar == 3)&&(!planarFit))
{
- delpar[2] = 0.000002;
- delpar[1] = 0.000001;
+ delpar[2] = 0.0002;
+ delpar[1] = 0.0001;
}
- for(nitt = 0; nitt<4; nitt++)
+ for(nitt = 0; nitt<6; nitt++)
{
if(debug) printRes = true;
setTrackPars();
getTrkInters();
getResiduals(residuals);
calculateChi2();
+ if(chi2<-0.001) System.out.println("Negative chi2 "+chi2+" at itteration "+nitt);
if(nitt == 0)
{
initpar[5]=chi2;
@@ -645,19 +746,18 @@
{
double maxdev = 0.1;
if((i==0) || (i==3)) maxdev = 3.;
- if(Math.abs(dpar[i]) > 0.001 * delpar[i]) nochange = false;
-// if((nFreePar==5) && (Math.abs(dpar[2]) > 0.02*Math.abs(initpar[2]))) dpar[2]=0.1*dpar[2];
-// if((nFreePar==3) && (Math.abs(dpar[1]) > 0.02*Math.abs(initpar[2]))) dpar[1]=0.1*dpar[1];
+ if(Math.abs(dpar[i]) > 0.01 * delpar[i]) nochange = false;
if(Math.abs(dpar[i]) > maxdev)
{
badTrack = true;
- System.out.println(" FitTrack: too large change in parameter "+i+" dpar="+dpar[i]);
+// System.out.println(" FitTrack: too large change in parameter "+i+" dpar="+dpar[i]);
return;
}
freePar[i]+=dpar[i];
}
if(nochange) break;
}
+ if(nitt > 4) {badTrack=true; }
setTrackPars();
getTrkInters();
getResiduals(residuals);
@@ -669,7 +769,7 @@
{
System.out.println("final chi2= "+df.format(chi2));
}
- if(NDoF > 1.)
+ if((NDoF > 1.)&& !badTrack)
{
if(Math.abs(chi2/NDoF) < GoodChi2) success = true;
else
@@ -680,7 +780,9 @@
if(success)
{
erma = wp.inverse();
- }
+ if(doHist) erma.print(14,12);
+ }
+ else erma = null;
}
}
lcsim/src/org/lcsim/contrib/NickSinev/tracking/wmfitter
diff -u -r1.2 -r1.3
--- WeightMatrix.java 9 Aug 2006 20:28:36 -0000 1.2
+++ WeightMatrix.java 8 Sep 2006 00:51:16 -0000 1.3
@@ -16,29 +16,28 @@
package org.lcsim.contrib.NickSinev.tracking.wmfitter;
import org.lcsim.contrib.NickSinev.tracking.util.*;
-import java.util.Collection;
-import java.util.List;
-import java.util.Map;
-import java.util.Set;
+import java.util.*;
import java.text.*;
import Jama.*;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
import hep.physics.vec.BasicHep3Vector;
import java.util.List;
-import org.lcsim.event.Track;
-import org.lcsim.event.TrackerHit;
+import org.lcsim.event.*;
import org.lcsim.util.swim.*;
import org.lcsim.geometry.FieldMap;
import org.lcsim.material.MaterialManager;
import org.lcsim.constants.Constants;
import org.lcsim.geometry.Detector;
import org.lcsim.util.aida.AIDA;
+import org.lcsim.contrib.NickSinev.event.base.*;
+import org.lcsim.contrib.NickSinev.event.util.*;
public class WeightMatrix
{
public final int MAXLRS = 30;
public boolean debug = false;
+ int debug_level = 0;
private boolean doHist = true;
private boolean equal_wts = false;
protected boolean[] active = new boolean[MAXLRS];
@@ -51,6 +50,7 @@
private boolean doeloss = true;
private boolean ignoreMS = false;
private DecimalFormat df = new DecimalFormat();
+ private DecimalFormat dfe = new DecimalFormat("+0.000E00;-0.000E00");
private Detector d;
// private Track tr;
private Matrix wm;
@@ -59,24 +59,32 @@
private Matrix erm;
private Hep3Vector[] trkpnts;
private int[] hitlrs;
+ private int[] masklrs = new int[MAXLRS];
+ private int nmasked;
private HelixSwimmer swmr = null;
private int tchg;
private int naclrs=0;
private int nlrswh = 0;
+ private int nmissing = 0;
private double Bfield = 5.;
private WMTrackPropagator mser = new WMTrackPropagator();
private CMTransform cmtr = new CMTransform();
private AIDA aida = AIDA.defaultInstance();
private Track track;
+ Random ran = new Random();
+ CMTransform cmt = new CMTransform();
public WeightMatrix()
{
initialized=false;
+ nmasked = 0;
df.setMaximumFractionDigits(5);
}
public void setEqualWeights(boolean eq) {equal_wts=eq;}
+ public void setDebugLevel(int dbl) { debug_level = dbl; }
+
public void setDetector(Detector det)
{
d=det;
@@ -103,7 +111,10 @@
public void setIncludeCorrelations(boolean incc) { inccor = incc; initialise(); }
public Hep3Vector[] getTrkPnts() {return trkPnts; }
public WMTrackPropagator getPropagator() { return mser; }
- public Matrix getErrorMatrix() { return ermats; }
+ public Matrix getErrorMatrix() { return ermats; }
+ public int getNMissingHits() { return nmissing; }
+ public void unmaskAll(){ nmasked = 0;}
+ public void maskLayer(int lrn) {masklrs[nmasked]=lrn; nmasked++; }
public void setTrack(Track tr)
{
@@ -128,6 +139,7 @@
public void assignHits()
{
+ if(debug_level==1) System.out.println("Assigning hits!");
wm = null;
nlrswh = 0;
if(naclrs < 3) return;
@@ -138,16 +150,38 @@
Matrix msm = new Matrix(mrank,mrank);
// double averme = erm.get(mrank/4,mrank/4);
List<TrackerHit> thits = track.getTrackerHits();
+ int thnhts = thits.size();
+ boolean[] htassnd = new boolean[thnhts];
+ for(int i=0; i<thnhts; i++)
+ htassnd[i]=false;
for(int i=0; i<mrank; i++)
for(int j=0; j<mrank; j++) msm.set(i,j,0.);
+ double resp05 = 0.;
+ double resp14 = 0.;
+ double resr05 = 0.;
+ double resr14 = 0.;
+ double resr23 = 0.;
+ boolean havh05 = false;
+ boolean havh14 = false;
+ boolean havh23 = false;
for(int i=0; i<naclrs; i++)
{
// System.out.println("Looking for close hits in layer index "+i+" which is layer "+hitlrs[i]);
// if(mser.isActiveCylinder(i)) System.out.println("Propagator declared layer cylindrical");
// if(!mser.isActiveCylinder(i)) System.out.println("Propagator declared layer disk");
+ boolean masked = false;
+ for(int j=0; j<nmasked; j++)
+ {
+ if(hitlrs[i] == masklrs[j]) masked = true;
+ }
+ active[i]=false;
+ closeHits[i] = null;
+ if(!masked)
+ {
Cylinder[i]=mser.isActiveCylinder(i);
double exper1s = erm.get(i,i);
double exper2s = erm.get(i+naclrs,i+naclrs);
+ double expcorers = erm.get(i,i+naclrs);
double tx = trkpnts[i].x();
double ty = trkpnts[i].y();
double tz = trkpnts[i].z();
@@ -161,27 +195,34 @@
// System.out.println("Expected position errors sq. from MS: "+df.format(exper1s)+" "+df.format(exper2s));
double cldists = 1000.;
TrackerHit clht=null;
- active[i]=false;
- closeHits[i] = null;
double[] hpo=null;
+ int hin = 0;
+ int clhin = -1;
for(TrackerHit th : thits)
{
boolean clrad = false;
boolean clsz = false;
+ hin++;
+ boolean disk = false;
+ if(th.getCovMatrix()[2] < 0.0000000001) disk=true;
+ if(!htassnd[hin-1])
+ {
hpo = th.getPosition();
double rht = Math.sqrt(hpo[0]*hpo[0]+hpo[1]*hpo[1]);
// System.out.println("See hit with radius "+df.format(rht));
- if(Math.abs(rtr-rht) < 1.)
+ if((Math.abs(rtr-rht) < 0.3)&& mser.isActiveCylinder(i))
{
// System.out.println("Hit with close rad: "+df.format(rht)+" X: "+df.format(hpo[0])+" Y: "+df.format(hpo[1])+
// " Z: "+df.format(hpo[2]));
clrad = true;
if(Math.abs(tz-hpo[2])<100.) clsz = true;
}
- if(Math.abs(tz-hpo[2])<1.) clsz=true;
+ if(!mser.isActiveCylinder(i))
+ {
+ if(Math.abs(tz-hpo[2])<0.3) clsz=true;
+ if(Math.abs(rtr-rht) < 200.) clrad = true;
+ }
double[] tcv=null;
- boolean disk = false;
- if(th.getCovMatrix()[2] < 0.0000000001) disk=true;
if(clsz)
{
// if(disk) System.out.println("hit belongs to disk surface");
@@ -217,7 +258,7 @@
if(dpsi < 0.) dphx+=2.*Math.PI;
}
}
- if(disk && (!mser.isActiveCylinder(i))&& clsz)
+ if(disk && (!mser.isActiveCylinder(i))&& clsz && clrad )
{
Hep3Vector ctrp = mser.getTrkPntAtZ(hpo[2]);
double txz = ctrp.x();
@@ -226,20 +267,30 @@
Hep3Vector trkpnt = new BasicHep3Vector(txz,tyz,tzz);
trkPnts[i]=trkpnt;
double rh = Math.sqrt(hpo[0]*hpo[0]+hpo[1]*hpo[1]);
+ double phih = Math.atan2(hpo[1], hpo[0]);
+ double phit = Math.atan2(tyz, txz);
+ double dphht = phih-phit;
+ if(dphht>Math.PI) dphht-=2.*Math.PI;
+ if(dphht<-Math.PI) dphht+=2.*Math.PI;
+ double rphdis = dphht*rh;
double rtz = Math.sqrt(txz*txz+tyz*tyz);
- dists = ((txz-hpo[0])*(txz-hpo[0])+(tyz-hpo[1])*(tyz-hpo[1]))/c1+(rh-rtz)*(rh-rtz)/c2;
+// dists = rphdis*rphdis/c1+(rh-rtz)*(rh-rtz)/c2+(tzz-hpo[2])*(tzz-hpo[2])/0.0001;
+ dists = rphdis*rphdis/c1+(rh-rtz)*(rh-rtz)/c2;
}
// if(clrad) System.out.println("Normalized distance squared: "+df.format(dists));
if(dists < cldists)
{
- clht = th;
- active[i]=true;
+ clht = th;
+ clhin = hin-1;
+ active[i]=true;
cldists = dists;
// System.out.println("Found close hit! Z="+df.format(clht.getPosition()[2]));
}
}
+ }
if(active[i])
{
+ htassnd[clhin]=true;
nlrswh++;
hpo = clht.getPosition();
closeHits[i]=clht;
@@ -282,26 +333,36 @@
double phixr = Math.atan2(tyr,txr);
double dphx = phixr-phih;
double dpsi = dphx;
+ double dphd = Math.sqrt((txr-hpo[0])*(txr-hpo[0])+(tyr-hpo[1])*(tyr-hpo[1]));
+ dphd = Math.signum(dpsi)*dphd;
double rt = Math.sqrt(txr*txr+tyr*tyr);
double c1 = cvm[0]+exper1s;
double c2 = cvm[1]+exper2s;
- if(Math.abs(dphx) > Math.PI)
+ if(dphx > Math.PI) dphx-=2.*Math.PI;
+ if(dphx < -Math.PI) dphx+=2.*Math.PI;
+/* if(Math.abs(dphx) > Math.PI)
{
if(dpsi > 0.) dphx-=2.*Math.PI;
if(dpsi < 0.) dphx+=2.*Math.PI;
- }
+ } */
if(hitlrs[i] < 14)
{
- double norres = dphx*rt/Math.sqrt(c1);
-// double norres = dphx*rt;
+// 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) < 30.)
aida.cloud1D("WeightMatrix/ residual XY in lr "+hitlrs[i]).fill(norres);
- norres = (tzr-hpo[2])/Math.sqrt(c2);
-// norres = (tzr-hpo[2]);
+// norres = (tzr-hpo[2])/Math.sqrt(c2);
+ norres = (tzr-hpo[2]);
if(Math.abs(norres) < 30.)
- aida.cloud1D("WeightMatrix/ residual Z in lr "+hitlrs[i]).fill(norres);
- }
+ aida.cloud1D("WeightMatrix/ residual Z in lr "+hitlrs[i]).fill(norres);
+// if(Math.abs(expcorers) > 1.E-18)
+// {
+// norres = ((tzr-hpo[2])*(dphd));
+// norres = Math.signum(norres)*Math.sqrt(Math.abs(norres));
+// aida.cloud1D("WeightMatrix/ Corr. term Z and RPhi res. lr "+hitlrs[i]).fill(norres);
+// }
+ }
}
else
{
@@ -323,29 +384,58 @@
double phixr = Math.atan2(tyr,txr);
double dphx = phixr-phih;
double dpsi = dphx;
+ double dphd = Math.sqrt((txr-hpo[0])*(txr-hpo[0])+(tyr-hpo[1])*(tyr-hpo[1]));
+ dphd = Math.signum(dpsi)*dphd;
double rt = Math.sqrt(txr*txr+tyr*tyr);
double c1 = cvm[0]+exper1s;
double c2 = cvm[1]+exper2s;
- if(Math.abs(dphx) > Math.PI)
+ if(dphx > Math.PI) dphx-=2.*Math.PI;
+ if(dphx < -Math.PI) dphx+=2.*Math.PI;
+/* if(Math.abs(dphx) > Math.PI)
{
if(dpsi > 0.) dphx-=2.*Math.PI;
if(dpsi < 0.) dphx+=2.*Math.PI;
- }
- if((hitlrs[i] > 14)||((hitlrs[i]<9)&&(hitlrs[i]>4)))
+ } */
+ if((hitlrs[i] > 13)||((hitlrs[i]<9)&&(hitlrs[i]>4)))
{
-// aida.cloud1D("WeightMatrix/ Z position of disk hits").fill(hpo[2]);
- 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));
+// double norres = dphx*rt/Math.sqrt(c1);
+ double norres = dphx*rt;
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);
+// norres = (rh-rt)/Math.sqrt(c2);
+ norres = (rh-rt);
+ if((hitlrs[i])==5)
+ {
+ resr05 = rh-rt;
+ resp05 = dphd;
+ havh05 = true;
+ }
+ if(hitlrs[i]==14)
+ {
+ resr14 = rh-rt;
+ resp14 = dphd;
+ havh14 = true;
+ }
+ if(hitlrs[i]==23)
+ {
+ resr23 = rh-rt;
+ havh23 = true;
+ }
+ if(hitlrs[i] > 13)
+ {
+ aida.cloud1D("WeightMatrix/ Assigned hit Z poz for layer "+hitlrs[i]).fill(Math.abs(hpo[2]));
+ }
if(Math.abs(norres) < 30.)
aida.cloud1D("WeightMatrix/ residual R in lr "+hitlrs[i]).fill(norres);
- }
+// if(Math.abs(expcorers) > 1.E-18)
+// {
+// norres = ((rh-rt)*(dphd));
+// norres = Math.signum(norres)*Math.sqrt(Math.abs(norres));
+// aida.cloud1D("WeightMatrix/ Corr. term R and RPhi res. lr "+hitlrs[i]).fill(norres);
+// }
}
}
+ }
if(cvm[0] < 0.) System.out.println("Negative cov.matrix element 0: "+df.format(cvm[0]));
if(cvm[1] < 0.) System.out.println("Negative cov.matrix element 1: "+df.format(cvm[1]));
msm.set(i, i, cvm[0]);
@@ -355,20 +445,9 @@
msm.set(i+naclrs,i, 0.);
msm.set(i+naclrs, i+naclrs, cvm[1]);
}
+ }
else closeHits[i]=null;
}
-/* double midra = (maxra+minra)/2.;
- double drmax = maxra-minra;
- for(int i=0; i<naclrs; i++)
- {
- double normdr = 100. * (lrsxra[i]-midra)*(lrsxra[i]-midra)/(drmax*drmax);
- double orme = erm.get(i, i);
- erm.set(i,i, orme/(normdr+1.));
- } */
-// int m=msm.getRowDimension();
-// int n=msm.getColumnDimension();
-// System.out.println("rank "+mrank+" msm matrix dimensions "+m+" : "+n);
-
if(!ignoreMS) ermat = erm.plus(msm);
// if(!ignoreMS) ermat = erm;
else
@@ -403,6 +482,7 @@
nlrswh = 0;
for(int i=0; i<naclrs; i++)
if(active[i]) nlrswh++;
+ nmissing = naclrs - nlrswh;
int mranks = 2*nlrswh;
ermats = new Matrix(mranks,mranks);
int in=0;
@@ -422,9 +502,17 @@
if(active[j]) jn++;
}
if(active[i]) in++;
- }
-
- wm = ermats.inverse();
+ }
+ double corres = resr14*resr23;
+ corres = Math.signum(corres) * Math.sqrt(Math.abs(corres));
+ if(doHist && havh14 && havh23) aida.cloud1D("WeightMatrix/ Corr. btw lr14 and lr23 R res").fill(corres);
+ corres = resr05*resr23;
+ corres = Math.signum(corres) * Math.sqrt(Math.abs(corres));
+ if(doHist && havh05 && havh23) aida.cloud1D("WeightMatrix/ Corr. btw lr05 and lr23 R res").fill(corres);
+ corres = resp05*resr23;
+ corres = Math.signum(corres) * Math.sqrt(Math.abs(corres));
+ if(doHist && havh05 && havh23) aida.cloud1D("WeightMatrix/ Corr. btw lr05 phi and lr23 R res").fill(corres);
+ wm = ermats.inverse();
in = 0;
for(int i=0; i<naclrs; i++)
{
@@ -440,6 +528,23 @@
}
}
for(int i=nlrswh; i<MAXLRS; i++) active[i]=false;
+ calculateTrkInt();
+ if(debug_level == 1)
+ {
+ System.out.println("Sqrt of diagonal elements of error matrix: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(ermats.get(i,i))));
+ }
+ System.out.println(" ");
+ df.setMaximumFractionDigits(2);
+ System.out.println("Sqrt of diagonal elements of weight matrix: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(wm.get(i,i))));
+ }
+ System.out.println(" ");
+ }
boolean negel = false;
for(int i=0; i<mranks; i++)
{
@@ -470,6 +575,182 @@
}
// System.out.println("Weight matrix has assigned "+nlrswh+" closest hits to track");
}
+
+ public void simulateHits(MCParticle mcp)
+ {
+ if(debug_level==1) System.out.println("Simulating hits");
+ SmearTrackerHits smth = SmearTrackerHits.getInstance();
+ double[] smearXY = new double[4];
+ double[] smearZR = new double[4];
+ smearXY[0]=smth.getVtxBarrRPResolution();
+ smearZR[0]=smth.getVtxBarrZResolution();
+ smearXY[1]=smth.getVtxECRPResolution();
+ smearZR[1]=smth.getVtxECRResolution();
+ smearXY[2]=smth.getTrkBarrRPResolution();
+ smearZR[2]=smth.getTrkBarrZResolution();
+ smearXY[3]=smth.getTrkECRPResolution();
+ smearZR[3]=smth.getTrkECRResolution();
+ wm = null;
+ nlrswh = 0;
+ if(naclrs < 3) return;
+ int mrank=2*naclrs;
+ Matrix msm = new Matrix(mrank,mrank);
+ for(int i=0; i<mrank; i++)
+ for(int j=0; j<mrank; j++) msm.set(i,j,0.);
+ if(debug_level == 1) System.out.println("Subdetectors: ");
+ df.setMaximumFractionDigits(5);
+ for(int i=0; i<naclrs; i++)
+ {
+ boolean masked = false;
+ for(int j=0; j<nmasked; j++)
+ {
+ if(hitlrs[i] == masklrs[j]) masked = true;
+ }
+ active[i]=false;
+ closeHits[i] = null;
+ if(!masked)
+ {
+ Cylinder[i]=mser.isActiveCylinder(i);
+ double exper1s = erm.get(i,i);
+ double exper2s = erm.get(i+naclrs,i+naclrs);
+ double tx = trkpnts[i].x();
+ double ty = trkpnts[i].y();
+ double tz = trkpnts[i].z();
+ double[] pnt=new double[3];
+ pnt[0]=tx;
+ pnt[1]=ty;
+ pnt[2]=tz;
+ SimTrackerHit simth = new BaseSimTrackerHit(pnt, mcp);
+ List<SimTrackerHit> lst = new ArrayList();
+ lst.add(simth);
+ double hitrad = Math.sqrt(tx*tx+ty*ty);
+ if(doHist)
+ {
+ if(Cylinder[i])
+ {
+ aida.cloud1D("WeightMatrix/ Simulated hit Z poz for layer "+hitlrs[i]).fill(Math.abs(tz));
+ }
+ else
+ {
+ aida.cloud1D("WeightMatrix/ Simulated hit R for layer "+hitlrs[i]).fill(hitrad);
+ }
+ }
+ double[] cvm = new double[6];
+ SmearedTrackerHit smht = new SmearedTrackerHit(pnt,cvm,0.,0.,1,lst);
+ int dind = mser.getActLrsSubDet()[hitlrs[i]];
+// double sm1 = smearXY[dind]*ran.nextGaussian();
+// double sm2 = smearZR[dind]*ran.nextGaussian();
+ double sm1 = 0.;
+ double sm2 = 0.;
+ if(debug_level == 1) System.out.print(" "+dind);
+ double[] V = new double[3];
+ double[] tcm;
+ V[0] = smearXY[dind]*smearXY[dind];
+ V[1] = smearZR[dind]*smearZR[dind];
+ V[2] = 0.;
+ if((dind==0)||(dind==2))
+ {
+ smht.smear(sm1,sm2,0.);
+ tcm = cmt.cyl2Cart(V,smht.getPoint());
+ smht.setCovMatrix(tcm);
+ }
+ if((dind==1)||(dind==3))
+ {
+ smht.smear(sm1, 0., sm2);
+ tcm = cmt.disk2Cart(V,smht.getPoint());
+ smht.setCovMatrix(tcm);
+ }
+ closeHits[i]=smht;
+ trkPnts[i]=trkpnts[i];
+ active[i]=true;
+ msm.set(i, i, V[0]);
+ msm.set(i+naclrs, i+naclrs, V[1]);
+ }
+ }
+ ermat = erm.plus(msm);
+ nlrswh = 0;
+ for(int i=0; i<naclrs; i++)
+ if(active[i]) nlrswh++;
+ nmissing = naclrs - 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;
+ calculateTrkInt();
+ if(debug_level == 1)
+ {
+ int nprnt = 12;
+ if(nlrswh < nprnt) nprnt = nlrswh;
+ System.out.println("Error matrix (first 12 elements):");
+ for(int i=0; i<nprnt; i++)
+ {
+ for(int j=0; j<nprnt; j++)
+ {
+ System.out.print(dfe.format(ermats.get(i,j))+" ");
+ }
+ System.out.println(" ");
+ }
+ System.out.println(" ");
+ System.out.println("Sqrt of diagonal elements of error matrix for first coord: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(ermats.get(i,i))));
+ }
+ System.out.println(" ");
+ System.out.println("Sqrt of diagonal elements of error matrix for second coord: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(ermats.get(i+nlrswh,i+nlrswh))));
+ }
+ df.setMaximumFractionDigits(2);
+ System.out.println("Sqrt of diagonal elements of weight matrix for first coord: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(wm.get(i,i))));
+ }
+ System.out.println(" ");
+ System.out.println("Sqrt of diagonal elements of weight matrix for second coord: ");
+ for(int i=0; i<nlrswh; i++)
+ {
+ System.out.print(" "+df.format(Math.sqrt(wm.get(i+nlrswh,i+nlrswh))));
+ }
+ System.out.println(" ");
+ }
+ }
+
+
public void calculateTrkInt()
{
for(int i=0; i< nlrswh; i++)
CVSspam 0.2.8
