Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl on MAIN | |||
GBLFileIO.java | +30 | -1 | 1.2 -> 1.3 |
GBLOutput.java | +120 | -12 | 1.3 -> 1.4 |
GBLOutputDriver.java | +4 | -1 | 1.1 -> 1.2 |
+154 | -14 |
Updated truth matching, included more info in output.
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())); + } +
}
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
diff -u -r1.1 -r1.2 --- GBLOutputDriver.java 14 Jul 2013 16:58:14 -0000 1.1 +++ GBLOutputDriver.java 21 Aug 2013 02:44:57 -0000 1.2 @@ -8,6 +8,7 @@
import java.util.logging.Logger; import org.lcsim.event.EventHeader; import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track; import org.lcsim.geometry.Detector; import org.lcsim.hps.recon.tracking.EventQuality;
@@ -113,7 +114,9 @@
mcParticles = event.get(MCParticle.class,this.MCParticleCollectionName); }
- gbl.printGBL(trk,mcParticles);
+ List<SimTrackerHit> simTrackerHits = event.getSimTrackerHits("TrackerHits"); + + gbl.printGBL(trk,mcParticles,simTrackerHits);
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1