Commit in hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl on MAIN | |||
TruthResiduals.java | +82 | -49 | 1.3 -> 1.4 |
Cleaning up.
diff -u -r1.3 -r1.4 --- TruthResiduals.java 17 Sep 2013 00:52:52 -0000 1.3 +++ TruthResiduals.java 10 Oct 2013 23:43:46 -0000 1.4 @@ -4,7 +4,11 @@
*/ package org.lcsim.hps.recon.tracking.gbl;
-import hep.aida.*;
+import hep.aida.IAnalysisFactory; +import hep.aida.IHistogram1D; +import hep.aida.IHistogram2D; +import hep.aida.IHistogramFactory; +import hep.aida.IPlotter;
import hep.physics.vec.Hep3Vector; import hep.physics.vec.VecOp;
@@ -13,17 +17,12 @@
import java.util.List; import java.util.Map;
-import org.lcsim.detector.tracker.silicon.ChargeCarrier; -import org.lcsim.detector.tracker.silicon.SiSensor;
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.SvtTrackExtrapolator;
import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.hps.recon.tracking.TrackerHitUtils;
import org.lcsim.util.aida.AIDA; /**
@@ -35,10 +34,7 @@
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;
@@ -46,12 +42,9 @@
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;
+ private IHistogram2D trkpos_y_vs_x; + private boolean firstWeirdTrack = true;
@@ -80,6 +73,14 @@
Map<Integer, List<SimTrackerHit>> simHitsLayerMap = new HashMap<Integer, List<SimTrackerHit> >(); Map<MCParticle, List<SimTrackerHit> > mcPartSimHitsMap = new HashMap<MCParticle, List<SimTrackerHit > >(); for(SimTrackerHit sh : simTrackerHits) {
+ Hep3Vector shpos = this._hpstrans.transformVectorToTracking(sh.getPositionVec()); + if(Math.abs(shpos.x()) < 50.0) { + System.out.printf("%s: Weird hit at %s (%s) in layer %d for MC part %d org %s p %s\n", + this.getClass().getSimpleName(),sh.getPositionVec().toString(),shpos.toString(),sh.getIdentifierFieldValue("layer"), + sh.getMCParticle().getPDGID(),sh.getMCParticle().getOrigin().toString(),sh.getMCParticle().getMomentum().toString()); + System.exit(1); + } +
int layer = sh.getIdentifierFieldValue("layer"); if(!simHitsLayerMap.containsKey(layer)) { simHitsLayerMap.put(layer, new ArrayList<SimTrackerHit>());
@@ -101,56 +102,53 @@
// Find the particle responsible for the hit in each layer and compute the residual for(int layer=1;layer<13;++layer) {
+ //System.out.printf("layer %d: \n",layer);
List<SimTrackerHit> simHitsLayer = simHitsLayerMap.get(layer);
+ +
if(simHitsLayer != null ) {
- - if(simHitsLayer.size()==2 && simHitsLayer.get(0).getMCParticle().equals(simHitsLayer.get(1).getMCParticle())) { - - 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); - if(this._debug>0) { - System.out.printf("%s: simHit at local %s and global %s\n", - this.getClass().getSimpleName(), - simHitPosLocal.toString(),simHit.getPositionVec().toString()); - - } - } - } -
+ + if(simHitsLayer.size()==2) continue; +
for(SimTrackerHit simHit : simHitsLayer) {
+ + // Find the MC particle + MCParticle mcp = simHit.getMCParticle(); + + if(mcp.getMomentum().magnitude()<0.5) continue; +
// 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(),
+ System.out.printf("%s: simHit for layer %d at %s (%s) from MC part %d org %s p %s\n", + this.getClass().getSimpleName(),layer,simHit.getPositionVec().toString(),simHitPosTracking.toString(),
simHit.getMCParticle().getPDGID(),simHit.getMCParticle().getOrigin().toString(),simHit.getMCParticle().getMomentum().toString());
+ + if(simHitPosTracking.x()<50.) System.exit(1); +
}
- // Find the MC particle - MCParticle mcp = simHit.getMCParticle();
+
// Get track parameters from MC particle //HelicalTrackFit htfTruth = TrackUtils.getHTF(mcp, -1*this._B.z()); HelicalTrackFit htfTruth = TrackUtils.getHTF(mcp, -1*this._B.z());
- // 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);
+ + SvtTrackExtrapolator extrapol = new SvtTrackExtrapolator(htfTruth.parameters()); + Hep3Vector trkposExtraPolator = extrapol.extrapolateTrack(simHitPosTracking.x()); + //System.out.printf("trkposextrapol (det) %s\n",trkposExtraPolator.toString()); + + trkposExtraPolator = this._hpstrans.transformVectorToTracking(trkposExtraPolator);
// Calculate residuals
- Hep3Vector res = VecOp.sub(simHitPosTracking, trkpos);
+ Hep3Vector res = VecOp.sub(simHitPosTracking, trkposExtraPolator);
+ //System.out.printf("trkpos %s trkposextrapol %s res %s\n",trkpos.toString(),trkposExtraPolator.toString(),res.toString()); +
// Fill residuals this.res_truthsimhit.get(layer).get(0).fill(res.y()); this.res_truthsimhit.get(layer).get(1).fill(res.z());
@@ -174,6 +172,28 @@
this.res_truthsimhit_bot_minus.get(layer).get(1).fill(res.z()); } }
+ + if(layer == 1 && res.y() > 0.1 && this.firstWeirdTrack) { + double dx = 1.0; + double xpos = mcp.getOriginZ(); + while(xpos< 100.) { + xpos += dx; + trkposExtraPolator = this._hpstrans.transformVectorToTracking(extrapol.extrapolateTrack(xpos)); + double ypos = trkposExtraPolator.y(); + trkpos_y_vs_x.fill(xpos,ypos); + } + + int idummy = 0; + while(idummy<2) { + trkpos_y_vs_x.fill(simHitPosTracking.x(),simHitPosTracking.y()); + idummy++; + //System.out.printf("weird simhit res pos %s \n", simHitPosTracking.toString()); + } + + this.firstWeirdTrack = false; + + } +
} }
@@ -233,9 +253,12 @@
res_truthsimhit_bot_minus.put(layer, new ArrayList<IHistogram1D>()); }
- if(layer<3) { - max = 0.07; - min = -0.07;
+ if(layer<2) { + max = 0.05;//0.07; + min = -0.05;//-0.07; + } else if(layer<3) { + max = 0.3;//0.07; + min = -0.3;//-0.07;
} else { max = 0.5 * layer; min = -1.0*max;
@@ -302,6 +325,16 @@
else pl_org.show();
+ trkpos_y_vs_x = hf.createHistogram2D("Track pos y vs x", 300, -150,150,100,-4,4); + IPlotter pl_pos_y_vs_x = af.createPlotterFactory().create("Track pos y vs x"); + pl_pos_y_vs_x.createRegions(1, 1); + pl_pos_y_vs_x.region(0).plot(trkpos_y_vs_x); + pl_pos_y_vs_x.region(0).style().setParameter("hist2DStyle", "colorMap"); + pl_pos_y_vs_x.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow"); + if(this._hideFrame) pl_pos_y_vs_x.hide(); + else pl_pos_y_vs_x.show(); + +
}
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