hps-java/src/main/java/org/lcsim/hps/recon/tracking/gbl
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();
+
+
}