 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl on MAIN | |||
| TruthResiduals.java | +20 | -294 | 1.1 -> 1.2 |
| GBLOutputDriver.java | +1 | -2 | 1.3 -> 1.4 |
| +21 | -296 | ||
Cleaning up.
diff -u -r1.1 -r1.2 --- TruthResiduals.java 29 Aug 2013 21:10:43 -0000 1.1 +++ TruthResiduals.java 30 Aug 2013 01:23:22 -0000 1.2 @@ -5,34 +5,22 @@
package org.lcsim.hps.recon.tracking.gbl; import hep.aida.*;
-import hep.aida.ref.histogram.Histogram2D; -import hep.physics.matrix.BasicMatrix; -import hep.physics.matrix.Matrix; -import hep.physics.matrix.MatrixOp; -import hep.physics.matrix.SymmetricMatrix; -import hep.physics.vec.*; -import java.io.FileWriter; -import java.io.PrintWriter;
+import hep.physics.vec.Hep3Vector; +import hep.physics.vec.VecOp;
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.detector.tracker.silicon.ChargeCarrier; import org.lcsim.detector.tracker.silicon.SiSensor;
-import org.lcsim.event.*; -import org.lcsim.fit.helicaltrack.*;
+import org.lcsim.event.MCParticle; +import org.lcsim.event.SimTrackerHit; +import org.lcsim.fit.helicaltrack.HelicalTrackFit; +import org.lcsim.fit.helicaltrack.HelixParamCalculator; +import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.geometry.Detector; import org.lcsim.hps.event.HPSTransformations;
-import org.lcsim.hps.recon.tracking.*; -import org.lcsim.hps.recon.tracking.MaterialSupervisor.DetectorPlane; -import org.lcsim.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume; -import org.lcsim.hps.recon.tracking.MultipleScattering.ScatterPoint; -import org.lcsim.hps.recon.tracking.MultipleScattering.ScatterPoints; -import org.lcsim.recon.tracking.seedtracker.ScatterAngle; -import org.lcsim.recon.tracking.seedtracker.SeedCandidate; -import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.hps.recon.tracking.TrackerHitUtils;
import org.lcsim.util.aida.AIDA; /**
@@ -84,9 +72,6 @@
if(res_truthsimhit == null) makePlots();
- // Use the sim tracker hits to find the truth particle - // this would give me the particle kinematics responsible for the track -
// map layer, mcparticles and sim hits
Map<Integer, List<SimTrackerHit>> simHitsLayerMap = new HashMap<Integer, List<SimTrackerHit> >();
@@ -97,9 +82,6 @@
simHitsLayerMap.put(layer, new ArrayList<SimTrackerHit>());
}
simHitsLayerMap.get(layer).add(sh);
-// if(!simHitsLayerMap.containsKey(layer) || (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength()) ) {
-// simHitsLayerMap.put(layer, sh);
-// }
MCParticle part = sh.getMCParticle();
if(!mcPartSimHitsMap.containsKey(part)) {
@@ -118,18 +100,22 @@
if(simHitsLayer.size()==2 && simHitsLayer.get(0).getMCParticle().equals(simHitsLayer.get(1).getMCParticle())) {
- System.out.printf("%s: There are two simhits in layer %d for MC part %d org %s p %s\n",
- this.getClass().getSimpleName(),layer,
- simHitsLayer.get(0).getMCParticle().getPDGID(),simHitsLayer.get(0).getMCParticle().getOrigin().toString(),simHitsLayer.get(0).getMCParticle().getMomentum().toString());
-
+ if(this._debug>0) {
+ System.out.printf("%s: There are two simhits in layer %d for MC part %d org %s p %s\n",
+ this.getClass().getSimpleName(),layer,
+ simHitsLayer.get(0).getMCParticle().getPDGID(),simHitsLayer.get(0).getMCParticle().getOrigin().toString(),simHitsLayer.get(0).getMCParticle().getMomentum().toString());
+ }
// find out what these hits represent
for(SimTrackerHit simHit : simHitsLayer) {
SiSensor sensor = (SiSensor)simHit.getDetectorElement();
Hep3Vector simHitPosLocal = simHit.getPositionVec();
sensor.getReadoutElectrodes(ChargeCarrier.HOLE).getGlobalToLocal().transform(simHitPosLocal);
- System.out.printf("%s: simHit at local %s and global %s\n",
- this.getClass().getSimpleName(),
- simHitPosLocal.toString(),simHit.getPositionVec().toString());
+ if(this._debug>0) {
+ System.out.printf("%s: simHit at local %s and global %s\n",
+ this.getClass().getSimpleName(),
+ simHitPosLocal.toString(),simHit.getPositionVec().toString());
+
+ }
}
}
@@ -189,204 +175,9 @@
- void processTrack(Track trk, List<MCParticle> mcParticles, List<SimTrackerHit> simTrackerHits) {
-
- SeedTrack st = (SeedTrack)trk;
- SeedCandidate seed = st.getSeedCandidate();
- HelicalTrackFit htf = seed.getHelix();
- List<HelicalTrackHit> hits = seed.getHits();
-
-
- // Use the sim tracker hits to find the truth particle
- // this would give me the particle kinematics responsible for the track
-
- List<MCParticle> mcp_pair = getAprimeDecayProducts(mcParticles);
- if(_debug>0) {
- double invMassTruth = Math.sqrt( Math.pow(mcp_pair.get(0).getEnergy()+mcp_pair.get(1).getEnergy(),2) - VecOp.add(mcp_pair.get(0).getMomentum(), mcp_pair.get(1).getMomentum()).magnitudeSquared() );
- double invMassTruthTrks = getInvMassTracks(getHTF(mcp_pair.get(0)),getHTF(mcp_pair.get(1)));
- System.out.printf("%s: invM = %f\n",this.getClass().getSimpleName(),invMassTruth);
- System.out.printf("%s: invMTracks = %f\n",this.getClass().getSimpleName(),invMassTruthTrks);
- }
-
- MCParticle mcp = getMatchedTruthParticle(trk);
-
- h_mcp_org.fill(mcp.getOriginX(), mcp.getOriginY());
-
- //printMCParticles(mcParticles);
-
- // cross-check
- if(!mcp_pair.contains(mcp)) {
- boolean hasBeamElectronParent = false;
- for(MCParticle parent : mcp.getParents()) {
- if(parent.getGeneratorStatus()!=MCParticle.FINAL_STATE && parent.getPDGID()==11 && parent.getMomentum().y()==0.0 && Math.abs(parent.getMomentum().magnitude() - _beamEnergy)<0.01) {
- hasBeamElectronParent = true;
- }
- }
- if(!hasBeamElectronParent) {
- System.out.printf("%s: the matched MC particle is not an A' daughter and not a the recoil electrons!?\n",this.getClass().getSimpleName());
- System.out.printf("%s: %s %d p %s org %s\n",this.getClass().getSimpleName(),mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString());
- printMCParticles(mcParticles);
- System.exit(1);
- } else {
- if(_debug>0) System.out.printf("%s: the matched MC particle is the recoil electron\n",this.getClass().getSimpleName());
- }
- }
-
- 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, List<SimTrackerHit> > simHitsLayerMap = new HashMap<Integer, List<SimTrackerHit> >();
- for(SimTrackerHit sh : simHits) {
- int layer = sh.getIdentifierFieldValue("layer");
- if(!simHitsLayerMap.containsKey(layer)) {
- //|| (sh.getPathLength() < simHitsLayerMap.get(layer).getPathLength())
- simHitsLayerMap.put(layer, new ArrayList<SimTrackerHit>() );
- }
- simHitsLayerMap.get(layer).add(sh);
- }
-
-
- // Get track paramteres from MC particle
- HelicalTrackFit htfTruth = getHTF(mcp);
-
- if(_debug>0) System.out.printf("%d hits\n",hits.size());
-
-
- int istrip = 0;
- for(int ihit=0;ihit!=hits.size();++ihit) {
-
- HelicalTrackHit hit = hits.get(ihit);
- HelicalTrackCross htc = (HelicalTrackCross) hit;
- List<HelicalTrackStrip> strips = htc.getStrips();
-
- for(HelicalTrackStrip strip : strips) {
-
- if(_debug>0) System.out.printf("%s: layer %d\n",this.getClass().getSimpleName(),strip.layer());
-
-
-
- // Find the sim tracker hit for this layer
- List<SimTrackerHit> simHitsLayer = simHitsLayerMap.get(strip.layer());
- for(SimTrackerHit simHit : simHitsLayer) {
- if(simHitsLayer!=null) {
-
- // Position in tracking coord
- Hep3Vector simHitPosTracking = this._hpstrans.transformVectorToTracking(simHit.getPositionVec());
-
-
- double s_truthSimHit = HelixUtils.PathToXPlane(htfTruth, simHitPosTracking.x(), 0, 0).get(0);
- Hep3Vector trkposTruthSimHit = HelixUtils.PointOnHelix(htfTruth, s_truthSimHit);
- Hep3Vector resTruthSimHit = VecOp.sub(simHitPosTracking,trkposTruthSimHit);
- if(this._debug>0) {
- System.out.printf("TruthSimHit residual %s for layer %d (pathLen=%f)\n",resTruthSimHit.toString(),strip.layer(),s_truthSimHit);
- {
- double ss[] = {0.32,1.0,5.0};
- for(int i=0;i<6;++i) {
- double s_ss = s_truthSimHit + (i<3?ss[i]:-1*ss[i-3]);
- Hep3Vector p1 = HelixUtils.PointOnHelix(htfTruth, s_ss);
- Hep3Vector resP1 = VecOp.sub(simHitPosTracking,p1);
- System.out.printf("TruthSimHit residual %s for layer %d (pathLen=%f)\n",resP1.toString(),strip.layer(),s_ss);
- }
- }
- }
- this.res_truthsimhit.get(strip.layer()).get(0).fill(resTruthSimHit.y());
- this.res_truthsimhit.get(strip.layer()).get(1).fill(resTruthSimHit.z());
- if(strip.layer()==1) this.h_mcp_orgy_resy.fill(mcp.getOriginY(),resTruthSimHit.z());
-
- if(simHitPosTracking.y()>0) {
- if(trk.getCharge()>0) {
- this.res_truthsimhit_top_plus.get(strip.layer()).get(0).fill(resTruthSimHit.y());
- this.res_truthsimhit_top_plus.get(strip.layer()).get(1).fill(resTruthSimHit.z());
- }
- else {
- this.res_truthsimhit_top_minus.get(strip.layer()).get(0).fill(resTruthSimHit.y());
- this.res_truthsimhit_top_minus.get(strip.layer()).get(1).fill(resTruthSimHit.z());
- }
- } else {
- if(trk.getCharge()>0) {
- this.res_truthsimhit_bot_plus.get(strip.layer()).get(0).fill(resTruthSimHit.y());
- this.res_truthsimhit_bot_plus.get(strip.layer()).get(1).fill(resTruthSimHit.z());
- }
- else {
- this.res_truthsimhit_bot_minus.get(strip.layer()).get(0).fill(resTruthSimHit.y());
- this.res_truthsimhit_bot_minus.get(strip.layer()).get(1).fill(resTruthSimHit.z());
- }
- }
- }
- }
-
-
- ++istrip;
-
-
-
- }
-
- }
-
-
-
- }
-
- 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 HelicalTrackFit getHTF(MCParticle mcp) {
@@ -406,71 +197,6 @@
-
- private List<MCParticle> getAprimeDecayProducts(List<MCParticle> mcParticles) {
- List<MCParticle> pair = new ArrayList<MCParticle>();
- for(MCParticle mcp : mcParticles) {
- if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
- boolean hasAprimeParent = false;
- for(MCParticle parent : mcp.getParents()) {
- if(Math.abs(parent.getPDGID())==622) hasAprimeParent = true;
- }
- if(hasAprimeParent) pair.add(mcp);
- }
- if(pair.size()!=2) {
- System.out.printf("%s: ERROR this event has %d mcp with 622 as parent!!?? \n",this.getClass().getSimpleName(),pair.size());
- this.printMCParticles(mcParticles);
- System.exit(1);
- }
- if( Math.abs(pair.get(0).getPDGID()) != 11 || Math.abs(pair.get(1).getPDGID()) != 11 ) {
- System.out.printf("%s: ERROR decay products are not e+e-? \n",this.getClass().getSimpleName());
- this.printMCParticles(mcParticles);
- System.exit(1);
- }
- if(pair.get(0).getPDGID()*pair.get(1).getPDGID() > 0) {
- System.out.printf("%s: ERROR decay products have the same sign? \n",this.getClass().getSimpleName());
- this.printMCParticles(mcParticles);
- System.exit(1);
- }
- return pair;
-
- }
-
- private void printMCParticles(List<MCParticle> mcParticles) {
- System.out.printf("%s: printMCParticles \n",this.getClass().getSimpleName());
- System.out.printf("%s: %d mc particles \n",this.getClass().getSimpleName(),mcParticles.size());
- for(MCParticle mcp : mcParticles) {
- if(mcp.getGeneratorStatus()!=MCParticle.FINAL_STATE) continue;
- System.out.printf("\n%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- mcp.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",mcp.getPDGID(),mcp.getMomentum().toString(),mcp.getOrigin().toString(),
- mcp.getParents().size()>0?"parents:":"");
- for(MCParticle parent : mcp.getParents()) {
- System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- parent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",parent.getPDGID(),parent.getMomentum().toString(),parent.getOrigin().toString(),
- parent.getParents().size()>0?"parents:":"");
- for(MCParticle grparent : parent.getParents()) {
- System.out.printf("%s: (%s) %d p %s org %s %s \n",this.getClass().getSimpleName(),
- grparent.getGeneratorStatus()==MCParticle.FINAL_STATE?"F":"I",grparent.getPDGID(),grparent.getMomentum().toString(),grparent.getOrigin().toString(),
- grparent.getParents().size()>0?"parents:":"");
- }
-
- }
- }
- return;
- }
-
- private double getInvMassTracks(HelicalTrackFit htf1, HelicalTrackFit htf2) {
- double p1 = htf1.p(this._B.magnitude());
- double p2 = htf2.p(this._B.magnitude());
- Hep3Vector p1vec = VecOp.mult(p1, HelixUtils.Direction(htf1, 0));
- Hep3Vector p2vec = VecOp.mult(p2, HelixUtils.Direction(htf2, 0));
- double me = 0.000510998910;
- double E1 = Math.sqrt(p1*p1 + me*me);
- double E2 = Math.sqrt(p2*p2 + me*me);
- //System.out.printf("p1 %f %s E1 %f\n",p1,p1vec.toString(),E1);
- //System.out.printf("p2 %f %s E2 %f\n",p2,p2vec.toString(),E2);
- return Math.sqrt( Math.pow(E1+E2,2) - VecOp.add(p1vec, p2vec).magnitudeSquared() );
- }
private void makePlots() {
diff -u -r1.3 -r1.4 --- GBLOutputDriver.java 29 Aug 2013 21:10:43 -0000 1.3 +++ GBLOutputDriver.java 30 Aug 2013 01:23:22 -0000 1.4 @@ -92,7 +92,7 @@
}
List<SimTrackerHit> simTrackerHits = event.getSimTrackerHits("TrackerHits");
- //truthRes.processSim(mcParticles, simTrackerHits);
+ truthRes.processSim(mcParticles, simTrackerHits);
@@ -127,7 +127,6 @@
gbl.printGBL(trk,mcParticles,simTrackerHits);
- truthRes.processTrack(trk,mcParticles,simTrackerHits);
++iTrack;
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