 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
| TwoTrackAnlysis.java | +143 | -19 | 1.1 -> 1.2 |
Added simple vtx plots and new variables to ntuple.
diff -u -r1.1 -r1.2 --- TwoTrackAnlysis.java 22 Dec 2012 20:53:33 -0000 1.1 +++ TwoTrackAnlysis.java 8 Jan 2013 00:37:47 -0000 1.2 @@ -9,15 +9,20 @@
import java.io.FileWriter; import java.io.IOException; import java.io.PrintWriter;
+import java.util.HashMap;
import java.util.List;
+import java.util.Map;
import java.util.logging.Level; import java.util.logging.Logger; import org.lcsim.event.EventHeader; import org.lcsim.event.MCParticle; import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
import org.lcsim.geometry.Detector; import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
import org.lcsim.hps.recon.tracking.SvtTrackExtrapolator; import org.lcsim.hps.recon.vertexing.SimpleVertexer; import org.lcsim.hps.recon.vertexing.TwoParticleVertexer;
@@ -42,11 +47,13 @@
private int totalTwoTrackMCEvents=0;
private boolean hideFrame = false;
private String outputPlotFileName;
+ private String ecalClusterCollectionName = "EcalClusters";
private boolean _debug;
private TwoTrackVertexer vertexer = new TwoTrackVertexer();
private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
private IPlotter _plotterParticleVertex;
private IPlotter _plotterTrackVertex;
+ private IPlotter _plotterNhits;
private IHistogram2D _vtxpos_xy;
private IHistogram1D _vtxpos_x;
private IHistogram1D _vtxpos_y;
@@ -55,6 +62,13 @@
private IHistogram1D _partvtxpos_x;
private IHistogram1D _partvtxpos_y;
private IHistogram1D _partvtxpos_z;
+ private IHistogram1D _nhits; + private IHistogram1D _nhits_vtxpospos; + private IHistogram1D _nhits_vtxposneg; + private IHistogram1D _layershit; + private IHistogram1D _layershit_vtxpospos; + private IHistogram1D _layershit_vtxposneg; +
public void setDebug(boolean v) {
@@ -106,6 +120,33 @@
}
}
+
+ if(tracklist.size()!=2) {
+ return;
+ }
+
+
+ Track trk1 = tracklist.get(0);
+ Track trk2 = tracklist.get(1);
+ List<TrackerHit> hitsOnTrack1 = trk1.getTrackerHits();
+ List<TrackerHit> hitsOnTrack2 = trk2.getTrackerHits();
+
+ this.vertexer.setTracks(trk1, trk2);
+ Hep3Vector vtxPos = this.vertexer.getVertex();
+ if(this._debug) System.out.printf("%s: vtxPos=%s\n", this.getClass().getSimpleName(),vtxPos.toString());
+ if(vtxPos.x() != vtxPos.x()) {
+ System.out.printf("%s: vtxPos is NaN -> Skip\n",this.getClass().getSimpleName());
+ return;
+ }
+
+
+
+
+
+ List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.class, ecalClusterCollectionName);
+
+
+
Hep3Vector vtxPosMC = null;
MCParticle electron=null;
MCParticle positron=null;
@@ -148,23 +189,42 @@
}
- if(tracklist.size()!=2) {
- return;
- }
- Track trk1 = tracklist.get(0);
- Track trk2 = tracklist.get(1);
-
- this.vertexer.setTracks(trk1, trk2);
- Hep3Vector vtxPos = this.vertexer.getVertex();
- if(this._debug) System.out.printf("%s: vtxPos=%s\n", this.getClass().getSimpleName(),vtxPos.toString());
+
this._vtxpos_xy.fill(vtxPos.x(), vtxPos.y());
this._vtxpos_x.fill(vtxPos.x());
this._vtxpos_y.fill(vtxPos.y());
this._vtxpos_z.fill(vtxPos.z());
+ this._nhits.fill(hitsOnTrack1.size());
+ this._nhits.fill(hitsOnTrack2.size());
+ if(vtxPos.x()<-350.0) {
+ this._nhits_vtxposneg.fill(hitsOnTrack1.size());
+ this._nhits_vtxposneg.fill(hitsOnTrack2.size());
+ } else {
+ this._nhits_vtxpospos.fill(hitsOnTrack1.size());
+ this._nhits_vtxpospos.fill(hitsOnTrack2.size());
+ }
+ for(TrackerHit hit : hitsOnTrack1) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ int l = (hth.Layer()+1)/2;
+ this._layershit.fill(l);
+ if(vtxPos.x()<-350.0) this._layershit_vtxposneg.fill(l);
+ else this._layershit_vtxpospos.fill(l);
+ }
+
+
+
+// for(HPSEcalCluster cl : clusters) {
+// int[] clusterPosIdx = new int[2];
+// clusterPosIdx[0] = cl.getSeedHit().getIdentifierFieldValue("ix");
+// clusterPosIdx[1] = cl.getSeedHit().getIdentifierFieldValue("iy");
+// }
+
+
+
totalTwoTrackEvents++;
- this.fillTextTuple(electron, positron, trk1, trk2, vtxPosMC, vtxPos, event);
+ this.fillTextTuple(electron, positron, trk1, trk2, vtxPosMC, vtxPos, clusters, event);
@@ -210,6 +270,14 @@
_partvtxpos_y = aida.histogram1D("Particle Vertex position Y", 100, -200, 200);
_partvtxpos_z = aida.histogram1D("Particle Vertex position Z", 100, -200, 200);
_partvtxpos_xy = aida.histogram2D("Particle Vertex position Y vs X", 50, -200, 200, 50, -200, 200);
+ _nhits = aida.histogram1D("# hits", 3, 4, 7);
+ _nhits_vtxpospos = aida.histogram1D("# hits for pos vtx", 3, 4, 7);
+ _nhits_vtxposneg = aida.histogram1D("# hits for neg vtx", 3, 4, 7);
+ _layershit = aida.histogram1D("Layers hit", 5, 1, 6);
+ _layershit_vtxpospos = aida.histogram1D("Layers hit pos vtx", 5, 1, 6);
+ _layershit_vtxposneg = aida.histogram1D("Layers hit neg vtx", 5, 1, 6);
+
+
_plotterTrackVertex = aida.analysisFactory().createPlotterFactory().create();
_plotterTrackVertex.createRegions(2,2);
_plotterTrackVertex.region(0).plot(_vtxpos_x);
@@ -226,10 +294,20 @@
_plotterParticleVertex.region(3).plot(_partvtxpos_xy);
_plotterParticleVertex.region(3).style().setParameter("hist2DStyle", "colorMap");
_plotterParticleVertex.region(3).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ _plotterNhits = aida.analysisFactory().createPlotterFactory().create(); + _plotterNhits.createRegions(1,2); + _plotterNhits.region(0).plot(_nhits); + _plotterNhits.region(0).plot(this._nhits_vtxposneg,"mode=overlay"); + _plotterNhits.region(0).plot(this._nhits_vtxpospos,"mode=overlay"); + _plotterNhits.region(0).style().dataStyle().fillStyle().setVisible(false); + _plotterNhits.region(1).plot(this._layershit); + _plotterNhits.region(1).plot(this._layershit_vtxposneg,"mode=overlay"); + _plotterNhits.region(1).plot(this._layershit_vtxpospos,"mode=overlay"); + _plotterNhits.region(1).style().dataStyle().fillStyle().setVisible(false);
_plotterParticleVertex.setTitle("MC particle Vertex");
_plotterTrackVertex.setTitle("Two Track Vertex");
-
+ _plotterNhits.setTitle("Hits on track");
for(int i=0;i<3;++i) {
((PlotterRegion) _plotterParticleVertex.region(i)).getPlot().setAllowUserInteraction(true);
@@ -242,6 +320,7 @@
if(!this.hideFrame) {
this._plotterParticleVertex.show();
this._plotterTrackVertex.show();
+ this._plotterNhits.show();
}
}
@@ -251,17 +330,20 @@
return f.length() == 0; //return zero also in case file doesn't exist
}
- private void fillTextTuple(MCParticle e, MCParticle p, Track trk1, Track trk2, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, EventHeader event) {
+ private void fillTextTuple(MCParticle e, MCParticle p, Track trk1, Track trk2, Hep3Vector vtxPosParticle, Hep3Vector vtxPos,List<HPSEcalCluster> clusters, EventHeader event) {
if(doPrintBranchInfoLine) {
String br_line = "";
br_line+="evtnr/I:";
br_line+="e_px/F:e_py/F:e_pz/F:p_px/F:p_py/F:p_pz/F:";
- br_line+="trk1_q/I:trk1_chi2/F:trk1_px/F:trk1_py/F:trk1_pz/F:"; - br_line+="trk2_q/I:trk2_chi2/F:trk2_px/F:trk2_py/F:trk2_pz/F:";
+ br_line+="trk1_q/I:trk1_chi2/F:trk1_px/F:trk1_py/F:trk1_pz/F:trk1_nhits/I:"; + for(int iLayer=1;iLayer<=5;++iLayer) br_line+="trk1_hit"+iLayer+"_x/F:"+"trk1_hit"+iLayer+"_y/F:"+"trk1_hit"+iLayer+"_z/F:"; + br_line+="trk2_q/I:trk2_chi2/F:trk2_px/F:trk2_py/F:trk2_pz/F:trk2_nhits/I:"; + for(int iLayer=1;iLayer<=5;++iLayer) br_line+="trk2_hit"+iLayer+"_x/F:"+"trk2_hit"+iLayer+"_y/F:"+"trk2_hit"+iLayer+"_z/F:";
br_line+="vtx_truth_x/F:vtx_truth_y/F:vtx_truth_z/F:";
br_line+="vtx_x/F:vtx_y/F:vtx_z/F:";
br_line+="trk1_conv_x/F:trk1_conv_y/F:trk1_conv_z/F:";
- br_line+="trk2_conv_x/F:trk2_conv_y/F:trk2_conv_z/F";
+ br_line+="trk2_conv_x/F:trk2_conv_y/F:trk2_conv_z/F:";
+ br_line+="ncl_top/I:ncl_bot/I";
printWriter.println(br_line);
doPrintBranchInfoLine = false;
}
@@ -271,10 +353,29 @@
//Truth
if(e!=null && p!=null) printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", e.getPX(),e.getPY(),e.getPZ(), p.getPX(),p.getPY(),p.getPZ() );
else printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", -9999999., -9999999., -9999999., -9999999., -9999999., -9999999. );
- System.out.printf("%s: fill txt tuple 0.1\n",this.getClass().getName());
//Track properties
- printWriter.format("%5d %5.5f %5.5f %5.5f %5.5f ",trk1.getCharge(),trk1.getChi2(), trk1.getTrackStates().get(0).getMomentum()[0],trk1.getTrackStates().get(0).getMomentum()[1],trk1.getTrackStates().get(0).getMomentum()[2]);
- printWriter.format("%5d %5.5f %5.5f %5.5f %5.5f ",trk2.getCharge(),trk2.getChi2(), trk2.getTrackStates().get(0).getMomentum()[0],trk2.getTrackStates().get(0).getMomentum()[1],trk2.getTrackStates().get(0).getMomentum()[2]);
+ List<TrackerHit> hitsOnTrack1 = trk1.getTrackerHits();
+ printWriter.format("%5d %5.5f %5.5f %5.5f %5.5f %5d ",trk1.getCharge(),trk1.getChi2(), trk1.getTrackStates().get(0).getMomentum()[0],trk1.getTrackStates().get(0).getMomentum()[1],trk1.getTrackStates().get(0).getMomentum()[2],hitsOnTrack1.size());
+ HashMap<Integer,HelicalTrackHit> hits1 = this.getHitMap(hitsOnTrack1);
+ for(TrackerHit hit : hitsOnTrack1) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ hits1.put(hth.Layer(), hth);
+ System.out.println("hit on layer " + hth.Layer());
+ }
+ // stupid but I want to keep one line per event so default in case there is not hits in all layers
+ for(int iLayer=0;iLayer<5;++iLayer) {
+ HelicalTrackHit hitOnLayer = hits1.get(iLayer*2+1);// = this.getHitOnLayer(iLayer, hitsOnTrack);
+ if (hitOnLayer != null) printWriter.format("%5.5f %5.5f %5.5f ", hitOnLayer.getPosition()[0],hitOnLayer.getPosition()[1],hitOnLayer.getPosition()[2]);
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ }
+ List<TrackerHit> hitsOnTrack2 = trk2.getTrackerHits();
+ printWriter.format("%5d %5.5f %5.5f %5.5f %5.5f %5d ",trk2.getCharge(),trk2.getChi2(), trk2.getTrackStates().get(0).getMomentum()[0],trk2.getTrackStates().get(0).getMomentum()[1],trk2.getTrackStates().get(0).getMomentum()[2],hitsOnTrack2.size());
+ HashMap<Integer,HelicalTrackHit> hits2 = this.getHitMap(hitsOnTrack2);
+ for(int iLayer=0;iLayer<5;++iLayer) {
+ HelicalTrackHit hitOnLayer = hits2.get(iLayer*2+1);// = this.getHitOnLayer(iLayer, hitsOnTrack);
+ if (hitOnLayer != null) printWriter.format("%5.5f %5.5f %5.5f ", hitOnLayer.getPosition()[0],hitOnLayer.getPosition()[1],hitOnLayer.getPosition()[2]);
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ }
//Particle vtx
if(vtxPosParticle!=null) printWriter.format("%5.5f %5.5f %5.5f ", vtxPosParticle.x(),vtxPosParticle.y(),vtxPosParticle.z() );
else printWriter.format("%5.5f %5.5f %5.5f ", -9999999., -9999999., -9999999. );
@@ -287,13 +388,36 @@
this.vertexer.extrapolator().setTrack(trk2);
posAtConverter = this.vertexer.extrapolator().extrapolateTrack(SvtTrackExtrapolator.HARP_POSITION);
printWriter.format("%5.5f %5.5f %5.5f ", posAtConverter.z(),posAtConverter.x(),posAtConverter.y()); //note rotation from JLab->tracking
+ int ncl_t=0; int ncl_b=0;
+ for(HPSEcalCluster cl : clusters) {
+ if(cl.getSeedHit().getIdentifierFieldValue("iy") > 0) ncl_t++;
+ else ncl_b++;
+ }
+ printWriter.format("%5d %5d",ncl_t,ncl_b);
printWriter.println();
}
+ private HashMap<Integer,HelicalTrackHit> getHitMap(List<TrackerHit> hits) {
+ HashMap<Integer,HelicalTrackHit> map = new HashMap<Integer,HelicalTrackHit>();
+ for(TrackerHit hit : hits) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ map.put(hth.Layer(), hth);
+ }
+ return map;
+ }
-
+ private HelicalTrackHit getHitOnLayer(int layer, List<TrackerHit> hits) {
+ HelicalTrackHit hitOnLayer = null;
+ for(TrackerHit hit : hits) {
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+ if( htc.Layer() == layer) {
+ hitOnLayer = htc;
+ }
+ }
+ return hitOnLayer;
+ }
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