hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -N TwoTrackAnlysis.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TwoTrackAnlysis.java 22 Dec 2012 20:53:33 -0000 1.1
@@ -0,0 +1,301 @@
+package org.lcsim.hps.users.phansson;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.IPlotter;
+import hep.aida.ref.plotter.PlotterRegion;
+import hep.physics.vec.Hep3Vector;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.List;
+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.util.ParticleTypeClassifier;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
+import org.lcsim.hps.recon.tracking.SvtTrackExtrapolator;
+import org.lcsim.hps.recon.vertexing.SimpleVertexer;
+import org.lcsim.hps.recon.vertexing.TwoParticleVertexer;
+import org.lcsim.hps.recon.vertexing.TwoTrackVertexer;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author phansson
+ */
+public class TwoTrackAnlysis extends Driver {
+
+ private FileWriter fileWriter;
+ private PrintWriter printWriter;
+ private String outputNameTextTuple = "twotrackAnlysisTuple.txt";
+ private boolean doPrintBranchInfoLine = true; //firs tline in text file
+ private AIDA aida = AIDA.defaultInstance();
+ private int totalEvents=0;
+ private int totalTwoTrackEvents=0;
+ private int totalMCEvents=0;
+ private int totalTwoTrackMCEvents=0;
+ private boolean hideFrame = false;
+ private String outputPlotFileName;
+ private boolean _debug;
+ private TwoTrackVertexer vertexer = new TwoTrackVertexer();
+ private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
+ private IPlotter _plotterParticleVertex;
+ private IPlotter _plotterTrackVertex;
+ private IHistogram2D _vtxpos_xy;
+ private IHistogram1D _vtxpos_x;
+ private IHistogram1D _vtxpos_y;
+ private IHistogram1D _vtxpos_z;
+ private IHistogram2D _partvtxpos_xy;
+ private IHistogram1D _partvtxpos_x;
+ private IHistogram1D _partvtxpos_y;
+ private IHistogram1D _partvtxpos_z;
+
+
+ public void setDebug(boolean v) {
+ this._debug = v;
+ }
+ public void setOutputPlotFileName(String filename) {
+ outputPlotFileName = filename;
+ }
+
+ public void setHideFrame(boolean hide) {
+ hideFrame = hide;
+ }
+
+
+
+ public TwoTrackAnlysis() {
+ }
+
+
+ @Override
+ public void detectorChanged(Detector detector) {
+
+ try {
+ fileWriter = new FileWriter(outputNameTextTuple);
+ } catch (IOException ex) {
+ Logger.getLogger(TwoTrackAnlysis.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ printWriter = new PrintWriter(fileWriter);
+
+
+ makePlots();
+
+
+ }
+
+
+
+ @Override
+ public void process(EventHeader event) {
+
+
+ totalEvents++;
+
+ List<Track> tracklist = null;
+ if(event.hasCollection(Track.class,"MatchedTracks")) {
+ tracklist = event.get(Track.class, "MatchedTracks");
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": Number of Tracks = " + tracklist.size());
+ }
+ }
+
+ Hep3Vector vtxPosMC = null;
+ MCParticle electron=null;
+ MCParticle positron=null;
+ if(event.hasCollection(MCParticle.class)) {
+ totalMCEvents++;
+ List<MCParticle> mcparticles = event.get(MCParticle.class).get(0);
+ List<MCParticle> fsParticles = HPSMCParticlePlotsDriver.makeGenFSParticleList(mcparticles);
+
+ for(MCParticle part : fsParticles) {
+ if(ParticleTypeClassifier.isElectron(part.getPDGID())) {
+ if(electron==null) {
+ electron = part;
+ } else {
+ if(part.getEnergy()>electron.getEnergy()) {
+ electron = part;
+ }
+ }
+ }
+ if(ParticleTypeClassifier.isPositron(part.getPDGID())) {
+ if(positron==null) {
+ positron = part;
+ } else {
+ if(part.getEnergy()>positron.getEnergy()) {
+ positron = part;
+ }
+ }
+ }
+ if(electron!=null && positron!=null) {
+ this.particleVertexer.setParticle(electron, positron);
+ vtxPosMC = this.particleVertexer.getVertex();
+ if(this._debug) System.out.printf("%s: vtxPosMC=%s\n", this.getClass().getSimpleName(),vtxPosMC.toString());
+ this._partvtxpos_xy.fill(vtxPosMC.x(), vtxPosMC.y());
+ this._partvtxpos_x.fill(vtxPosMC.x());
+ this._partvtxpos_y.fill(vtxPosMC.y());
+ this._partvtxpos_z.fill(vtxPosMC.z());
+ totalTwoTrackMCEvents++;
+ }
+
+ }
+ }
+
+
+ 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());
+ totalTwoTrackEvents++;
+
+ this.fillTextTuple(electron, positron, trk1, trk2, vtxPosMC, vtxPos, event);
+
+
+
+ }
+
+
+
+ @Override
+ public void endOfData() {
+
+ 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);
+
+ if (!"".equals(outputPlotFileName)) {
+ try {
+ aida.saveAs(outputPlotFileName);
+ } catch (IOException ex) {
+ Logger.getLogger(TrigRateDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
+ }
+ }
+
+
+
+ printWriter.close();
+ try {
+ fileWriter.close();
+ } catch (IOException ex) {
+ Logger.getLogger(TwoTrackAnlysis.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+
+ }
+
+
+ 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);
+ _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");
+
+ _plotterParticleVertex.setTitle("MC particle Vertex");
+ _plotterTrackVertex.setTitle("Two Track Vertex");
+
+
+ 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();
+ }
+ }
+
+
+ private boolean isFileEmpty(String fileName) {
+ File f = new File(fileName);
+ 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) {
+ 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+="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";
+ printWriter.println(br_line);
+ doPrintBranchInfoLine = false;
+ }
+
+ //Event info
+ printWriter.format("%5d ",event.getEventNumber());
+ //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]);
+ //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. );
+ //Track vtx
+ printWriter.format("%5.5f %5.5f %5.5f ", vtxPos.x(),vtxPos.y(),vtxPos.z() );
+ //Track at converter
+ this.vertexer.extrapolator().setTrack(trk1);
+ Hep3Vector 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
+ 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
+ printWriter.println();
+
+ }
+
+
+
+
+
+
+
+
+}