hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.8 -r1.9
--- GBLOutput.java 26 Aug 2013 21:40:03 -0000 1.8
+++ GBLOutput.java 5 Sep 2013 00:05:14 -0000 1.9
@@ -88,8 +88,24 @@
SeedTrack st = (SeedTrack)trk;
SeedCandidate seed = st.getSeedCandidate();
- HelicalTrackFit htf = seed.getHelix();
+ HelicalTrackFit htf = seed.getHelix();
+ double[] htf_params = new double[5];
+ System.arraycopy(htf.parameters(), 0, htf_params, 0, 5);
+ htf_params[0] *= -1;
+ HelicalTrackFit htf_dca_flip = new HelicalTrackFit(htf_params, htf.covariance(), htf.chisq(), htf.ndf(), htf.PathMap(), htf.ScatterMap());
+ //htf = htf_dca_flip;
+
List<HelicalTrackHit> hits = seed.getHits();
+
+ // Find the truth particle of the track
+ MCParticle mcp = getMatchedTruthParticle(trk);
+
+ // Get track paramteres from MC particle
+ HelicalTrackFit htfTruth = getHTF(mcp);
+
+ // Use the truth helix as the initial track for GBL?
+ //htf = htfTruth;
+
if (this._debug>1) {
System.out.printf("%s: find scatters\n",this.getClass().getSimpleName());
@@ -132,7 +148,7 @@
System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
}
- MCParticle mcp = getMatchedTruthParticle(trk);
+
// cross-check
if(!mcp_pair.contains(mcp)) {
@@ -177,8 +193,7 @@
}
- // Get track paramteres from MC particle
- HelicalTrackFit htfTruth = getHTF(mcp);
+
PerigeeParams perParTruth = new PerigeeParams(htfTruth);
file.printPerTrackParamTruth(perParTruth);
ClParams clParTruth = new ClParams(htfTruth);
@@ -234,11 +249,24 @@
//Find intercept point with sensor in tracking frame
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z()));
Hep3Vector trkposXPlane = TrackUtils.getHelixXPlaneIntercept(htf, strip.w(), origin);
+ Hep3Vector trkposXPlaneDCAFlip = TrackUtils.getHelixXPlaneIntercept(htf_dca_flip, strip.w(), origin);
+
if(_debug>0) {
System.out.printf("trkpos at intercept %s\n",trkpos.toString());
System.out.printf("trkposXPlane at intercept %s\n",trkposXPlane.toString());
+ System.out.printf("trkposXPlaneDCAFlip at intercept %s\n",trkposXPlaneDCAFlip.toString());
+ // is it on the plane?
+ double trkpos_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkpos, strip.origin())));
+ double trkposXPlane_d = VecOp.dot(strip.w(),VecOp.unit(VecOp.sub(trkposXPlane, strip.origin())));
+ System.out.printf("trkpos_d %f with diff %s trkpos %s org %s \n",trkpos_d,VecOp.unit(VecOp.sub(trkpos, strip.origin())), trkpos.toString(),strip.origin().toString());
+ System.out.printf("trkposXPlane_d %f with diff %s trkpos %s org %s \n",trkposXPlane_d,VecOp.unit(VecOp.sub(trkposXPlane, strip.origin())), trkposXPlane.toString(),strip.origin().toString());
}
+ //trkpos = trkposXPlane;
+ //trkpos = trkposXPlaneDCAFlip;
+
+
+
Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
if(_debug>0) {
@@ -282,10 +310,26 @@
}
//path length to intercept
- double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
+ double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
double s3D = s / Math.cos(Math.atan(htf.slope()));
if(_debug>0) {
+
+ //Calculate path length by hand for x plane
+ {
+ double phi_at_x = Math.asin( 1/htf.R()* ( -1*htf.dca()*Math.sin(htf.phi0()) + htf.R()*Math.sin(htf.phi0()) - trkpos.x() ) );
+ double dphi_at_x = phi_at_x - htf.phi0();
+ //if(dphi_at_x > Math.PI) dphi_at_x -= 2.0*Math.PI;
+ double s_at_x = -1.0*dphi_at_x*htf.R();
+ double y_at_x = htf.dca()*Math.cos(htf.phi0()) - htf.R()*Math.cos(htf.phi0()) + htf.R()*Math.cos(phi_at_x);
+ double z_at_x = htf.z0() + s_at_x*htf.slope();
+ Hep3Vector trkposmanual_at_x = new BasicHep3Vector(trkpos.x(),y_at_x,z_at_x);
+ System.out.printf("s = %f trkpos %s (initial manual xplane)\n", s_at_x , trkposmanual_at_x.toString());
+
+ }
+
+
+
double s_XPlane = HelixUtils.PathToXPlane(htf,trkposXPlane.x(),0,0).get(0);
double s_truth = HelixUtils.PathToXPlane(htfTruth, trkposTruth.x(), 0, 0).get(0);
double s3D_truth = s_truth / Math.cos(Math.atan(htfTruth.slope()));
@@ -648,7 +692,9 @@
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);
+ SymmetricMatrix cov = new SymmetricMatrix(5);
+ for(int i=0;i<cov.getNRows();++i) cov.setElement(i, i, 1.);
+ HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
return htf;
}
@@ -811,8 +857,8 @@
Hep3Matrix perToClPrj = GBLOutput.this.getPerToClPrj(htf);
- double d0 = htf.dca();
- double z0 = htf.dca();
+ double d0 = -1 * htf.dca(); //sign convention for curvilinear frame
+ double z0 = htf.z0();
Hep3Vector vecPer = new BasicHep3Vector(0.,d0,z0);
//System.out.printf("%s: vecPer=%s\n",this.getClass().getSimpleName(),vecPer.toString());