6 modified files
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFConstants.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFConstants.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -3,20 +3,22 @@
public class KFConstants
{
- public KFConstants()
- {
- }
+ public KFConstants()
+ {
+ }
- public double Pi() {return pi;}
- public double twoPi() {return twopi;}
- public double vtxResolution() {return vtxres;}
- public double sodResolution() {return sodres;}
- public double pointThreeTimesB() {return ptconst;}
+ public double Pi() {return pi;}
+ public double twoPi() {return twopi;}
+ public double vtxResolution() {return vtxres;}
+ public double sodResolution() {return sodres;}
+ public double pointThreeTimesB() {return ptconst;}
+ public double pionMass() {return pimass;}
- private static double pi = 3.14159265359;
- private static double twopi = 6.283185307;
- private static double vtxres = 0.0005;
- private static double sodres = 0.0007;
- private static double ptconst = 0.015; // = 0.3*B
+ private static double pi = 3.14159265359;
+ private static double twopi = 6.283185307;
+ private static double vtxres = 0.0005;
+ private static double sodres = 0.0007;
+ private static double ptconst = 0.015; // = 0.3*B
+ private static double pimass = 0.13957018;
}
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFFitterDriver.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFFitterDriver.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -28,7 +28,10 @@
{
super.process(event); // this takes care that the child Drivers are loaded and processed.
+ // 1. Extract list of SODTracks from event
List<SODTrack> SODTrackList = (List<SODTrack>) event.get("SODTracks");
+
+ // 2. Loop over SODTracks and fit them
for ( SODTrack SODtrk : SODTrackList ) {
double trkPT = SODtrk.getPT();
@@ -41,6 +44,7 @@
System.out.print(SODtrk.getTrackParameter(4)); System.out.print(" ");
System.out.println(SODtrk.getCharge());
+ // create new KFTrack and initialize it to the SODTrack fit parameters
KFTrack track = new KFTrack();
track.defineSeedTrack(SODtrk);
@@ -49,22 +53,21 @@
System.out.print("Driver: SOD hit list size = "); System.out.println(nbHits);
double doca = 0.1*SODtrk.getTrackParameter(0);
- System.out.print("Driver: SOD doca = "); System.out.print(nbHits); System.out.println(" cm");
+ System.out.print("Driver: SOD doca = "); System.out.print(doca); System.out.println(" cm");
double chi2 = SODtrk.getChi2();
- System.out.print("Driver: SOD chi2 = "); System.out.print(chi2);
+ System.out.print("Driver: SOD chi2 = "); System.out.println(chi2);
aida.cloud1D("chi2").fill(chi2);
double ndf = nbHits;
if(ndf>0.) { double r = chi2/ndf; aida.cloud1D("chi2 per ndf").fill(r); }
- //if(Math.abs(doca)<0.01) {
- if(chi2<100.) {
+ if(chi2<100.) { // run fitter only if good track seed
for ( SimTrackerHit hit : hits ) {
- track.addHit(hit);
+ track.addHit(hit); // assign hits to KFTrack
}
System.out.print("Driver: track params before fit = ");
track.getTrackParameters().printFitParams();
- track.fit();
+ track.fit(); // run Kalman Filter Fit
System.out.println();
double KFpt = track.getFittedPtAtPOCA();
@@ -87,7 +90,7 @@
} else {
System.out.println("Driver: skip this track");
System.out.println();
- }
- }
+ } // if(chi2<100.) {
+ } // for ( SODTrack SODtrk : SODTrackList ) {
}
}
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFPoint.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFPoint.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -29,7 +29,10 @@
_phi = Math.atan2(_y,_x);
_rphiErr = 0.0007; if(_r<10.){_rphiErr = 0.0005;}
- _zErr = 0.005; //if(_r>10.){_zErr = 1000.;}
+ _zErr = 0.0005;
+ if(_r>10.){_zErr = 1000.;}
+
+ _zErr = 1e+6; // fix me
//System.out.print("KFPoint: phi, rphiErr = ");
//System.out.print(_phi); System.out.print(" ");
@@ -42,6 +45,7 @@
public double getZ() { return _z; }
public double getR() { return _r; }
public double getPhi() { return _phi; }
+ public double getDistanceFromOrigin() {return Math.sqrt(_r*_r+_z*_z);}
public double getRPhiErr() { return _rphiErr;}
public double getZErr() { return _zErr;}
@@ -49,7 +53,6 @@
public Matrix getPhiZErrorMatrix()
{
double[][] vals = {{_rphiErr*_rphiErr/(_r*_r),0.},{0.,_zErr*_zErr}};
- //double[][] vals = {{_rphiErr*_rphiErr,0.},{0.,_zErr*_zErr}};
Matrix errMat = new Matrix(vals);
return errMat;
}
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFSite.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFSite.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -13,9 +13,6 @@
public KFSite()
{
- _siteType = 1; // 1=hit
- //_siteType = 2; // 2=material
- //_siteType = 3; // 3=B field discontinuity
}
public void setLocation(SimTrackerHit hit)
@@ -38,115 +35,302 @@
//double[][] vals = {{rphierr,0.},{0.,zerr}};
_PhiZErrorMatrix = _siteLocation.getPhiZErrorMatrix();
- System.out.println("KFSite: PhiZ cov matrix:");
- _PhiZErrorMatrix.print(12,10);
-
_PhiZErrorMatrixInverse = _PhiZErrorMatrix.inverse();
- System.out.println("KFSite: PhiZ cov matrix inverted:");
- _PhiZErrorMatrixInverse.print(12,10);
+ System.out.println("KFSite: PhiZ cov matrix:");
+ _PhiZErrorMatrix.print(12,10);
- }
+ //System.out.println("KFSite: PhiZ cov matrix inverted:");
+ //_PhiZErrorMatrixInverse.print(12,10);
- public void setTrackParOut(KFTrackParameters trk)
- {
- _trackParamOutward = trk;
- }
-
- public void setTrackParIn(KFTrackParameters trk)
- {
- _trackParamInward = trk;
}
public KFTrackParameters filter(KFTrackParameters trk)
{
- // define matrix relating the measured parameters (phi,z) at radius R and the
- // state vector (fit parameters). In the current setup, the state vector is (phi,z,px,py,pz)
- double[][] vals = {{1.,0.,0.,0.,0.},{0.,1.,0.,0.,0.}};
- Matrix H = new Matrix(vals);
-
- Matrix HT = H.transpose();
- System.out.println("KFSite: H:");
- H.print(12,10);
- System.out.println("KFSite: HT:");
- HT.print(12,10);
+ double[] tpp = trk.getFitParameters();
+ Matrix stateVect = new Matrix(tpp,5);
+ System.out.println("KFSite: initial stateVect:");
+ stateVect.print(15,10);
- Matrix trkCov = trk.getCovMatrixAtRadiusR();
+ double d0 = tpp[0];
+ double phi0 = tpp[1];
+ double omega = tpp[2];
+ double z0 = tpp[3];
+ double tanl = tpp[4];
+ double l = 0.;
+
+ Matrix trkCovNoNoise = trk.getErrorMatrix();
+ System.out.println("KFSite: trkCovNoNoise:");
+ trkCovNoNoise.print(15,12);
+
+ Matrix Q = new Matrix(5,5,0.);
+
+ if(_materialSite) {
+ //add here process noise due to dE/dx and scattering effects
+
+ // Consider 3 unit vectors:
+ // e1 = (px,py,pz)/p
+ // e2 = (-py,px,0)/pt
+ // e3 = (-px*pz,-py*pz,px*px+py*py)/(p*pt)
+ // Original direction = e1
+ // Scattered direction =
+ // sqrt(1-sin^2(t2)*sin^2(t3))*e1 + sin(t2)*e2 + sin(t3)*e3
+ //
+ // 1. Determine d(d0,phi0,omega,z0,s)/d(theta1) and d(d0,phi0,omega,z0,s)/d(theta2)
+ // Create tracks with momentum rotated by small amounts about e2 and e3
+ // Measure difference in fit parameters due to these variations
+ // --> get 5x2 matrix SCT = dp/dtheta
+ // 2. Determine average scattering angle theta -> sigmaSCT
+ // --> 2x2 covariance matrix sigSCT = ( sigmaSCT^2 0 )
+ // ( 0 sigmaSCT^2 )
+ // 3. Q_scattering = SCT x sigSCT x SCT^T
+
+ // 1. Determine derivatives
+ double sctVar = 0.000001; // variation of scattering angle
+ double[] ptPar = trk.getPointParameters();
+ double x = ptPar[0];
+ double y = ptPar[1];
+ double z = ptPar[2];
+ double px = ptPar[3];
+ double py = ptPar[4];
+ double pz = ptPar[5];
+ int chg = trk.getCharge();
+
+ double pt = Math.sqrt(px*px+py*py);
+ double p = Math.sqrt(px*px+py*py+pz*pz);
+
+ double pxt1p = Math.cos(sctVar)*px - Math.sin(sctVar)*py*p/pt;
+ double pyt1p = Math.cos(sctVar)*py + Math.sin(sctVar)*px*p/pt;
+ double pzt1p = Math.cos(sctVar)*pz;
+
+ double pxt1m = Math.cos(sctVar)*px + Math.sin(sctVar)*py*p/pt;
+ double pyt1m = Math.cos(sctVar)*py - Math.sin(sctVar)*px*p/pt;
+ double pzt1m = Math.cos(sctVar)*pz;
+
+ double pxt2p = Math.cos(sctVar)*px - Math.sin(sctVar)*px*pz/pt;
+ double pyt2p = Math.cos(sctVar)*py - Math.sin(sctVar)*py*pz/pt;
+ double pzt2p = Math.cos(sctVar)*pz + Math.sin(sctVar)*pt;
+
+ double pxt2m = Math.cos(sctVar)*px + Math.sin(sctVar)*px*pz/pt;
+ double pyt2m = Math.cos(sctVar)*py + Math.sin(sctVar)*py*pz/pt;
+ double pzt2m = Math.cos(sctVar)*pz - Math.sin(sctVar)*pt;
+
+
+ KFTrackParameters trkVar1P = new KFTrackParameters();
+ trkVar1P.setTrackPointParameters(chg,x,y,z,pxt1p,pyt1p,pzt1p);
+ double[] fit1P = trkVar1P.getFitParameters();
+ KFTrackParameters trkVar1M = new KFTrackParameters();
+ trkVar1M.setTrackPointParameters(chg,x,y,z,pxt1m,pyt1m,pzt1m);
+ double[] fit1M = trkVar1M.getFitParameters();
+ double[] v1 = new double [5];
+ for(int i1=0;i1<5;i1++) {v1[i1]=(fit1P[i1]-fit1M[i1])/(2.*sctVar);}
+
+ KFTrackParameters trkVar2P = new KFTrackParameters();
+ trkVar2P.setTrackPointParameters(chg,x,y,z,pxt2p,pyt2p,pzt2p);
+ double[] fit2P = trkVar2P.getFitParameters();
+ KFTrackParameters trkVar2M = new KFTrackParameters();
+ trkVar2M.setTrackPointParameters(chg,x,y,z,pxt2m,pyt2m,pzt2m);
+ double[] fit2M = trkVar2M.getFitParameters();
+ double[] v2 = new double [5];
+ for(int i2=0;i2<5;i2++) {v2[i2]=(fit2P[i2]-fit2M[i2])/(2.*sctVar);}
+
+ double[][] vals = {{v1[0],v1[1],v1[2],v1[3],v1[4]},{v2[0],v2[1],v2[2],v2[3],v2[4]}};
+ Matrix SCTderiv = new Matrix(vals);
+ Matrix SCTderivT = SCTderiv.transpose();
+ System.out.println("KFSite: SCTderiv:");
+ SCTderiv.print(15,12);
+
+ // 2. Average scattering angle
+ double mass = 0.13957; // assume pion hypothesis
+ double beta = p/Math.sqrt(p*p+mass*mass); // use particle momentum
+ double sctRMS = (0.0136/(beta*p))*Math.sqrt(_radThickness)*(1+0.038*Math.log(_radThickness));
+
+ double[][] valsSCT = {{sctRMS*sctRMS,0.},{0.,sctRMS*sctRMS}};
+ Matrix sigSCT = new Matrix(valsSCT);
+ System.out.println("KFSite: sigSCT:");
+ sigSCT.print(15,12);
+
+ // 3. determine Q matrix for scattering
+ Q = SCTderivT.times(sigSCT).times(SCTderiv);
+ System.out.println("KFSite: Q:");
+ Q.print(15,12);
+ }
+ Matrix trkCov = trkCovNoNoise.plus(Q);
System.out.println("KFSite: trkCov:");
- trkCov.print(15,4);
-
- Matrix trkCovInverse = trkCov.inverse();
-
- System.out.println("KFSite: trkCovInverse:");
- trkCovInverse.print(15,4);
+ trkCov.print(15,12);
- // Check matrix inversion
- Matrix checkInv = trkCov.times(trkCovInverse);
- System.out.println("KFSite: trkCov*trkCovInverse:");
- checkInv.print(15,4);
+ Matrix filteredCov = new Matrix(5,5);
- double[] tpp = trk.getPFitParameters();
- Matrix stateVect = new Matrix(tpp,5);
-
- System.out.println("KFSite: stateVect:");
- stateVect.print(12,10);
+ if(_measurementSite) {
+ // define matrix relating the measured parameters (phi,z) at radius R and the
+ // state vector (fit parameters). In the current setup, the state vector is
+ // given by the fit parameters (d0,phi0,omega,z0,s)
+
+ double[] vP = new double [5];
+ double[] vZ = new double [5];
+
+ double[] variations = new double[] {1e-05,1e-05,1e-10,1e-05,1e-06};
+
+ double[] fitPars = trk.getFitParameters();
+ double pthlength = 0.;
+
+ for(int j=0; j<5; j++) {
+ double p1p = fitPars[0]; double p1m = fitPars[0];
+ if(j==0) { p1p += 0.5*variations[0]; p1m -= 0.5*variations[0]; }
+ double p2p = fitPars[1]; double p2m = fitPars[1];
+ if(j==1) { p2p += 0.5*variations[1]; p2m -= 0.5*variations[1]; }
+ double p3p = fitPars[2]; double p3m = fitPars[2];
+ if(j==2) { p3p += 0.5*variations[2]; p3m -= 0.5*variations[2]; }
+ double p4p = fitPars[3]; double p4m = fitPars[3];
+ if(j==3) { p4p += 0.5*variations[3]; p4m -= 0.5*variations[3]; }
+ double p5p = fitPars[4]; double p5m = fitPars[4];
+ if(j==4) { p5p += 0.5*variations[4]; p5m -= 0.5*variations[4]; }
+
+ KFTrackParameters trkVarPlus = new KFTrackParameters();
+ trkVarPlus.setTrackFitParameters(p1p, p2p, p3p, p4p, p5p,pthlength);
+ KFTrackParameters trkVarPlusPropagated =
+ trkVarPlus.makePredictionAtRadius(_siteLocation.getR());
+ KFTrackParameters trkVarMinus = new KFTrackParameters();
+ trkVarMinus.setTrackFitParameters(p1m, p2m, p3m, p4m, p5m,pthlength);
+ KFTrackParameters trkVarMinusPropagated =
+ trkVarMinus.makePredictionAtRadius(_siteLocation.getR());
+
+ vP[j] = (trkVarPlusPropagated.getPhi()-trkVarMinusPropagated.getPhi())/variations[j];
+ vZ[j] = (trkVarPlusPropagated.getZ()-trkVarMinusPropagated.getZ())/variations[j];
+
+ //System.out.print("KFSite: vP = ");
+ //System.out.print(vP[j]);
+ //System.out.print(" ; vZ = ");
+ //System.out.println(vZ[j]);
+ }
+
+ double[][] vals = {{vP[0],vP[1],vP[2],vP[3],vP[4]},{vZ[0],vZ[1],vZ[2],vZ[3],vZ[4]}};
+ Matrix H = new Matrix(vals);
+
+ Matrix HT = H.transpose();
+ System.out.println("KFSite: H:");
+ H.print(12,10);
+ System.out.println("KFSite: HT:");
+ HT.print(12,10);
+
+ Matrix trkCovInverse = trkCov.inverse();
+
+ System.out.println("KFSite: trkCovInverse:");
+ trkCovInverse.print(15,12);
+
+ double[] mp = {_siteLocation.getPhi(),_siteLocation.getZ()};
+ Matrix measVect = new Matrix(mp,2);
+
+ System.out.println("KFSite: measVect:");
+ measVect.print(12,10);
+
+ System.out.println("KFSite: PhiZ cov matrix inverted:");
+ _PhiZErrorMatrixInverse.print(12,10);
+
+ Matrix cv1 = _PhiZErrorMatrixInverse.times(H);
+ Matrix cv2 = HT.times(cv1);
+ System.out.println("KFSite: measCovInverse:");
+ cv2.print(15,12);
+ Matrix filteredCovInv = trkCovInverse.plus(cv2);
+ filteredCov = filteredCovInv.inverse();
+
+ System.out.println("KFSite: Filtered Cov:");
+ filteredCov.print(15,12);
+
+ Matrix tmp1 = trkCovInverse.times(stateVect);
+ Matrix tmp2 = HT.times(_PhiZErrorMatrixInverse).times(measVect);
+ Matrix tmp3 = tmp1.plus(tmp2);
+ Matrix filteredTrkPar = filteredCov.times(tmp3);
+
+ System.out.println("KFSite: Filtered State Vector:");
+ filteredTrkPar.print(6,6);
+
+ double[][] par = filteredTrkPar.getArrayCopy();
+ d0 = par[0][0];
+ phi0 = par[1][0];
+ omega = par[2][0];
+ z0 = par[3][0];
+ tanl = par[4][0];
+
+ System.out.print("KFSite: Fitted Pt = "); System.out.println(0.015/Math.abs(omega));
+
+ } //if(_measurementSite) {
+
+ KFTrackParameters tmpTrk = new KFTrackParameters();
+ tmpTrk.setTrackFitParameters(d0,phi0,omega,z0,tanl,l);
+ tmpTrk.setPointToRadius(_siteLocation.getR());
+ double tmpPtot = tmpTrk.getPtot();
+ double dEdx = dEdx(tmpPtot);
+ System.out.print("KFSite: Ptot = "); System.out.println(tmpPtot);
+ System.out.print("KFSite: dEdx = "); System.out.println(dEdx);
+
+ double tmpEtot = Math.sqrt(_cst.pionMass()*_cst.pionMass()+tmpPtot*tmpPtot);
+ double tmpEtotCorr = tmpEtot+0.01*dEdx; // (fix me) change sign depending on fit direction
+ double ratioCorr = tmpEtotCorr/tmpEtot;
- //double[] mp = {_siteLocation.getPhi()*_siteLocation.getR(),_siteLocation.getZ()};
- double[] mp = {_siteLocation.getPhi(),_siteLocation.getZ()};
- //double[] mp = {_siteLocation.getPhi(),0.}; // fix me
- Matrix measVect = new Matrix(mp,2);
-
- System.out.println("KFSite: measVect:");
- measVect.print(12,10);
-
- Matrix cv1 = _PhiZErrorMatrixInverse.times(H);
- Matrix cv2 = HT.times(cv1);
- Matrix filteredCovInv = trkCovInverse.plus(cv2);
- Matrix filteredCov = filteredCovInv.inverse();
-
- System.out.println("KFSite: Filtered Cov:");
- filteredCov.print(15,12);
-
- Matrix tmp1 = trkCovInverse.times(stateVect);
- Matrix tmp2 = HT.times(_PhiZErrorMatrixInverse).times(measVect);
- Matrix tmp3 = tmp1.plus(tmp2);
- Matrix filteredTrkPar = filteredCov.times(tmp3);
-
- System.out.println("KFSite: Filtered State Vector:");
- filteredTrkPar.print(6,6);
-
- double[][] par = filteredTrkPar.getArrayCopy();
- int charge = trk.getCharge();
- double phi = par[0][0];
- double x1 = _siteLocation.getR()*Math.cos(phi);
- double x2 = _siteLocation.getR()*Math.sin(phi);
- double x3 = par[1][0];
- double p1 = par[2][0];
- double p2 = par[3][0];
- double p3 = par[4][0];
-
- System.out.print("KFSite: Fitted Pt = "); System.out.println(Math.sqrt(p1*p1+p2*p2));
-
- //double p1 = par[0][0];
- //double p2 = par[1][0];
- //double p3 = par[2][0];
- //double p4 = par[3][0];
- //double p5 = par[4][0];
+ double[] tmpPar = tmpTrk.getPointParameters();
+ int tmpCh = tmpTrk.getCharge();
KFTrackParameters filtTrk = new KFTrackParameters();
- //filtTrk.setTrackFitParameters(charge,p1,p2,p3,p4,p5);
- //filtTrk.setPointToRadius(_siteLocation.getR());
- filtTrk.setTrackPointParameters(charge,x1,x2,x3,p1,p2,p3);
- filtTrk.setErrorMatrixAtRadiusR(filteredCov);
+ filtTrk.setTrackPointParameters(tmpCh,
+ tmpPar[0],
+ tmpPar[1],
+ tmpPar[2],
+ ratioCorr*tmpPar[3],
+ ratioCorr*tmpPar[4],
+ ratioCorr*tmpPar[5]);
+
+ //KFTrackParameters filtTrk = new KFTrackParameters();
+ //filtTrk.setTrackFitParameters(d0,phi0,omega,z0,tanl,l);
+
+ if(_measurementSite) {
+ filtTrk.setErrorMatrix(filteredCov);
+ } else {
+ filtTrk.setErrorMatrix(trkCov);
+ }
+ filtTrk.setPointToRadius(_siteLocation.getR());
return filtTrk;
}
+ public void setTrackParOut(KFTrackParameters trk) {_trackParamOutward = trk;}
+ public void setTrackParIn(KFTrackParameters trk) {_trackParamInward = trk;}
+
public KFPoint getSiteLocation() {return _siteLocation;}
+ public void isMeasurementSite(boolean bool) {_measurementSite=bool;}
+ public void isMaterialSite(boolean bool,double thickness) {
+ _materialSite=bool;
+ _radThickness = thickness;
+ }
+ public void isBFieldDiscontinuitySite(boolean bool) {_BFieldDiscontinuitySite=bool;}
+
+ private double dEdx(double p)
+ {
+ double dedx = 0.0;
+ double rho = 2.33; // g/cm^3 (Si)
+ double Z = 14.; // Si
+ double A = 28.0855; // Si
+ double beta = p/Math.sqrt(p*p+_cst.pionMass()*_cst.pionMass());
+ double gamma = 1./Math.sqrt(1-beta*beta);
+ double eta = beta*gamma;
+ double s = 0.510998/_cst.pionMass();
+ double Wmax = 2.*0.510998*eta*eta/(1.+2.*s*Math.sqrt(1+eta*eta)+s*s);
+ double I = 172.;
+
+ double frac = 2.*0.510998*eta*eta*Wmax/(I*I);
+
+ double delta = 0.;
+ double C = 0.;
+
+ dedx = -0.1535*rho*(Z/A)*(Math.log(frac)-2*beta*beta-delta-2*C/Z)/(beta*beta);
+
+ return dedx;
+ }
+
private KFTrackParameters _trackParamOutward;
private KFTrackParameters _trackParamOutwardCombined;
private KFTrackParameters _trackParamInward;
@@ -154,9 +338,14 @@
private KFPoint _siteLocation = new KFPoint();
- private Matrix _PhiZErrorMatrix = new Matrix(3,3);
- private Matrix _PhiZErrorMatrixInverse = new Matrix(3,3);
+ private Matrix _PhiZErrorMatrix = new Matrix(2,2);
+ private Matrix _PhiZErrorMatrixInverse = new Matrix(2,2);
+
+ private boolean _measurementSite = false;
+ private boolean _materialSite = false;
+ private boolean _BFieldDiscontinuitySite = false;
- private int _siteType;
+ private double _radThickness = 0.;
+ private KFConstants _cst = new KFConstants();
}
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFTrack.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFTrack.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -16,6 +16,7 @@
public KFTrack()
{
_siteList = new java.util.LinkedList();
+ _orderedSiteList = new java.util.LinkedList();
}
public void addHit(SimTrackerHit hit)
@@ -23,25 +24,55 @@
// define site for this hit
KFSite kfhit = new KFSite();
kfhit.setLocation(hit);
+ kfhit.isMeasurementSite(true);
+ double radius = kfhit.getSiteLocation().getR();
+ double radThickness = 0.00107;
+ if (radius>10.) {radThickness=0.00532;}
+ kfhit.isMaterialSite(true,radThickness);
_siteList.add(kfhit);
//System.out.print("KFTrack: read back radius = ");
//System.out.println(kfhit.getSiteLocation().getR());
}
+ private void sortSiteList()
+ {
+ // loop over _siteList
+ ListIterator hits = _siteList.listIterator();
+ int hitCounter = 0;
+ double[] dist = new double[100];
+ for (; hits.hasNext();) {
+ KFSite hit = (KFSite)hits.next();
+ // record ordering parameter (e.g. distance from origin)
+ dist[hitCounter] = hit.getSiteLocation().getDistanceFromOrigin();
+ hitCounter++;
+ }
+ // loop over _siteList and fill _orderedSiteList
+ double currentDist = 0.;
+ int index = 0;
+ for (int j=0;j<hitCounter;j++) {
+ double nextSmallestDist = 100000.;
+ for (int i=0;i<hitCounter;i++) {
+ if(dist[i]<nextSmallestDist&&dist[i]>currentDist) {
+ index=i;
+ nextSmallestDist=dist[i];
+ }
+ }
+ _orderedSiteList.add(_siteList.get(index));
+ currentDist = nextSmallestDist;
+ }
+ }
+
public void defineSeedTrack(SODTrack seedTrack)
{
// load track parameters d0,omega,...
double d = seedTrack.getTrackParameter(0)/10.; // divided by 10 to get cm
double phi = seedTrack.getTrackParameter(1);
double omega = seedTrack.getTrackParameter(2)*(-10.); // mult. by -10 to get cm-1
- //omega = 2.*omega; // fix me
double z = seedTrack.getTrackParameter(3)/10.; // divided by 10 to get cm
double s = seedTrack.getTrackParameter(4);
- //double z = 0.; // fix me
- //double s = 0.; // fix me
- int charge = 1; if(omega>0.) {charge = -1;}
+ double l = 0.0; // path length set to zero (=POCA)
System.out.print("KFTrack: define track par = ");
System.out.print(d); System.out.print(" ");
@@ -50,21 +81,25 @@
System.out.print(z); System.out.print(" ");
System.out.println(s);
- _seedTrack.setTrackFitParameters(charge,d,phi,omega,z,s);
+ _seedTrack.setTrackFitParameters(d,phi,omega,z,s,l);
}
public void fit()
{
+ // sort list of hits by distance from origin
+ sortSiteList();
+
// initialize track param covariance matrix
double[][] diag = new double[5][5];
- diag[0][0] = 0.1;
- diag[1][1] = 10.;
- for(int i=2; i<5; i++) {diag[i][i] = 100.;}
+ diag[0][0] = 10.; // d0
+ diag[1][1] = 1.; // phi0
+ diag[2][2] = 100.; // omega
+ diag[3][3] = 10.; // z0
+ diag[4][4] = 10.; // tanl
Matrix initCov = new Matrix(diag);
System.out.println("KFTrack: init cov matrix");
initCov.print(6,2);
- //_seedTrack.setErrorMatrixAtRadiusR(initCov);
System.out.print("KFTrack: ");
_seedTrack.printFitParams();
@@ -75,8 +110,10 @@
//_seedTrack.setPointToPOCA();
double refRadius = 1.0;
currentTrackOut.setPointToRadius(refRadius);
+ System.out.print("KFTrack: set radius to 1 : ");
+ currentTrackOut.printFitParams();
// initialize error matrix
- currentTrackOut.setErrorMatrixAtRadiusR(initCov);
+ currentTrackOut.setErrorMatrix(initCov);
// initialize track parameters for inward fit
@@ -84,10 +121,12 @@
// determine track params at fixed small radius
refRadius = 130.0;
currentTrackIn.setPointToRadius(refRadius);
+ System.out.print("KFTrack: set radius to 130 : ");
+ currentTrackIn.printFitParams();
// initialize error matrix
- currentTrackIn.setErrorMatrixAtRadiusR(initCov);
+ currentTrackIn.setErrorMatrix(initCov);
- ListIterator hits = _siteList.listIterator();
+ ListIterator hits = _orderedSiteList.listIterator();
//***************
// Fit inside-out
@@ -97,25 +136,27 @@
for (; hits.hasNext();){
KFSite hit = (KFSite)hits.next();
- KFPoint siteLocation = new KFPoint();
- siteLocation = hit.getSiteLocation();
- double R = siteLocation.getR();
+ //KFPoint siteLocation = new KFPoint();
+ //siteLocation = hit.getSiteLocation();
+ //double R = siteLocation.getR();
+ double R = hit.getSiteLocation().getR();
System.out.print("KFTrack: hit radius = "); System.out.println(R);
System.out.print("KFTrack: current pars = ");
- currentTrackOut.printPFitParams();
+ currentTrackOut.printFitParams();
// propagate track to radius of hit
KFTrackParameters predictedTrack = currentTrackOut.makePredictionAtRadius(R);
System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printPFitParams();
+ predictedTrack.printFitParams();
+ predictedTrack.printPointParams();
// filter predicted track parameters with hit at radius R
currentTrackOut = hit.filter(predictedTrack);
- // record track parameters in KFSite (fix me)
-
+ // record track parameters in KFSite
+ hit.setTrackParOut(currentTrackOut);
}
//***************
@@ -126,24 +167,26 @@
for (; hits.hasPrevious();){
KFSite hit = (KFSite)hits.previous();
- KFPoint siteLocation = new KFPoint();
- siteLocation = hit.getSiteLocation();
- double R = siteLocation.getR();
+ //KFPoint siteLocation = new KFPoint();
+ //siteLocation = hit.getSiteLocation();
+ //double R = siteLocation.getR();
+ double R = hit.getSiteLocation().getR();
System.out.print("KFTrack: hit radius = "); System.out.println(R);
System.out.print("KFTrack: current pars = ");
- currentTrackIn.printPFitParams();
+ currentTrackIn.printFitParams();
// propagate track to radius of hit
KFTrackParameters predictedTrack = currentTrackIn.makePredictionAtRadius(R);
System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printPFitParams();
+ predictedTrack.printFitParams();
// filter predicted track parameters with hit at radius R
currentTrackIn = hit.filter(predictedTrack);
- // record track parameters in KFSite (fix me)
+ // record track parameters in KFSite
+ hit.setTrackParIn(currentTrackIn);
}
@@ -156,6 +199,7 @@
public double getFittedPtAtPOCA() { return _fittedPtAtPOCA; }
private java.util.LinkedList _siteList;
+ private java.util.LinkedList _orderedSiteList;
private KFTrackParameters _seedTrack = new KFTrackParameters();
private double _fittedPtAtPOCA;
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFTrackParameters.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFTrackParameters.java 30 Jan 2007 17:01:58 -0000 1.2
@@ -13,43 +13,30 @@
KFTrackParameters(KFTrackParameters trk)
{
double[] par = trk.getFitParameters();
- int chg = trk.getCharge();
- this.setTrackFitParameters(chg, par[0], par[1], par[2], par[3], par[4]);
+ double l = trk.getPathLength();
+ this.setTrackFitParameters(par[0], par[1], par[2], par[3], par[4], l);
- double[] ptp = trk.getPointParameters();
- for(int i=0; i<6; i++) { this._pointParams[i] = ptp[i]; }
- _x = ptp[0]; _y = ptp[1]; _z = ptp[2]; _px = ptp[3]; _py = ptp[4]; _pz = ptp[5];
- _r = Math.sqrt(_x*_x+_y*_y);
- //_phi = 0.; if(_x!=0.&&_y!=0) { _phi = _r*Math.atan2(_y,_x); }
- _phi = 0.; if(_x!=0.&&_y!=0) { _phi = Math.atan2(_y,_x); }
- this._pFitParams[0] = _phi;
- for(int i=2; i<6; i++) { this._pFitParams[i-1] = ptp[i]; }
-
- this._FitErrorMatrix = trk.getFitCovMatrix();
- this._ErrorMatrixAtRadiusR = trk.getCovMatrixAtRadiusR();
- this._ErrorMatrixAtRadiusRprime = trk.getCovMatrixAtRadiusRprime();
+ this._FitErrorMatrix = trk.getErrorMatrix();
}
private int _charge = 0;
private double _d0 = 0., _phi0 = 0., _omega = 0., _z0 = 0., _s = 0.;
+ private double _l = 0., _alpha = 0.;
private double[] _fitParams = new double[] {0.,0.,0.,0.,0.};
- private double _p = 0., _pt = 0., _px0 = 0., _py0 = 0., _pz0 = 0.;
- private double _x0 = 0., _y0 = 0.;
- private double[] _pointOriginParams = new double[] {0.,0.,0.,0.,0.,0.};
-
- private double _xc = 0., _yc = 0., _rc = 0.;
+ private double _xc = 0., _yc = 0., _r0 = 0., _rc = 0, _pt = 0.;
private double[] _helixParams = new double[] {0.,0.,0.,0.,0.};
- private double _px = 0., _py = 0., _pz = 0.;
- private double _x = 0., _y = 0., _z = 0., _r, _phi = 0.;
+ private double _px = 0., _py = 0., _pz = 0., _p = 0.;
+ private double _x = 0., _y = 0., _z = 0.;
private double[] _pointParams = new double[] {0.,0.,0.,0.,0.,0.};
- private double[] _pFitParams = new double[] {0.,0.,0.,0.,0.};
- public void setTrackFitParameters(int ch, double dzero, double phizero, double omeg, double zzero, double spar)
+ public void setTrackFitParameters(double dzero, double phizero, double omeg,
+ double zzero, double spar, double pathLength)
{
- _charge = ch;
+ _charge = (int) (-1*Math.signum(omeg));
+ //System.out.print("charge = "); System.out.println(_charge);
_d0 = dzero;
_phi0 = phizero;
_omega = omeg; if(_omega==0.) {System.out.println("Null curvature!");}
@@ -62,11 +49,18 @@
_fitParams[3] = _z0;
_fitParams[4] = _s;
- System.out.print("KFTrackParameters: setting fit pars = ");
- this.printFitParams();
+ _l = pathLength;
+ _alpha = _l*_omega;
+
+ //System.out.print("KFTrackParameters: setting fit pars = ");
+ //this.printFitParams();
makePointFromFitRepresentation();
- makeHelixFromPointRepresentation();
+ makeHelixFromFitRepresentation();
+
+ //System.out.print("KFTrackParameters: control fit pars = ");
+ //this.printFitParams();
+
}
public void setTrackPointParameters(int ch, double x, double y, double z, double px, double py, double pz)
@@ -74,14 +68,13 @@
_charge = ch;
_x = x;
_y = y;
- _r = Math.sqrt(_x*_x+_y*_y);
- //_phi = 0.; if(_x!=0.&&_y!=0) { _phi = _r*Math.atan2(_y,_x); }
- _phi = 0.; if(_x!=0.&&_y!=0) { _phi = Math.atan2(_y,_x); }
_z = z;
_px = px;
_py = py;
_pz = pz;
+ _p = Math.sqrt(_px*_px + _py*_py + _pz*_pz);
+
_pointParams[0] = _x;
_pointParams[1] = _y;
_pointParams[2] = _z;
@@ -89,30 +82,45 @@
_pointParams[4] = _py;
_pointParams[5] = _pz;
- _pFitParams[0] = _phi;
- _pFitParams[1] = _z;
- _pFitParams[2] = _px;
- _pFitParams[3] = _py;
- _pFitParams[4] = _pz;
+ makeFitFromPointRepresentation();
+ makeHelixFromPointRepresentation();
+
+ //System.out.print("KFTrackParameters: control fit pars = ");
+ //this.printFitParams();
- _pt = Math.sqrt(_px*_px + _py*_py);
+ }
+
+ private void makePointFromFitRepresentation()
+ {
+ _pt = _cst.pointThreeTimesB()/Math.abs(_omega);
+ _px = _pt*(Math.cos(_alpha)*Math.cos(_phi0)-Math.sin(_alpha)*Math.sin(_phi0));
+ _py = _pt*(Math.sin(_alpha)*Math.cos(_phi0)+Math.cos(_alpha)*Math.sin(_phi0));
+ _pz = _pt*_s;
_p = Math.sqrt(_px*_px + _py*_py + _pz*_pz);
- _xc = _x+_charge*_py/_cst.pointThreeTimesB();
- _yc = _y-_charge*_px/_cst.pointThreeTimesB();
- _rc = Math.sqrt((_x-_xc)*(_x-_xc)+(_y-_yc)*(_y-_yc));
- _pt = Math.sqrt(_px*_px+_py*_py); // = _rc/_cst.pointThreeTimesB();
- _s = _pz/_pt;
+ _x = _d0*_charge*Math.sin(_phi0)+(_charge/Math.abs(_omega))*
+ (Math.sin(_phi0)-Math.cos(_alpha)*Math.sin(_phi0)-Math.sin(_alpha)*Math.cos(_phi0));
+ _y = -_d0*_charge*Math.cos(_phi0)-(_charge/Math.abs(_omega))*
+ (Math.cos(_phi0)+Math.sin(_alpha)*Math.sin(_phi0)-Math.cos(_alpha)*Math.cos(_phi0));
+ _z = _z0+_alpha*_s/Math.abs(_omega);
- double r0 = Math.sqrt(_xc*_xc+_yc*_yc);
- double alpha = Math.acos((r0*r0-_x*_xc-_y*_yc)/(r0*_rc));
- double z0Plus = z+alpha*_pz/_cst.pointThreeTimesB();
- double z0Minus = z-alpha*_pz/_cst.pointThreeTimesB();
- if((_z-z0Plus)*_pz>=0.) {
- _z0 = z0Plus;
- } else {
- _z0 = z0Minus;
- }
+ _pointParams[0] = _x;
+ _pointParams[1] = _y;
+ _pointParams[2] = _z;
+ _pointParams[3] = _px;
+ _pointParams[4] = _py;
+ _pointParams[5] = _pz;
+ }
+
+ private void makeHelixFromFitRepresentation()
+ {
+ _xc = _charge*Math.sin(_phi0)*(_d0+1./Math.abs(_omega));
+ _yc = -_charge*Math.cos(_phi0)*(_d0+1./Math.abs(_omega));
+ _r0 = Math.sqrt(_xc*_xc+_yc*_yc);
+ _rc = 1./Math.abs(_omega);
+ _pt = _cst.pointThreeTimesB()/Math.abs(_omega);
+ //_s = _s;
+ //_z0 = _z0;
_helixParams[0] = _xc;
_helixParams[1] = _yc;
@@ -120,37 +128,23 @@
_helixParams[3] = _z0;
_helixParams[4] = _s;
- makeFitFromHelixRepresentation();
- makePointFromFitRepresentation();
- }
-
- private void makePointFromFitRepresentation()
- {
- _pt = _cst.pointThreeTimesB()/Math.abs(_omega);
- _px0 = _pt*Math.cos(_phi0);
- _py0 = _pt*Math.sin(_phi0);
- _pz0 = _pt*_s;
- _p = Math.sqrt(_px0*_px0 + _py0*_py0 + _pz0*_pz0);
-
- _x0 = _d0*_charge*Math.sin(_phi0);
- _y0 = -_d0*_charge*Math.cos(_phi0);
-
- _pointOriginParams[0] = _x0;
- _pointOriginParams[1] = _y0;
- _pointOriginParams[2] = _z0;
- _pointOriginParams[3] = _px0;
- _pointOriginParams[4] = _py0;
- _pointOriginParams[5] = _pz0;
+ //System.out.print("KFTrackParameters: helix par = ");
+ //System.out.print(_xc); System.out.print(" ");
+ //System.out.print(_yc); System.out.print(" ");
+ //System.out.print(_rc); System.out.print(" ");
+ //System.out.print(_z0); System.out.print(" ");
+ //System.out.println(_s);
}
private void makeHelixFromPointRepresentation()
{
- _xc = _x0+_charge*_py0/_cst.pointThreeTimesB();
- _yc = _y0-_charge*_px0/_cst.pointThreeTimesB();
- _rc = Math.sqrt((_x0-_xc)*(_x0-_xc)+(_y0-_yc)*(_y0-_yc));
- _pt = Math.sqrt(_px0*_px0+_py0*_py0); // = _rc/_cst.pointThreeTimesB();
- _s = _pz0/_pt;
- //_z0 = _z0; // valid only if point is POCA => fix me!
+ _xc = _x+_charge*_py/_cst.pointThreeTimesB();
+ _yc = _y-_charge*_px/_cst.pointThreeTimesB();
+ _r0 = Math.sqrt(_xc*_xc+_yc*_yc);
+ _rc = Math.sqrt((_x-_xc)*(_x-_xc)+(_y-_yc)*(_y-_yc));
+ _pt = Math.sqrt(_px*_px+_py*_py); // = _rc/_cst.pointThreeTimesB();
+ _s = _pz/_pt;
+ _z0 = _z+Math.atan2(_y*_xc-_x*_yc,_r0*_r0-_x*_xc-_y*_yc)*_pz/_cst.pointThreeTimesB();
_helixParams[0] = _xc;
_helixParams[1] = _yc;
@@ -158,12 +152,38 @@
_helixParams[3] = _z0;
_helixParams[4] = _s;
+ //System.out.print("KFTrackParameters: helix par = ");
+ //System.out.print(_xc); System.out.print(" ");
+ //System.out.print(_yc); System.out.print(" ");
+ //System.out.print(_rc); System.out.print(" ");
+ //System.out.print(_z0); System.out.print(" ");
+ //System.out.println(_s);
+ }
+
+ private void makeFitFromPointRepresentation()
+ {
+ double tmpXc = _x+_charge*_py/_cst.pointThreeTimesB();
+ double tmpYc = _y-_charge*_px/_cst.pointThreeTimesB();
+ double tmpPt = Math.sqrt(_px*_px+_py*_py);
+ _d0 = Math.sqrt(tmpXc*tmpXc+tmpYc*tmpYc)-tmpPt/_cst.pointThreeTimesB();
+ _phi0 = Math.atan2(_charge*tmpXc,-_charge*tmpYc);
+ _omega = -_charge*_cst.pointThreeTimesB()/tmpPt;
+ _alpha = -Math.atan2(_y*_py+_x*_px,_charge*tmpPt*tmpPt/_cst.pointThreeTimesB()+_x*_py-_y*_px);
+ _z0 = _z-_alpha*_pz/_cst.pointThreeTimesB();
+ _s = _pz/tmpPt;
+ _l = _alpha/_omega;
+
+ _fitParams[0] = _d0;
+ _fitParams[1] = _phi0;
+ _fitParams[2] = _omega;
+ _fitParams[3] = _z0;
+ _fitParams[4] = _s;
}
private void makeFitFromHelixRepresentation()
{
_d0 = Math.sqrt(_xc*_xc+_yc*_yc)-_rc;
- _phi0 = Math.atan2(_xc,-_yc);
+ _phi0 = Math.atan2(_charge*_xc,-_charge*_yc);
_omega = -_charge/_rc;
//_z0 = _z0;
//_s = _s;
@@ -177,233 +197,73 @@
public KFTrackParameters makePredictionAtRadius(double newRadius) {
- System.out.print("KFTrackParameters: pre-propagate "); this.printFitParams();
+ // in fit representation, only the path length changes when propagating
+ // to a specific radius
- KFTrackParameters newTrack = new KFTrackParameters();
- newTrack = propagate(newRadius);
-
- System.out.print("KFTrackParameters: propagated "); newTrack.printFitParams();
-
- // propagate error matrix
- Matrix jacobian = getJacobian(newRadius);
- Matrix transposedJacobian = jacobian.transpose();
- this._ErrorMatrixAtRadiusRprime = jacobian.times(_ErrorMatrixAtRadiusR).times(transposedJacobian);
- newTrack.setErrorMatrixAtRadiusR(_ErrorMatrixAtRadiusRprime);
-
- System.out.println("KFTrackParameters: Error matrix at radius R:");
- _ErrorMatrixAtRadiusR.print(15,12);
-
- System.out.println("KFTrackParameters: Error matrix at radius Rprime:");
- _ErrorMatrixAtRadiusRprime.print(15,12);
+ double newPathLength = propagate(newRadius);
- System.out.println("KFTrackParameters: Jacobian matrix at radius R:");
- jacobian.print(15,12);
+ KFTrackParameters newTrack = new KFTrackParameters();
+ newTrack.setTrackFitParameters(_d0,_phi0,_omega,_z0,_s,newPathLength);
+ newTrack.setErrorMatrix(_FitErrorMatrix);
- // test linearity of jacobian
- Matrix stateVec = new Matrix(_pFitParams,5);
- Matrix pred = jacobian.times(stateVec);
- pred.print(10,5);
+ //System.out.print("KFTrackParameters: propagated "); newTrack.printFitParams();
return newTrack;
}
- public KFTrackParameters propagate(double r)
+ public double propagate(double r)
{
- // propagate trajectory to radius r
- double r0Sqr = _xc*_xc+_yc*_yc;
- double deltaSqr = 0.5*(r*r+r0Sqr-_rc*_rc);
- //double Delta = deltaSqr*deltaSqr*(_xc*_xc-r0Sqr)+r0Sqr*r*r*_yc*_yc;
- double Delta = (r0Sqr*r*r-deltaSqr*deltaSqr)*_yc*_yc;
- double sqrtDelta = 0;
- if(Delta>=0) {
- sqrtDelta = Math.sqrt(Delta);
- } else {
- System.out.println("Negative sqrt");
- }
- double xPlus = (deltaSqr*_xc + sqrtDelta)/r0Sqr;
- double xMinus = (deltaSqr*_xc - sqrtDelta)/r0Sqr;
-
- double yPlus = (deltaSqr-xPlus*_xc)/_yc;
- double yMinus = (deltaSqr-xMinus*_xc)/_yc;
-
- double rCurrent = Math.sqrt(_x*_x+_y*_y);
- double testPlus = _px*(xPlus-_x)+_py*(yPlus-_y);
- double testMinus = _px*(xMinus-_x)+_py*(yMinus-_y);
-
- double xNew = _x;
- double yNew = _y;
- if(testPlus>testMinus) {
- xNew = xPlus;
- yNew = yPlus;
- } else {
- xNew = xMinus;
- yNew = yMinus;
- }
-
- //System.out.print("Old radius = "); System.out.println(_rc);
- //double newRC = Math.sqrt((xNew-_xc)*(xNew-_xc)+(yNew-_yc)*(yNew-_yc));
- //System.out.print("New radius = "); System.out.println(newRC);
-
- double cosAlpha = ((_x-_xc)*(xNew-_xc)+(_y-_yc)*(yNew-_yc))/(_rc*_rc);
- double sinAlpha = ((_x-_xc)*(yNew-_yc)-(_y-_yc)*(xNew-_xc))/(_rc*_rc);
- double alpha = Math.atan2(sinAlpha,cosAlpha);
-
- //System.out.print("xc,yc = "); System.out.print(_xc); System.out.print(" "); System.out.println(_yc);
- //System.out.print("x,y = "); System.out.print(_x); System.out.print(" "); System.out.println(_y);
- //System.out.print("xNew,yNew = "); System.out.print(xNew); System.out.print(" "); System.out.println(yNew);
- //System.out.print("cos,sin = "); System.out.print(cosAlpha); System.out.print(" "); System.out.println(sinAlpha);
- //System.out.print("alpha = "); System.out.println(alpha);
-
- double zNew = _z + alpha*_pz/_cst.pointThreeTimesB();
-
- double pxNew = _px*Math.cos(alpha) - _py*Math.sin(alpha);
- double pyNew = _px*Math.sin(alpha) + _py*Math.cos(alpha);
- double pzNew = _pz;
-
- KFTrackParameters newTrackPoint = new KFTrackParameters();
- newTrackPoint.setTrackPointParameters(_charge, xNew, yNew, zNew, pxNew, pyNew, pzNew);
-
- return newTrackPoint;
-
- }
-
- private Matrix getJacobian(double r)
- {
- // calculate the jacobian (dphi/dphi, dphi/dz, ...)
- double[][] Jacobian = new double[5][5];
-
- //double[] variations = new double[] {0.00001,0.001,0.0000001,0.001,0.001};
- double[] variations = new double[] {0.001,0.001,0.0001,0.0001,0.0001};
-
- double[] fitPars = new double[5];
- for(int k=0; k<5; k++) { fitPars[k] = _pFitParams[k]; }
-
- for(int j=0; j<5; j++) {
-
- double p1p = fitPars[0]; double p1m = fitPars[0];
- if(j==0) { p1p += 0.5*variations[0]; p1m -= 0.5*variations[0]; }
- double p2p = fitPars[1]; double p2m = fitPars[1];
- if(j==1) { p2p += 0.5*variations[1]; p2m -= 0.5*variations[1]; }
- double p3p = fitPars[2]; double p3m = fitPars[2];
- if(j==2) { p3p += 0.5*variations[2]; p3m -= 0.5*variations[2]; }
- double p4p = fitPars[3]; double p4m = fitPars[3];
- if(j==3) { p4p += 0.5*variations[3]; p4m -= 0.5*variations[3]; }
- double p5p = fitPars[4]; double p5m = fitPars[4];
- if(j==4) { p5p += 0.5*variations[4]; p5m -= 0.5*variations[4]; }
-
- double xp = _r*Math.cos(p1p); double yp = _r*Math.sin(p1p);
- double xm = _r*Math.cos(p1m); double ym = _r*Math.sin(p1m);
-
- System.out.print("KFTrackParameters d0/phi0/omega...Plus= ");
- System.out.print(p1p); System.out.print(" (");
- System.out.print(xp); System.out.print(" , ");
- System.out.print(yp); System.out.print(") ");
- System.out.print(p2p); System.out.print(" ");
- System.out.print(p3p); System.out.print(" ");
- System.out.print(p4p); System.out.print(" ");
- System.out.println(p5p);
-
- System.out.print("KFTrackParameters d0/phi0/omega...Minus= ");
- System.out.print(p1m); System.out.print(" (");
- System.out.print(xm); System.out.print(" , ");
- System.out.print(ym); System.out.print(") ");
- System.out.print(p2m); System.out.print(" ");
- System.out.print(p3m); System.out.print(" ");
- System.out.print(p4m); System.out.print(" ");
- System.out.println(p5m);
-
- KFTrackParameters trkVarPlus = new KFTrackParameters();
- trkVarPlus.setTrackPointParameters(_charge, xp, yp, p2p, p3p, p4p, p5p);
- KFTrackParameters trkVarPlusPropagated = trkVarPlus.propagate(r);
- double[] trkVarPlusParam = trkVarPlusPropagated.getPFitParameters();
- KFTrackParameters trkVarMinus = new KFTrackParameters();
- trkVarMinus.setTrackPointParameters(_charge, xm, ym, p2m, p3m, p4m, p5m);
- KFTrackParameters trkVarMinusPropagated = trkVarMinus.propagate(r);
- double[] trkVarMinusParam = trkVarMinusPropagated.getPFitParameters();
-
- for(int i=0; i<5; i++) {
- //Jacobian[j][i] = (trkVarPlusParam[i]-trkVarMinusParam[i])/variations[j];
- Jacobian[i][j] = (trkVarPlusParam[i]-trkVarMinusParam[i])/variations[j];
- }
- }
-
-
- Matrix jacobianMatrix = new Matrix(Jacobian);
-
- return jacobianMatrix;
-
+ double par = 1.-((r*r-_d0*_d0)*_omega*_omega)/(2.*(1.+_d0*Math.abs(_omega)));
+ double newl = Math.abs(Math.acos(par)/_omega);
+ return newl;
}
public void setPointToRadius(double r) {
- this.setPointToPOCA();
- KFTrackParameters trk = propagate(r);
- double[] nP = trk.getPointParameters();
- for(int i=0; i<6; i++) { _pointParams[i] = nP[i]; }
- _x = nP[0]; _y = nP[1]; _z = nP[2]; _px = nP[3]; _py = nP[4]; _pz = nP[5];
- _r = Math.sqrt(_x*_x+_y*_y);
- //_phi = 0.; if(_x!=0.&&_y!=0) { _phi = _r*Math.atan2(_y,_x); }
- _phi = 0.; if(_x!=0.&&_y!=0) { _phi = Math.atan2(_y,_x); }
- _pFitParams[0] = _phi;
- for(int i=2; i<6; i++) { _pFitParams[i-1] = nP[i]; }
+ double newPathLength = propagate(r);
+ _l = newPathLength;
+ makePointFromFitRepresentation();
+ makeHelixFromFitRepresentation();
}
public void setPointToPOCA() {
- for(int i=0; i<6; i++) { _pointParams[i] = _pointOriginParams[i]; }
- _x = _x0; _y = _y0; _z = _z0; _px = _px0; _py = _py0; _pz = _pz0;
- _r = Math.sqrt(_x*_x+_y*_y);
- //_phi = 0.; if(_x!=0.&&_y!=0) { _phi = _r*Math.atan2(_y,_x); }
- _phi = 0.; if(_x!=0.&&_y!=0) { _phi = Math.atan2(_y,_x); }
- _pFitParams[0] = _phi;
- for(int i=2; i<6; i++) { _pFitParams[i-1] = _pointOriginParams[i]; }
+ _l = 0.;
+ makePointFromFitRepresentation();
+ makeHelixFromFitRepresentation();
}
public double[] getFitParameters() {return _fitParams;}
- public double[] getPointOriginParameters() {return _pointOriginParams;}
public double[] getPointParameters() {return _pointParams;}
- public double[] getPFitParameters() {return _pFitParams;}
public double[] getHelixParameters() {return _helixParams;}
public int getCharge() {return _charge;}
+ public double getPathLength() {return _l;}
public double getPt() {return _pt;}
+ public double getPtot() {return _p;}
+ public double getPhi() {return Math.atan2(_y,_x);}
+ public double getZ() {return _z;}
- public void setErrorMatrixAtRadiusR(Matrix cov) { _ErrorMatrixAtRadiusR = cov; }
+ public void setErrorMatrix(Matrix cov) { _FitErrorMatrix = cov; }
- public Matrix getFitCovMatrix() {return _FitErrorMatrix;}
- public Matrix getCovMatrixAtRadiusR() {return _ErrorMatrixAtRadiusR;}
- public Matrix getCovMatrixAtRadiusRprime() {return _ErrorMatrixAtRadiusRprime;}
+ public Matrix getErrorMatrix() {return _FitErrorMatrix;}
private Matrix _FitErrorMatrix = new Matrix(5,5);
- private Matrix _ErrorMatrixAtRadiusR = new Matrix(5,5);
- private Matrix _ErrorMatrixAtRadiusRprime = new Matrix(5,5);
+
private KFConstants _cst = new KFConstants();
public void resetTrackParameters() {
int i;
_charge = 0;
_d0 = 0.; _phi0 = 0.; _omega = 0.; _z0 = 0.; _s = 0.;
+ _l = 0; _alpha = 0.;
for(i=0;i<5;i++){_fitParams[i]=0.;};
- _p = 0.; _pt = 0.; _px0 = 0.; _py0 = 0.; _pz0 = 0.;
- _x0 = 0.; _y0 = 0.;
- for(i=0;i<6;i++){_pointOriginParams[i]=0.;};
-
- _xc = 0.; _yc = 0.; _rc = 0.;
- for(i=0;i<5;i++){_helixParams[i]=0.;};
-
- _px = 0.; _py = 0.; _pz = 0.;
- _x = 0.; _y = 0.; _z = 0.; _r = 0.; _phi = 0.;
+ _p = 0.; _pt = 0.; _px = 0.; _py = 0.; _pz = 0.;
+ _x = 0.; _y = 0.; _z = 0.;
for(i=0;i<6;i++){_pointParams[i]=0.;};
- for(i=0;i<5;i++){_pFitParams[i]=0.;};
- }
+ _xc = 0.; _yc = 0.; _rc = 0.; _r0 = 0.;
+ for(i=0;i<5;i++){_helixParams[i]=0.;};
- public void printPFitParams() {
- System.out.print("PFit params = ");
- System.out.print(_pFitParams[0]); System.out.print(" ");
- System.out.print(_pFitParams[1]); System.out.print(" ");
- System.out.print(_pFitParams[2]); System.out.print(" ");
- System.out.print(_pFitParams[3]); System.out.print(" ");
- System.out.println(_pFitParams[4]);
}
public void printFitParams() {
@@ -412,6 +272,20 @@
System.out.print(_fitParams[1]); System.out.print(" ");
System.out.print(_fitParams[2]); System.out.print(" ");
System.out.print(_fitParams[3]); System.out.print(" ");
- System.out.println(_fitParams[4]);
+ System.out.print(_fitParams[4]); System.out.print(" ");
+ //System.out.print(_alpha); System.out.print(" ");
+ System.out.println(_l);
+ }
+
+ public void printPointParams() {
+ System.out.print("Point params = ");
+ System.out.print(_pointParams[0]); System.out.print(" ");
+ System.out.print(_pointParams[1]); System.out.print(" ");
+ System.out.print(_pointParams[2]); System.out.print(" ");
+ System.out.print(_pointParams[3]); System.out.print(" ");
+ System.out.print(_pointParams[4]); System.out.print(" ");
+ System.out.print(_pointParams[5]); System.out.print(" ");
+ //System.out.print(_alpha); System.out.print(" ");
+ System.out.println(getPhi());
}
}
CVSspam 0.2.8