 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl on MAIN | |||
| TruthResiduals.java | +606 | added 1.1 | |
| GBLFileIO.java | +1 | -1 | 1.5 -> 1.6 |
| GBLOutputDriver.java | +16 | -6 | 1.2 -> 1.3 |
| +623 | -7 | ||
Added 3D path length and truth residual class.
diff -N TruthResiduals.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TruthResiduals.java 29 Aug 2013 21:10:43 -0000 1.1 @@ -0,0 +1,606 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+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 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.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.util.aida.AIDA;
+
+/**
+ *
+ * @author phansson
+ */
+public class TruthResiduals {
+
+ private int _debug;
+ private boolean _hideFrame = false;
+ private Hep3Vector _B;
+ private TrackerHitUtils _trackerHitUtils = new TrackerHitUtils();
+ private Detector detector;
+ private HPSTransformations _hpstrans = new HPSTransformations();
+ private double _beamEnergy = 2.2; //GeV
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ private Map<Integer, List<IHistogram1D> > res_truthsimhit = null;
+ private Map<Integer, List<IHistogram1D> > res_truthsimhit_top_plus = null;
+ private Map<Integer, List<IHistogram1D> > res_truthsimhit_bot_plus = null;
+ private Map<Integer, List<IHistogram1D> > res_truthsimhit_top_minus = null;
+ private Map<Integer, List<IHistogram1D> > res_truthsimhit_bot_minus = null;
+
+
+
+
+ private IHistogram2D h_mcp_org;
+ private IHistogram2D h_mcp_orgy_resy;
+
+
+
+ /*
+ * file name
+ * Bz in Tesla
+ */
+ TruthResiduals(Hep3Vector bfield) {
+ _B = _hpstrans.transformVectorToTracking(bfield);
+ System.out.printf("%s: B field %s\n",this.getClass().getSimpleName(),_B.toString());
+ }
+ public void setDebug(int debug) {
+ _debug = debug;
+ }
+ public void setHideFrame(boolean hide) {
+ _hideFrame = hide;
+ }
+
+
+ void processSim(List<MCParticle> mcParticles, List<SimTrackerHit> simTrackerHits) {
+
+ 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> >();
+ Map<MCParticle, List<SimTrackerHit> > mcPartSimHitsMap = new HashMap<MCParticle, List<SimTrackerHit > >();
+ for(SimTrackerHit sh : simTrackerHits) {
+ int layer = sh.getIdentifierFieldValue("layer");
+ if(!simHitsLayerMap.containsKey(layer)) {
+ 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)) {
+ mcPartSimHitsMap.put(part, new ArrayList<SimTrackerHit>());
+ }
+ mcPartSimHitsMap.get(part).add(sh);
+ }
+
+
+ // Find the particle responsible for the hit in each layer and compute the residual
+
+ for(int layer=1;layer<13;++layer) {
+
+ List<SimTrackerHit> simHitsLayer = simHitsLayerMap.get(layer);
+ if(simHitsLayer != null ) {
+
+ 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());
+
+ // 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());
+ }
+ }
+
+ for(SimTrackerHit simHit : simHitsLayer) {
+ // Position in tracking coord
+ Hep3Vector simHitPosTracking = this._hpstrans.transformVectorToTracking(simHit.getPositionVec());
+ if(_debug>0) {
+ System.out.printf("%s: simHit for layer %d at %s from MC part %d org %s p %s\n",
+ this.getClass().getSimpleName(),layer,simHitPosTracking.toString(),
+ simHit.getMCParticle().getPDGID(),simHit.getMCParticle().getOrigin().toString(),simHit.getMCParticle().getMomentum().toString());
+ }
+
+ // Find the MC particle
+ MCParticle mcp = simHit.getMCParticle();
+
+ // Get track parameters from MC particle
+ HelicalTrackFit htfTruth = getHTF(mcp);
+
+ // Find pathlength to sim hit
+ double s = HelixUtils.PathToXPlane(htfTruth, simHitPosTracking.x(), 0., 0).get(0);
+
+ // Find position on helix
+ Hep3Vector trkpos = HelixUtils.PointOnHelix(htfTruth, s);
+
+ // Calculate residuals
+ Hep3Vector res = VecOp.sub(simHitPosTracking, trkpos);
+
+ // Fill residuals
+ this.res_truthsimhit.get(layer).get(0).fill(res.y());
+ this.res_truthsimhit.get(layer).get(1).fill(res.z());
+ if(simHit.getPositionVec().y()>0) {
+ if(simHit.getMCParticle().getPDGID()<0) {
+ this.res_truthsimhit_top_plus.get(layer).get(0).fill(res.y());
+ this.res_truthsimhit_top_plus.get(layer).get(1).fill(res.z());
+ }
+ else {
+ this.res_truthsimhit_top_minus.get(layer).get(0).fill(res.y());
+ this.res_truthsimhit_top_minus.get(layer).get(1).fill(res.z());
+ }
+ }
+ else {
+ if(simHit.getMCParticle().getPDGID()<0) {
+ this.res_truthsimhit_bot_plus.get(layer).get(0).fill(res.y());
+ this.res_truthsimhit_bot_plus.get(layer).get(1).fill(res.z());
+ }
+ else {
+ this.res_truthsimhit_bot_minus.get(layer).get(0).fill(res.y());
+ this.res_truthsimhit_bot_minus.get(layer).get(1).fill(res.z());
+ }
+ }
+
+ }
+ }
+ }
+
+ }
+
+
+
+ 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) {
+ Hep3Vector org = this._hpstrans.transformVectorToTracking(mcp.getOrigin());
+ Hep3Vector p = this._hpstrans.transformVectorToTracking(mcp.getMomentum());
+ HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1*((int)mcp.getCharge()), -1.0*this._B.z());
+ double par[] = new double[5];
+ par[HelicalTrackFit.dcaIndex] = helixParamCalculator.getDCA();
+ par[HelicalTrackFit.slopeIndex] = helixParamCalculator.getSlopeSZPlane();
+ 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);
+ return htf;
+ }
+
+
+
+
+
+ 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() {
+
+
+ res_truthsimhit = new HashMap<Integer, List<IHistogram1D> >();
+ res_truthsimhit_top_plus = new HashMap<Integer, List<IHistogram1D> >();
+ res_truthsimhit_bot_plus = new HashMap<Integer, List<IHistogram1D> >();
+ res_truthsimhit_top_minus = new HashMap<Integer, List<IHistogram1D> >();
+ res_truthsimhit_bot_minus = new HashMap<Integer, List<IHistogram1D> >();
+
+
+
+
+
+ IHistogramFactory hf = aida.histogramFactory();
+ double min=-1.;
+ double max=1.;
+
+ List<IPlotter> plotter1 = new ArrayList<IPlotter>();
+ List<IPlotter> plotter2 = new ArrayList<IPlotter>();
+ List<IPlotter> plotter3 = new ArrayList<IPlotter>();
+ List<IPlotter> plotter4 = new ArrayList<IPlotter>();
+ List<IPlotter> plotter5 = new ArrayList<IPlotter>();
+
+
+ for(int idir=0;idir<2;++idir) {
+ String dir = idir==0?"x":"y";
+ IPlotter pl1 = af.createPlotterFactory().create(String.format("SimHit-Truth Track Residual %s",dir));
+ pl1.createRegions(3, 4);
+ IPlotter pl2 = af.createPlotterFactory().create(String.format("SimHit-Truth Track Residual %s",dir));
+ pl2.createRegions(3, 4);
+ IPlotter pl3 = af.createPlotterFactory().create(String.format("SimHit-Truth Track Residual %s",dir));
+ pl3.createRegions(3, 4);
+ IPlotter pl4 = af.createPlotterFactory().create(String.format("SimHit-Truth Track Residual %s",dir));
+ pl4.createRegions(3, 4);
+ IPlotter pl5 = af.createPlotterFactory().create(String.format("SimHit-Truth Track Residual %s",dir));
+ pl5.createRegions(3, 4);
+
+ for(int layer=1;layer<13;++layer) {
+
+ if(!res_truthsimhit.containsKey(layer)) {
+ res_truthsimhit.put(layer, new ArrayList<IHistogram1D>());
+ res_truthsimhit_top_plus.put(layer, new ArrayList<IHistogram1D>());
+ res_truthsimhit_bot_plus.put(layer, new ArrayList<IHistogram1D>());
+ res_truthsimhit_top_minus.put(layer, new ArrayList<IHistogram1D>());
+ res_truthsimhit_bot_minus.put(layer, new ArrayList<IHistogram1D>());
+ }
+
+ if(layer<3) {
+ max = 0.07;
+ min = -0.07;
+ } else {
+ max = 0.5 * layer;
+ min = -1.0*max;
+ }
+
+ IHistogram1D h = hf.createHistogram1D(String.format("dres_truthsimhit_layer%d_%s",layer,dir),50, min, max);
+ h.setTitle(String.format("L%d SimHit-Truth Track Residual in %s",layer , dir));
+ res_truthsimhit.get(layer).add(h);
+ pl1.region(layer-1).plot(h);
+
+ h = hf.createHistogram1D(String.format("res_truthsimhit_top_plus_layer%d_%s",layer,dir),50, min, max);
+ h.setTitle(String.format("L%d SimHit-Truth Track (top,q=+1) Residual in %s",layer , dir));
+ res_truthsimhit_top_plus.get(layer).add(h);
+ pl2.region(layer-1).plot(h);
+
+ h = hf.createHistogram1D(String.format("res_truthsimhit_top_minus_layer%d_%s",layer,dir),50, min, max);
+ h.setTitle(String.format("L%d SimHit-Truth Track (top,q=-1) Residual in %s",layer , dir));
+ res_truthsimhit_top_minus.get(layer).add(h);
+ pl3.region(layer-1).plot(h);
+
+ h = hf.createHistogram1D(String.format("res_truthsimhit_bot_minus_layer%d_%s",layer,dir),50, min, max);
+ h.setTitle(String.format("L%d SimHit-Truth Track (bot,q=-1) Residual in %s",layer , dir));
+ res_truthsimhit_bot_minus.get(layer).add(h);
+ pl4.region(layer-1).plot(h);
+
+ h = hf.createHistogram1D(String.format("res_truthsimhit_bot_plus_layer%d_%s",layer,dir),50, min, max);
+ h.setTitle(String.format("L%d SimHit-Truth Track (bot,q=+1) Residual in %s",layer , dir));
+ res_truthsimhit_bot_plus.get(layer).add(h);
+ pl5.region(layer-1).plot(h);
+
+
+ }
+ plotter1.add(pl1);
+ plotter2.add(pl2);
+ plotter3.add(pl3);
+ plotter4.add(pl4);
+ plotter5.add(pl5);
+
+ if(!this._hideFrame) {
+ pl1.show();
+ pl2.show();
+ pl3.show();
+ pl4.show();
+ pl5.show();
+ }
+ else {
+ pl1.hide();
+ pl2.hide();
+ pl3.hide();
+ pl4.hide();
+ pl5.hide();
+ }
+
+ }
+
+
+ this.h_mcp_org = hf.createHistogram2D("MC particle origin", 50, -0.2,0.2,50,-0.2,0.2);
+ IPlotter pl_org = af.createPlotterFactory().create("MC particle origin");
+ pl_org.createRegions(1, 1);
+ pl_org.region(0).plot(h_mcp_org);
+ pl_org.region(0).style().setParameter("hist2DStyle", "colorMap");
+ pl_org.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ if(this._hideFrame) pl_org.hide();
+ else pl_org.show();
+
+ this.h_mcp_orgy_resy = hf.createHistogram2D("MC particle res y vs origin y ", 50, -0.2,0.2,50,-0.1,0.1);
+ IPlotter pl_org_2 = af.createPlotterFactory().create("MC particle res y vs origin y");
+ pl_org_2.createRegions(1, 1);
+ pl_org_2.region(0).plot(this.h_mcp_orgy_resy);
+ pl_org_2.region(0).style().setParameter("hist2DStyle", "colorMap");
+ pl_org_2.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ if(this._hideFrame) pl_org_2.hide();
+ else pl_org_2.show();
+
+
+ }
+
+
+
+
+}
diff -u -r1.5 -r1.6 --- GBLFileIO.java 25 Aug 2013 21:42:19 -0000 1.5 +++ GBLFileIO.java 29 Aug 2013 21:10:43 -0000 1.6 @@ -191,7 +191,7 @@
}
void printStripPathLen3D(double s) {
- addLine(String.format("Strip pathLen 3D %f", s));
+ addLine(String.format("Strip pathLen3D %f", s));
}
void printStereoAngle(double stereoAngle) {
diff -u -r1.2 -r1.3 --- GBLOutputDriver.java 21 Aug 2013 02:44:57 -0000 1.2 +++ GBLOutputDriver.java 29 Aug 2013 21:10:43 -0000 1.3 @@ -27,6 +27,7 @@
String[] detNames = {"Tracker"};
int nevt = 0;
GBLOutput gbl;
+ TruthResiduals truthRes;
private String gblFile = "gblinput.txt";
int totalTracks=0;
int totalTracksProcessed=0;
@@ -64,6 +65,9 @@
gbl.setDebug(_debug);
gbl.setHideFrame(hideFrame);
gbl.buildModel(detector);
+ truthRes = new TruthResiduals(bfield); + truthRes.setDebug(_debug); + truthRes.setHideFrame(hideFrame);
}
@@ -81,6 +85,17 @@
}
}
+
+ List<MCParticle> mcParticles = null;
+ if(event.hasCollection(MCParticle.class,this.MCParticleCollectionName)) {
+ mcParticles = event.get(MCParticle.class,this.MCParticleCollectionName);
+ }
+
+ List<SimTrackerHit> simTrackerHits = event.getSimTrackerHits("TrackerHits");
+ //truthRes.processSim(mcParticles, simTrackerHits);
+
+
+
//gbl.printNewEvent(event.getEventNumber());
gbl.printNewEvent(iEvent);
@@ -109,15 +124,10 @@
gbl.printTrackID(iTrack);
- List<MCParticle> mcParticles = null;
- if(event.hasCollection(MCParticle.class,this.MCParticleCollectionName)) {
- mcParticles = event.get(MCParticle.class,this.MCParticleCollectionName);
- }
- List<SimTrackerHit> simTrackerHits = event.getSimTrackerHits("TrackerHits");
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