hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
diff -u -r1.5 -r1.6
--- GBLOutput.java 21 Aug 2013 20:13:58 -0000 1.5
+++ GBLOutput.java 22 Aug 2013 00:29:41 -0000 1.6
@@ -12,14 +12,12 @@
import java.io.FileWriter;
import java.io.PrintWriter;
import java.util.ArrayList;
+import java.util.Collections;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
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.event.*;
import org.lcsim.fit.helicaltrack.*;
import org.lcsim.geometry.Detector;
import org.lcsim.hps.event.HPSTransformations;
@@ -120,95 +118,47 @@
file.printClTrackParam(clPar);
- 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);
- }
+
// Use the sim tracker hits to find the truth particle
// this would give me the particle kinematics responsible for the track
- MCParticle mcp = null;
- for(Map.Entry<MCParticle, List<SimTrackerHit>> part : MCParticleHitMap.entrySet()) {
- if(_debug>0) {
- 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
- // Note that charge is screwy because of lcsim wanting positive z bfield in the fit
- 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);
- }
-
- // 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));
- }
+ MCParticle mcp = getMatchedTruthParticle(trk);
-
- // use only the truth particle that should've given the track
- List<SimTrackerHit> simHits = MCParticleHitMap.get(mcp);
+ List<SimTrackerHit> simHits = new ArrayList<SimTrackerHit>();
+ for (SimTrackerHit hit : simTrackerHits) {
+ if(hit.getMCParticle()==mcp) simHits.add(hit);
+ }
if(_debug>0) {
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());
}
+ if(mcp==null) {
+ System.out.printf("%s: no truth particle found!\n",this.getClass().getSimpleName());
+ System.exit(1);
+ }
+
// 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;
- }
+ if(!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength()) ) {
+ simHitsLayerMap.put(layer, sh);
}
- simHitsLayerMap.put(layer, sh);
}
+ // 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);
+
+
// find the projection from the I,J,K to U,V,T coordinates
Hep3Matrix perToClPrj = this.getPerToClPrj(htf);
@@ -270,15 +220,30 @@
// 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());
+ System.out.printf("%s: it as %d mc particles associated with it:\n",this.getClass().getSimpleName(),hit.getMCParticles().size());
+ for (MCParticle particle : hit.getMCParticles()) {
+ System.out.printf("%s: %d p %s \n",this.getClass().getSimpleName(),particle.getPDGID(),particle.getMomentum().toString());
+ }
+ System.out.printf("%s: these are sim hits in the layer map:\n",this.getClass().getSimpleName());
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());
}
+ System.out.printf("%s: these are sim hits in the event:\n",this.getClass().getSimpleName());
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.out.printf("%s: these are all the MC particles in the event:\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");
+ }
//System.exit(1);
}
@@ -433,7 +398,51 @@
}
return fsParticles;
}
-
+
+ MCParticle getMatchedTruthParticle(Track track) {
+ boolean debug = false;
+
+ Map<MCParticle,Integer> particlesOnTrack = new HashMap<MCParticle,Integer>();
+
+ if(debug) System.out.printf("getmatched\n");
+
+ for(TrackerHit hit : track.getTrackerHits()) {
+ List<MCParticle> mcps = ((HelicalTrackHit)hit).getMCParticles();
+ if(mcps==null) {
+ System.out.printf("%s: warning, this hit (layer %d pos=%s) has no mc particles.\n",this.getClass().getSimpleName(),((HelicalTrackHit)hit).Layer(),((HelicalTrackHit)hit).getCorrectedPosition().toString());
+ }
+ else {
+ for(MCParticle mcp : mcps) {
+ if( !particlesOnTrack.containsKey(mcp) ) {
+ particlesOnTrack.put(mcp, 0);
+ }
+ int c = particlesOnTrack.get(mcp);
+ particlesOnTrack.put(mcp, c+1);
+ }
+ }
+ }
+ if(debug) {
+ System.out.printf("Track p=[ %f, %f, %f] \n",track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[1]);
+ System.out.printf("FOund %d particles\n",particlesOnTrack.size());
+ for(Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
+ System.out.printf("%d hits assigned to %d p=%s \n",entry.getValue(),entry.getKey().getPDGID(),entry.getKey().getMomentum().toString());
+ }
+ }
+ Map.Entry<MCParticle,Integer> maxEntry = null;
+ for(Map.Entry<MCParticle,Integer> entry : particlesOnTrack.entrySet()) {
+ if ( maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0 ) maxEntry = entry;
+ //if ( maxEntry != null ) {
+ // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
+ //}
+ //maxEntry = entry;
+ }
+ if(debug) {
+ System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
+ maxEntry.getKey().getPDGID(),maxEntry.getKey().getMomentum().toString(),
+ track.getCharge(),track.getTrackStates().get(0).getMomentum()[0],track.getTrackStates().get(0).getMomentum()[1],track.getTrackStates().get(0).getMomentum()[2]);
+ }
+ return maxEntry.getKey();
+ }
private BasicMatrix getPerParVector(HelicalTrackFit htf) {