hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.2 -r1.3
--- TwoTrackAnlysis.java 8 Jan 2013 00:37:47 -0000 1.2
+++ TwoTrackAnlysis.java 10 Jan 2013 18:35:53 -0000 1.3
@@ -9,6 +9,7 @@
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
+import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
@@ -19,9 +20,12 @@
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;
import org.lcsim.event.util.ParticleTypeClassifier;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
import org.lcsim.geometry.Detector;
import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.evio.TriggerData;
import org.lcsim.hps.recon.ecal.HPSEcalCluster;
import org.lcsim.hps.recon.tracking.SvtTrackExtrapolator;
import org.lcsim.hps.recon.vertexing.SimpleVertexer;
@@ -48,6 +52,7 @@
private boolean hideFrame = false;
private String outputPlotFileName;
private String ecalClusterCollectionName = "EcalClusters";
+ private String stereoHitCollectionName = "RotatedHelicalTrackHits";
private boolean _debug;
private TwoTrackVertexer vertexer = new TwoTrackVertexer();
private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
@@ -109,6 +114,7 @@
@Override
public void process(EventHeader event) {
+ if(this._debug) System.out.println(this.getClass().getSimpleName() + ": processing event "+totalEvents);
totalEvents++;
@@ -121,7 +127,7 @@
}
- if(tracklist.size()!=2) {
+ if(tracklist.size()<2) {
return;
}
@@ -226,6 +232,7 @@
this.fillTextTuple(electron, positron, trk1, trk2, vtxPosMC, vtxPos, clusters, event);
+ if(this._debug) System.out.println(this.getClass().getSimpleName() + ": # two track events so far = "+totalTwoTrackEvents);
}
@@ -238,7 +245,7 @@
System.out.println(this.getClass().getSimpleName() + ": Total Number of Events = "+this.totalEvents);
System.out.println(this.getClass().getSimpleName() + ": Total Number of Two-Track events Processed = "+this.totalTwoTrackEvents);
System.out.println(this.getClass().getSimpleName() + ": Total Number of MCEvents = "+this.totalMCEvents);
- System.out.println(this.getClass().getSimpleName() + ": Total Number of Two-Track events Processed = "+this.totalTwoTrackMCEvents);
+ System.out.println(this.getClass().getSimpleName() + ": Total Number of Two-Track MC events Processed = "+this.totalTwoTrackMCEvents);
if (!"".equals(outputPlotFileName)) {
try {
@@ -261,69 +268,7 @@
}
- private void makePlots() {
- _vtxpos_x = aida.histogram1D("Vertex position X", 100, -1000, 1000);
- _vtxpos_y = aida.histogram1D("Vertex position Y", 100, -200, 200);
- _vtxpos_z = aida.histogram1D("Vertex position Z", 100, -200, 200);
- _vtxpos_xy = aida.histogram2D("Vertex position Y vs X", 50, -200, 200, 50, -200, 200);
- _partvtxpos_x = aida.histogram1D("Particle Vertex position X", 100, -1000, 1000);
- _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);
- _plotterTrackVertex.region(1).plot(_vtxpos_y);
- _plotterTrackVertex.region(2).plot(_vtxpos_z);
- _plotterTrackVertex.region(3).plot(_vtxpos_xy);
- _plotterTrackVertex.region(3).style().setParameter("hist2DStyle", "colorMap");
- _plotterTrackVertex.region(3).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
- _plotterParticleVertex = aida.analysisFactory().createPlotterFactory().create();
- _plotterParticleVertex.createRegions(2,2);
- _plotterParticleVertex.region(0).plot(_partvtxpos_x);
- _plotterParticleVertex.region(1).plot(_partvtxpos_y);
- _plotterParticleVertex.region(2).plot(_partvtxpos_z);
- _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);
- ((PlotterRegion) _plotterParticleVertex.region(i)).getPlot().setAllowPopupMenus(true);
- ((PlotterRegion) _plotterTrackVertex.region(i)).getPlot().setAllowUserInteraction(true);
- ((PlotterRegion) _plotterTrackVertex.region(i)).getPlot().setAllowPopupMenus(true);
-
- }
-
- if(!this.hideFrame) {
- this._plotterParticleVertex.show();
- this._plotterTrackVertex.show();
- this._plotterNhits.show();
- }
- }
-
+
private boolean isFileEmpty(String fileName) {
File f = new File(fileName);
@@ -337,13 +282,21 @@
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: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:";
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk1_strip"+iLayer+"_u/F:";
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk1_strip"+iLayer+"_n/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:";
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk2_strip"+iLayer+"_u/F:";
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk2_strip"+iLayer+"_n/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+="ncl_top/I:ncl_bot/I";
+ br_line+="cl1_E/F:cl1_ix/I:cl1_iy/I:";
+ br_line+="cl2_E/F:cl2_ix/I:cl2_iy/I:";
+ br_line+="cl3_E/F:cl3_ix/I:cl3_iy/I:";
+ br_line+="ncl_top/I:ncl_bot/I:";
+ br_line+="trig_top/I:trig_bot/I";
printWriter.println(br_line);
doPrintBranchInfoLine = false;
}
@@ -360,7 +313,7 @@
for(TrackerHit hit : hitsOnTrack1) {
HelicalTrackHit hth = (HelicalTrackHit) hit;
hits1.put(hth.Layer(), hth);
- System.out.println("hit on layer " + hth.Layer());
+ //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) {
@@ -368,14 +321,44 @@
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<HelicalTrackHit> stereoHits = new ArrayList<HelicalTrackHit>();
+ if(event.hasCollection(HelicalTrackHit.class, stereoHitCollectionName)) {
+ stereoHits = event.get(HelicalTrackHit.class, stereoHitCollectionName);
+ }
+ HashMap<Integer,List<HelicalTrackStrip>> striphits1 = this.getStripHitsMap(hitsOnTrack1);
+ for(int iLayer=1;iLayer<=10;++iLayer) {
+ HelicalTrackStrip strip=null;
+ if(striphits1.containsKey(iLayer)) strip = striphits1.get(iLayer).get(0);
+ if(strip!=null) printWriter.format("%5.5f ", strip.umeas());
+ else printWriter.format("%5.5f ", -99999999.9);
+ }
+ HashMap<Integer,List<HelicalTrackStrip>> allstriphits1 = this.getAllStripHitsMap(stereoHits);
+ for(int iLayer=1;iLayer<=10;++iLayer) {
+ if(striphits1.containsKey(iLayer)) printWriter.format("%5d ", allstriphits1.get(iLayer).size());
+ else printWriter.format("%5d ", -99999999);
+ }
+
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,List<HelicalTrackStrip>> striphits2 = this.getStripHitsMap(hitsOnTrack2);
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);
}
+ for(int iLayer=1;iLayer<=10;++iLayer) {
+ HelicalTrackStrip strip=null;
+ if(striphits2.containsKey(iLayer)) strip = striphits2.get(iLayer).get(0);
+ if(strip!=null) printWriter.format("%5.5f ", strip.umeas());
+ else printWriter.format("%5.5f ", -99999999.9);
+ }
+ HashMap<Integer,List<HelicalTrackStrip>> allstriphits2 = this.getAllStripHitsMap(stereoHits);
+ for(int iLayer=1;iLayer<=10;++iLayer) {
+ if(striphits2.containsKey(iLayer)) printWriter.format("%5d ", striphits2.get(iLayer).size());
+ else printWriter.format("%5d ", -99999999);
+ }
+
//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. );
@@ -389,15 +372,47 @@
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++;
+ for(int i=0;i<3;++i) {
+ if(clusters.size()<=i) {
+ printWriter.format("%5.5f %5d %5d ",-999999.9,-999999,-999999);
+ } else {
+ //for(HPSEcalCluster cl : clusters) {
+ int iy = clusters.get(i).getSeedHit().getIdentifierFieldValue("iy");
+ int ix = clusters.get(i).getSeedHit().getIdentifierFieldValue("ix");
+ double E = clusters.get(i).getEnergy();
+ printWriter.format("%5.5f %5d %5d ",E,ix,iy);
+ if( iy > 0) ncl_t++;
+ else ncl_b++;
+ }
}
- printWriter.format("%5d %5d",ncl_t,ncl_b);
+ printWriter.format("%5d %5d ",ncl_t,ncl_b);
+ printWriter.format("%5d %5d ",ncl_t,ncl_b);
+
+ TriggerData trigger = getTriggerInfo(event);
+ if(trigger==null) printWriter.format("%5d %5d",0,0);
+ else printWriter.format("%5d %5d",trigger.getTopTrig(),trigger.getBotTrig());
printWriter.println();
}
+
+ private TriggerData getTriggerInfo(EventHeader event) {
+ String triggerDecisionCollectionName = "TriggerBank";
+ List<TriggerData> t;
+ if(!event.hasCollection(TriggerData.class, triggerDecisionCollectionName)) {
+ if(_debug) System.out.println( "Event has NO trigger bank");
+ return null;
+ } else {
+ List<TriggerData> triggerDataList = event.get(TriggerData.class, "TriggerBank");
+ if(triggerDataList.isEmpty()) {
+ if(_debug) System.out.println( "Event has trigger bank exists but is empty");
+ return null;
+ } else {
+ return triggerDataList.get(0);
+ }
+ }
+ }
+
private HashMap<Integer,HelicalTrackHit> getHitMap(List<TrackerHit> hits) {
HashMap<Integer,HelicalTrackHit> map = new HashMap<Integer,HelicalTrackHit>();
@@ -408,6 +423,49 @@
return map;
}
+ private HashMap<Integer,HelicalTrackStrip> getStripHitMap(List<TrackerHit> hits) {
+ HashMap<Integer,HelicalTrackStrip> map = new HashMap<Integer,HelicalTrackStrip>();
+ for(TrackerHit hit : hits) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+ HelicalTrackStrip s1 = htc.getStrips().get(0);
+ HelicalTrackStrip s2 = htc.getStrips().get(1);
+ map.put(s1.layer(), s1);
+ map.put(s2.layer(), s2);
+ }
+ return map;
+ }
+
+ private HashMap<Integer,List<HelicalTrackStrip>> getStripHitsMap(List<TrackerHit> hits) {
+ HashMap<Integer,List<HelicalTrackStrip>> map = new HashMap<Integer,List<HelicalTrackStrip>>();
+ for(TrackerHit hit : hits) {
+ HelicalTrackHit hth = (HelicalTrackHit) hit;
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+ HelicalTrackStrip s1 = htc.getStrips().get(0);
+ HelicalTrackStrip s2 = htc.getStrips().get(1);
+ if(!map.containsKey(s1.layer())) map.put(s1.layer(), new ArrayList<HelicalTrackStrip>());
+ if(!map.containsKey(s2.layer())) map.put(s2.layer(), new ArrayList<HelicalTrackStrip>());
+ map.get(s1.layer()).add(s1);
+ map.get(s2.layer()).add(s2);
+ }
+ return map;
+ }
+
+ private HashMap<Integer,List<HelicalTrackStrip>> getAllStripHitsMap(List<HelicalTrackHit> stereoHits) {
+ HashMap<Integer,List<HelicalTrackStrip>> map = new HashMap<Integer,List<HelicalTrackStrip>>();
+ for(HelicalTrackHit hth : stereoHits) {
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+ HelicalTrackStrip s1 = htc.getStrips().get(0);
+ HelicalTrackStrip s2 = htc.getStrips().get(1);
+ if(!map.containsKey(s1.layer())) map.put(s1.layer(), new ArrayList<HelicalTrackStrip>());
+ if(!map.containsKey(s2.layer())) map.put(s2.layer(), new ArrayList<HelicalTrackStrip>());
+ map.get(s1.layer()).add(s1);
+ map.get(s2.layer()).add(s2);
+
+ }
+ return map;
+ }
+
private HelicalTrackHit getHitOnLayer(int layer, List<TrackerHit> hits) {
HelicalTrackHit hitOnLayer = null;
for(TrackerHit hit : hits) {
@@ -419,6 +477,68 @@
return hitOnLayer;
}
+ private void makePlots() {
+ _vtxpos_x = aida.histogram1D("Vertex position X", 100, -1000, 1000);
+ _vtxpos_y = aida.histogram1D("Vertex position Y", 100, -200, 200);
+ _vtxpos_z = aida.histogram1D("Vertex position Z", 100, -200, 200);
+ _vtxpos_xy = aida.histogram2D("Vertex position Y vs X", 50, -200, 200, 50, -200, 200);
+ _partvtxpos_x = aida.histogram1D("Particle Vertex position X", 100, -1000, 1000);
+ _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);
+ _plotterTrackVertex.region(1).plot(_vtxpos_y);
+ _plotterTrackVertex.region(2).plot(_vtxpos_z);
+ _plotterTrackVertex.region(3).plot(_vtxpos_xy);
+ _plotterTrackVertex.region(3).style().setParameter("hist2DStyle", "colorMap");
+ _plotterTrackVertex.region(3).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ _plotterParticleVertex = aida.analysisFactory().createPlotterFactory().create();
+ _plotterParticleVertex.createRegions(2,2);
+ _plotterParticleVertex.region(0).plot(_partvtxpos_x);
+ _plotterParticleVertex.region(1).plot(_partvtxpos_y);
+ _plotterParticleVertex.region(2).plot(_partvtxpos_z);
+ _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);
+ ((PlotterRegion) _plotterParticleVertex.region(i)).getPlot().setAllowPopupMenus(true);
+ ((PlotterRegion) _plotterTrackVertex.region(i)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) _plotterTrackVertex.region(i)).getPlot().setAllowPopupMenus(true);
+
+ }
+
+ if(!this.hideFrame) {
+ this._plotterParticleVertex.show();
+ this._plotterTrackVertex.show();
+ this._plotterNhits.show();
+ }
+ }