hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.1 -r1.2
--- GBLFileIO.java 14 Jul 2013 16:58:13 -0000 1.1
+++ GBLFileIO.java 13 Aug 2013 22:24:16 -0000 1.2
@@ -74,6 +74,14 @@
addLine(this.getPerTrackParamStr(perPar));
}
+ String getPerTrackParamTruthStr(PerigeeParams perPar) {
+ return String.format("Truth perPar (kappa theta phi d0 z0) %f %f %f %f %f",perPar.getKappa(),perPar.getTheta(),perPar.getPhi(),perPar.getD0(),perPar.getZ0());
+ }
+
+ void printPerTrackParamTruth(PerigeeParams perPar) {
+ addLine(this.getPerTrackParamTruthStr(perPar));
+ }
+
String getClTrackParamStr(ClParams perPar) {
return String.format("Track clPar (q/p lambda phi xT yT) %f %f %f %f %f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
}
@@ -131,11 +139,12 @@
addLine(str);
}
- void printPerTrackCov(SymmetricMatrix cov) {
- String str = "Track perCov ";
+ void printPerTrackCov(HelicalTrackFit htf) {
+ String str = "Track perCov (idx: dca,phi0,curv,z0,slope) ";
+ SymmetricMatrix cov = htf.covariance();
for(int irow=0;irow<cov.getNRows();++irow) {
for(int icol=0;icol<cov.getNColumns();++icol) {
- str += String.format("%f ", cov.e(irow, icol));
+ str += String.format("%e ", cov.e(irow, icol));
}
}
addLine(str);
@@ -173,6 +182,9 @@
addLine(String.format("Strip ures %f %f", ures, uresErr));
}
+ void printStripMeasTruth(double ures, double uresErr) {
+ addLine(String.format("Strip truth ures %f %f", ures, uresErr));
+ }
void printStripScat(double scatAngle) {
addLine(String.format("Strip scatangle %f",scatAngle));
}
@@ -181,6 +193,15 @@
addLine(String.format("Strip meas dir %f %f %f", u.x(), u.y(), u.z()));
}
+ void printNonMeasDir(Hep3Vector u) {
+ addLine(String.format("Strip non-meas dir %f %f %f", u.x(), u.y(), u.z()));
+ }
+
+ void printNormalDir(Hep3Vector u) {
+ addLine(String.format("Strip normal dir %f %f %f", u.x(), u.y(), u.z()));
+ }
+
+
void printMomentum(double p, double p_truth) {
addLine(String.format("Track mom %f %f", p, p_truth));
}
@@ -195,5 +216,10 @@
addLine(str);
}
+ void printChi2(double[] chisq, int[] ndf) {
+ addLine(String.format("Track chi2/ndf (circle,zfit) %f %d %f %d",chisq[0],ndf[0],chisq[1],ndf[1]));
+ }
+
+
}
hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.1 -r1.2
--- GBLOutput.java 14 Jul 2013 16:58:14 -0000 1.1
+++ GBLOutput.java 13 Aug 2013 22:24:16 -0000 1.2
@@ -104,6 +104,9 @@
//BasicMatrix clPar = (BasicMatrix) MatrixOp.mult(jacPerToCl, MatrixOp.transposed(perPar));
//clPar = (BasicMatrix) MatrixOp.transposed(clPar);
+ // print chi2 of fit
+ file.printChi2(htf.chisq(),htf.ndf());
+
// print track parameters for easier debugging later
file.printOldPerTrackParam(htf);
if(_debug>0) {
@@ -114,25 +117,43 @@
file.printClTrackParam(clPar);
List<MCParticle> truthParticles = getTruthParticle(trk.getCharge()>0?11:-11,mcParticles);
- double p_truth = 0.;
double dp_min = 99999999.9;
+ MCParticle mcp = null;
for(MCParticle p : truthParticles) {
double dp_loop = Math.abs(p.getMomentum().magnitude() - getMomentum(trk).magnitude());
if(dp_loop < dp_min) {
dp_min = dp_loop;
- p_truth = p.getMomentum().magnitude();
+ mcp = p;
}
}
- file.printMomentum(htf.p(this._B.magnitude()),p_truth);
+ file.printMomentum(htf.p(this._B.magnitude()),mcp.getMomentum().magnitude());
+
+ // Get track paramteres from MC particle
+ HelicalTrackFit htfTruth = getHTF(mcp);
+ PerigeeParams perParTruth = new PerigeeParams(htfTruth);
+ file.printPerTrackParamTruth(perParTruth);
+
+ if(_debug>0) {
+ System.out.printf("%s\n",file.getPerTrackParamTruthStr(perParTruth));
+ }
+
// find the projection from the I,J,K to U,V,T coordinates
Hep3Matrix perToClPrj = this.getPerToClPrj(htf);
file.printPerToClPrj(perToClPrj);
- // covariance matrix
- SymmetricMatrix perCov = htf.covariance();
- file.printPerTrackCov(perCov);
+ // covariance matrix from the helix fit
+ file.printPerTrackCov(htf);
+
+ // calculate truth chi2
+ double chi2 = truthTrackFitChi2(perPar,perParTruth,htf.covariance());
+
+ //if(_debug>0)
+ System.out.printf("%s: perPar covariance matrix\n%s\n",this.getClass().getSimpleName(),htf.covariance().toString());
+ System.out.printf("%s: truth perPar chi2 %f\n",this.getClass().getSimpleName(),chi2);
+
+
BasicMatrix clCov = new BasicMatrix(5,5);
this.initUnit(clCov);
clCov = (BasicMatrix) MatrixOp.mult(0.1*0.1,clCov);
@@ -160,6 +181,16 @@
//Find intercept point with sensor in tracking frame
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z()));
+ if(_debug>0) {
+ System.out.printf("trkpos at intercept %s\n",trkpos.toString());
+ }
+
+ Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
+
+ if(_debug>0) {
+ System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
+ }
+
//path length to intercept
double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
double s_hit = HelixUtils.PathLength(htf, hit); //cross-check
@@ -168,9 +199,10 @@
file.printStripPathLen(s);
//print stereo angle in YZ plane
- if(_debug>0) System.out.printf("u = %s v = %s \n",strip.u().toString(),strip.v().toString());
- double stereoAngle = Math.atan(strip.v().z());
+ if(_debug>0) System.out.printf("u = %s v = %s w = %s \n",strip.u().toString(),strip.v().toString(), strip.w().toString());
file.printMeasDir(strip.u());
+ file.printNonMeasDir(strip.v());
+ file.printNormalDir(strip.w());
//Print track direction at intercept
Hep3Vector tDir = HelixUtils.Direction(htf, s);
@@ -182,10 +214,12 @@
// start by find the distance vector between the center and the track position
Hep3Vector vdiffTrk = VecOp.sub(trkpos, origin);
+ Hep3Vector vdiffTrkTruth = VecOp.sub(trkposTruth, origin);
// then find the rotation from tracking to measurement frame
Hep3Matrix trkToStripRot = _trackerHitUtils.getTrackToStripRotation(strip);
// then rotate that vector into the measurement frame to get the predicted measurement position
Hep3Vector trkpos_meas = VecOp.mult(trkToStripRot, vdiffTrk);
+ Hep3Vector trkposTruth_meas = VecOp.mult(trkToStripRot, vdiffTrkTruth);
// hit uncertainty in measurement frame
double umeas = strip.umeas();
@@ -200,12 +234,15 @@
// residual in measurement frame
Hep3Vector res_meas = VecOp.sub(m_meas, trkpos_meas);
+ Hep3Vector resTruth_meas = VecOp.sub(m_meas, trkposTruth_meas);
// residual uncertainty in measurement frame
Hep3Vector res_err_meas = new BasicHep3Vector(uError,vError,wError);
+ Hep3Vector resTruth_err_meas = new BasicHep3Vector(uError,vError,wError);
// print measurement to file
file.printStripMeas(res_meas.x(),res_err_meas.x());
+ file.printStripMeasTruth(resTruth_meas.x(),resTruth_err_meas.x());
// print scattering angle
@@ -426,6 +463,42 @@
}
}
+ private HelicalTrackFit getHTF(MCParticle mcp) {
+ Hep3Vector org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
+ Hep3Vector p = this._hpstrans.transformVectorToTracking(mcp.getMomentum());
+ HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1*((int)mcp.getCharge()), -1.0*this._B.z());
+ double par[] = new double[5];
+ par[HelicalTrackFit.dcaIndex] = helixParamCalculator.getDCA();
+ par[HelicalTrackFit.slopeIndex] = helixParamCalculator.getSlopeSZPlane();
+ par[HelicalTrackFit.phi0Index] = helixParamCalculator.getPhi0();
+ par[HelicalTrackFit.curvatureIndex] = 1.0/helixParamCalculator.getRadius();
+ par[HelicalTrackFit.z0Index] = helixParamCalculator.getZ0();
+ HelicalTrackFit htf = new HelicalTrackFit(par, null, new double[2], new int[2], null, null);
+ return htf;
+ }
+
+ private double truthTrackFitChi2(PerigeeParams perPar, PerigeeParams perParTruth, SymmetricMatrix covariance) {
+ //re-shuffle the param vector to match the covariance order of parameters
+ BasicMatrix p = new BasicMatrix(1,5);
+ p.setElement(0, 0, perPar.getD0());
+ p.setElement(0, 1, perPar.getPhi());
+ p.setElement(0, 2, perPar.getKappa());
+ p.setElement(0, 0, perPar.getZ0());
+ p.setElement(0, 4, Math.tan(Math.PI/2.0-perPar.getTheta()));
+ BasicMatrix pt = new BasicMatrix(1,5);
+ pt.setElement(0, 0, perParTruth.getD0());
+ pt.setElement(0, 1, perParTruth.getPhi());
+ pt.setElement(0, 2, perParTruth.getKappa());
+ pt.setElement(0, 0, perParTruth.getZ0());
+ pt.setElement(0, 4, Math.tan(Math.PI/2.0-perParTruth.getTheta()));
+ BasicMatrix res = (BasicMatrix) MatrixOp.sub(p, pt);
+ BasicMatrix chi2 = (BasicMatrix) MatrixOp.mult(res,MatrixOp.mult(covariance, MatrixOp.transposed(res)));
+ if(chi2.getNColumns()!=1 ||chi2.getNRows()!=1) {
+ throw new RuntimeException("matrix dim is screwed up!");
+ }
+ return chi2.e(0, 0);
+ }
+