hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.2 -r1.3
--- GBLFileIO.java 13 Aug 2013 22:24:16 -0000 1.2
+++ GBLFileIO.java 21 Aug 2013 02:44:57 -0000 1.3
@@ -81,7 +81,15 @@
void printPerTrackParamTruth(PerigeeParams perPar) {
addLine(this.getPerTrackParamTruthStr(perPar));
}
-
+
+ String getClTrackParamTruthStr(ClParams perPar) {
+ return String.format("Truth clPar (q/p lambda phi xT yT) %f %f %f %f %f",perPar.getQoverP(),perPar.getLambda(),perPar.getPhi(),perPar.getXt(),perPar.getYt());
+ }
+
+ void printClTrackParamTruth(ClParams perPar) {
+ addLine(this.getClTrackParamTruthStr(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());
}
@@ -166,6 +174,14 @@
addLine(str);
}
+ void printStripTrackDirFull(Hep3Vector dir) {
+ addLine(String.format("Strip track dir %f %f %f",dir.x(),dir.y(),dir.z()));
+ }
+
+ void printStripTrackPos(Hep3Vector pos) {
+ addLine(String.format("Strip track pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ }
+
void printStrip(int id, int layer) {
addLine(String.format("New Strip id layer %d %d", id,layer));
}
@@ -185,6 +201,11 @@
void printStripMeasTruth(double ures, double uresErr) {
addLine(String.format("Strip truth ures %f %f", ures, uresErr));
}
+
+ void printStripMeasSimHit(double ures, double uresErr) {
+ addLine(String.format("Strip simhit ures %f %f", ures, uresErr));
+ }
+
void printStripScat(double scatAngle) {
addLine(String.format("Strip scatangle %f",scatAngle));
}
@@ -220,6 +241,14 @@
addLine(String.format("Track chi2/ndf (circle,zfit) %f %d %f %d",chisq[0],ndf[0],chisq[1],ndf[1]));
}
+ void printOrigin(Hep3Vector pos) {
+ addLine(String.format("Strip origin pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ }
+
+ void printHitPos3D(Hep3Vector pos) {
+ addLine(String.format("Strip 3D hit pos %f %f %f",pos.x(),pos.y(),pos.z()));
+ }
+
}
hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.3 -r1.4
--- GBLOutput.java 14 Aug 2013 00:47:33 -0000 1.3
+++ GBLOutput.java 21 Aug 2013 02:44:57 -0000 1.4
@@ -18,6 +18,7 @@
import org.lcsim.constants.Constants;
import org.lcsim.event.MCParticle;
import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track;
import org.lcsim.fit.helicaltrack.*;
import org.lcsim.geometry.Detector;
@@ -62,7 +63,7 @@
_scattering = new MultipleScattering(_materialmanager);
_B = _hpstrans.transformVectorToTracking(bfield);
System.out.printf("%s: B field %s\n",this.getClass().getSimpleName(),_B.toString());
- _scattering.setBField(_B.z());
+ _scattering.setBField(Math.abs(_B.z())); // only absolute of B is needed as it's used for momentum calculation only
}
public void setDebug(int debug) {
_debug = debug;
@@ -84,7 +85,7 @@
}
- void printGBL(Track trk, List<MCParticle> mcParticles) {
+ void printGBL(Track trk, List<MCParticle> mcParticles, List<SimTrackerHit> simTrackerHits) {
SeedTrack st = (SeedTrack)trk;
SeedCandidate seed = st.getSeedCandidate();
@@ -95,6 +96,7 @@
System.out.printf("%s: find scatters\n",this.getClass().getSimpleName());
}
//find the scattering angle from material manager
+ _scattering.setDebug(this._debug>0?true:false);
ScatterPoints scatters = _scattering.FindHPSScatterPoints(htf);
@@ -117,26 +119,92 @@
file.printPerTrackParam(perPar);
file.printClTrackParam(clPar);
- List<MCParticle> truthParticles = getTruthParticle(trk.getCharge()>0?11:-11,mcParticles);
- double dp_min = 99999999.9;
+ // Use the sim tracker hits to find the truth particle
+ // this would give me the particle kinematics responsible for the track
+ // but it might have scattered in the target
+
+ Map<MCParticle, List<SimTrackerHit>> MCParticleHitMap = new HashMap<MCParticle, List<SimTrackerHit>>();
+ for (SimTrackerHit hit : simTrackerHits) {
+ List hitList = MCParticleHitMap.get(hit.getMCParticle());
+ if (hitList == null) {
+ hitList = new ArrayList<SimTrackerHit>();
+ MCParticleHitMap.put(hit.getMCParticle(), hitList);
+ }
+ hitList.add(hit);
+ }
+
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;
- mcp = p;
+ for(Map.Entry<MCParticle, List<SimTrackerHit>> part : MCParticleHitMap.entrySet()) {
+
+ System.out.printf("%s: %d p=%s parents ",this.getClass().getSimpleName(),part.getKey().getPDGID(),part.getKey().getMomentum().toString());
+ for(MCParticle parent : part.getKey().getParents()) {
+ System.out.printf(" [%d p=%s]",parent.getPDGID(),parent.getMomentum().toString());
+ }
+ System.out.printf(" with simhits:\n");
+ for(SimTrackerHit simhit : part.getValue()) {
+ System.out.printf("%s: sim hit at %s\n",this.getClass().getSimpleName(),simhit.getPositionVec().toString());
+ }
+ // check charge and type
+ if(Math.abs(part.getKey().getPDGID()) == 11) {
+ if((trk.getCharge() * part.getKey().getPDGID()) > 0) {
+ if(mcp != null) {
+ //use the one with closest momentum ... -> FIX THIS?!
+ double p = htf.p(this._B.magnitude());
+ if( Math.abs(p - part.getKey().getMomentum().magnitude()) > Math.abs(p - mcp.getMomentum().magnitude())) {
+ continue;
+ }
+ }
+ mcp = part.getKey();
+ }
+ }
+ }
+
+ if(mcp==null) {
+ System.out.printf("%s: no truth particle found!\n",this.getClass().getSimpleName());
+ for(MCParticle loop_p : mcParticles) {
+ System.out.printf("%s: %d p=%s parents ",this.getClass().getSimpleName(),loop_p.getPDGID(),loop_p.getMomentum().toString());
+ for(MCParticle parent : loop_p.getParents()) {
+ System.out.printf(" [%d p=%s]",parent.getPDGID(),parent.getMomentum().toString());
+ }
+ System.out.printf("\n");
}
+ for (SimTrackerHit hit : simTrackerHits) {
+ System.out.printf("%s sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),hit.getPositionVec().toString(),hit.getMCParticle().getPDGID(),hit.getMCParticle().getMomentum().toString());
+ }
+ System.exit(1);
}
- 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);
+ ClParams clParTruth = new ClParams(htfTruth);
+ file.printClTrackParamTruth(clParTruth);
if(_debug>0) {
System.out.printf("%s\n",file.getPerTrackParamTruthStr(perParTruth));
}
+
+
+ // use only the truth particle that should've given the track
+ List<SimTrackerHit> simHits = MCParticleHitMap.get(mcp);
+
+ if(_debug>-1) {
+ System.out.printf("%s: %d sim tracker hits to MC particle with p = %f (trk p = %f q=%f):\n",this.getClass().getSimpleName(),simHits.size(),mcp.getMomentum().magnitude(),htf.p(this._B.magnitude()),Math.signum(htf.R()));
+ for(SimTrackerHit sh : simHits) System.out.printf("%s: sim hit at layer %d and pos %s \n",this.getClass().getSimpleName(),sh.getIdentifierFieldValue("layer"),sh.getPositionVec().toString());
+ }
+
+ // map layer and sim hits
+ Map<Integer, SimTrackerHit> simHitsLayerMap = new HashMap<Integer, SimTrackerHit >();
+ for(SimTrackerHit sh : simHits) {
+ int layer = sh.getIdentifierFieldValue("layer");
+ if(simHitsLayerMap.containsKey(layer)) {
+ if(simHitsLayerMap.get(layer).getPathLength() < sh.getPathLength()) {
+ continue;
+ }
+ }
+ simHitsLayerMap.put(layer, sh);
+ }
// find the projection from the I,J,K to U,V,T coordinates
@@ -180,6 +248,11 @@
//Center of the sensor
Hep3Vector origin = strip.origin();
+ file.printOrigin(origin);
+
+ // associated 3D position of the crossing of this and it's stereo partner sensor
+ file.printHitPos3D(hit.getCorrectedPosition());
+
//Find intercept point with sensor in tracking frame
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(htf, strip, Math.abs(_B.z()));
@@ -188,11 +261,26 @@
}
Hep3Vector trkposTruth = TrackUtils.getHelixPlaneIntercept(htfTruth, strip, Math.abs(_B.z()));
-
+
if(_debug>0) {
System.out.printf("trkposTruth at intercept %s\n",trkposTruth.toString());
}
+ // Find the sim tracker hit for this layer
+ SimTrackerHit simHit = simHitsLayerMap.get(strip.layer());
+
+ if(simHit==null) {
+ System.out.printf("%s: no sim hit for strip hit at layer %d\n",this.getClass().getSimpleName(),strip.layer());
+ for(Map.Entry<Integer,SimTrackerHit> entry : simHitsLayerMap.entrySet()) {
+ SimTrackerHit simhit = entry.getValue();
+ System.out.printf("%s layer %d sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),entry.getKey(),simhit.getPositionVec().toString(),simhit.getMCParticle().getPDGID(),simhit.getMCParticle().getMomentum().toString());
+ }
+ for (SimTrackerHit simhit : simTrackerHits) {
+ System.out.printf("%s sim hit at %s with MC particle pdgid %d with p %s \n",this.getClass().getSimpleName(),simhit.getPositionVec().toString(),simhit.getMCParticle().getPDGID(),simhit.getMCParticle().getMomentum().toString());
+ }
+ //System.exit(1);
+ }
+
//path length to intercept
double s = HelixUtils.PathToXPlane(htf,trkpos.x(),0,0).get(0);
double s_hit = HelixUtils.PathLength(htf, hit); //cross-check
@@ -211,18 +299,24 @@
double phi = htf.phi0() - s/htf.R();
double lambda = Math.atan(htf.slope());
file.printStripTrackDir(Math.sin(phi),Math.sin(lambda));
-
+ file.printStripTrackDirFull(tDir);
+ file.printStripTrackPos(trkpos);
+
//Print residual in measurment system
// 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();
double uError = strip.du();
@@ -238,13 +332,27 @@
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);
+ Hep3Vector resSimHit_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());
+ if(simHit!=null) {
+ Hep3Vector simHitPos = null;simHitPos = this._hpstrans.transformVectorToTracking(simHit.getPositionVec());
+ if(_debug>0) {
+ System.out.printf("simHitPos %s\n",simHitPos.toString());
+ }
+ Hep3Vector vdiffSimHit = VecOp.sub(simHitPos, trkpos);
+ Hep3Vector simHitPos_meas = VecOp.mult(trkToStripRot, vdiffSimHit);
+ Hep3Vector resSimHit_meas = simHitPos_meas; //VecOp.sub(m_meas, simHitPos_meas);
+ file.printStripMeasSimHit(resSimHit_meas.x(),resSimHit_err_meas.x());
+ } else {
+ file.printStripMeasSimHit(-999999.9,-999999.9);
+ }
// print scattering angle