hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.13 -r1.14
--- TwoTrackAnlysis.java 18 Mar 2013 22:33:44 -0000 1.13
+++ TwoTrackAnlysis.java 27 Mar 2013 01:15:25 -0000 1.14
@@ -9,10 +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;
+import java.util.*;
import java.util.logging.Level;
import java.util.logging.Logger;
import org.lcsim.detector.IDetectorElement;
@@ -24,7 +21,6 @@
import org.lcsim.event.util.ParticleTypeClassifier;
import org.lcsim.fit.helicaltrack.*;
import org.lcsim.geometry.Detector;
-import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
import org.lcsim.hps.analysis.ecal.HPSMCParticlePlotsDriver;
import org.lcsim.hps.evio.TriggerData;
import org.lcsim.hps.recon.ecal.HPSEcalCluster;
@@ -32,10 +28,8 @@
import org.lcsim.hps.recon.tracking.SvtTrackExtrapolator;
import org.lcsim.hps.recon.tracking.SvtUtils;
import org.lcsim.hps.recon.tracking.TrackUtils;
-import org.lcsim.hps.recon.vertexing.SimpleVertexer;
import org.lcsim.hps.recon.vertexing.TwoParticleVertexer;
import org.lcsim.hps.recon.vertexing.TwoTrackVertexer;
-import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.Driver;
@@ -70,14 +64,25 @@
private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
private IPlotter _plotterParticleVertex;
private IPlotter _plotterTrackVertex;
+ private IPlotter _plotterTrackMult;
private IHistogram1D _vtxpos_x;
private IHistogram1D _vtxpos_y;
private IHistogram1D _vtxpos_z;
private IHistogram1D _partvtxpos_x;
private IHistogram1D _partvtxpos_y;
private IHistogram1D _partvtxpos_z;
+ private IHistogram2D _ntrks_px;
-
+ public class CmpTrack implements Comparable<CmpTrack> {
+ private Track _track;
+ public CmpTrack(Track track) {
+ _track = track;
+ }
+ @Override
+ public int compareTo(CmpTrack t) {
+ return _track.getTrackStates().get(0).getMomentum()[0]>t._track.getTrackStates().get(0).getMomentum()[0]?0:1;
+ }
+ }
public void setDebug(boolean v) {
@@ -138,7 +143,7 @@
if(event.hasCollection(Track.class,trackCollectionName)) {
tracklist = event.get(Track.class, trackCollectionName);
if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": Number of Tracks = " + tracklist.size());
+ System.out.println(this.getClass().getSimpleName() + ": Number of Tracks = " + tracklist.size() + " in event " + event.getEventNumber());
}
} else {
if(_debug) {
@@ -147,43 +152,51 @@
return;
}
- if(_debug) {
- for(int i=0;i<tracklist.size();++i) {
- Track trk = tracklist.get(i);
-
- System.out.printf("%s: trk momentum (%.3f,%.3f,%.3f) chi2=%.3f\n",this.getClass().getSimpleName(),trk.getTrackStates().get(0).getMomentum()[0],trk.getTrackStates().get(0).getMomentum()[1],trk.getTrackStates().get(0).getMomentum()[2],trk.getChi2());
-
- }
- }
- if(tracklist.size()==0) {
- if(_debug) {
- System.out.printf("%s: event %d has only %d tracks \n",this.getClass().getSimpleName(),event.getEventNumber(),tracklist.size());
- }
- return;
- }
-
+ //if(tracklist.size()<2) {
+ //if(_debug)
+ // System.out.printf("%s: event %d has only %d tracks \n",this.getClass().getSimpleName(),event.getEventNumber(),tracklist.size());
+ // System.exit(1);
+ // return;
+ //}
+ TrackUtils trkutil = new TrackUtils();
+ ArrayList<CmpTrack> tracks = new ArrayList<CmpTrack>();
for(int i=0;i<tracklist.size();++i) {
Track trk = tracklist.get(i);
- if(trk.getTrackStates().get(0).getMomentum()[0]<0.1) {
- if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": trk failed momentum cut " + event.getEventNumber() + "\n" + trk.toString());
- }
- return;
+ if(_debug)
+ System.out.printf("%s: trk momentum (%.3f,%.3f,%.3f) chi2=%.3f\n",this.getClass().getSimpleName(),trk.getTrackStates().get(0).getMomentum()[0],trk.getTrackStates().get(0).getMomentum()[1],trk.getTrackStates().get(0).getMomentum()[2],trk.getChi2());
+ if(trk.getTrackStates().get(0).getMomentum()[0]>0.) {
+ tracks.add(new CmpTrack(trk));
+ } else {
+ if(_debug) System.out.println(this.getClass().getSimpleName() + ": trk failed momentum cut " + event.getEventNumber() + "\n" + trk.toString());
}
}
- List<Track> tracks = new ArrayList<Track>();
+ Collections.sort(tracks);
+
+ if(_debug) {
+ CmpTrack trk_prev = null;
+ for(CmpTrack trk:tracks) {
+ if(trk_prev!=null) {
+ if(trk_prev._track.getTrackStates().get(0).getMomentum()[0]<trk._track.getTrackStates().get(0).getMomentum()[0]) {
+ throw new RuntimeException(String.format("%s: ERROR prev px=%f trk=%f",this.getClass().getSimpleName(),trk_prev._track.getTrackStates().get(0).getMomentum()[0],trk._track.getTrackStates().get(0).getMomentum()[0]));
+ }
+ }
+ trk_prev = trk;
+
+ }
+ }
+
Hep3Vector vtxPos = null;
Hep3Vector vtxPosFringe = null;
- if(tracklist.size()>1) {
- Track trk1 = tracklist.get(0);
- Track trk2 = tracklist.get(1);
+ if(tracks.size()>1) {
+ Track trk1 = tracks.get(0)._track;
+ Track trk2 = tracks.get(1)._track;
- this.vertexer.setTracks(tracklist.get(0), trk2);
+ this.vertexer.setTracks(trk1, trk2);
vtxPos = this.vertexer.getVertex();
if(this._debug)
System.out.printf("%s: vtxPos=%s\n", this.getClass().getSimpleName(),vtxPos.toString());
@@ -266,7 +279,7 @@
totalTwoTrackEvents++;
- this.fillTextTuple(electron, positron, tracklist, vtxPosMC, vtxPos, vtxPosFringe, clusters, event);
+ this.fillTextTuple(electron, positron, tracks, vtxPosMC, vtxPos, vtxPosFringe, clusters, event);
if(this._debug) System.out.println(this.getClass().getSimpleName() + ": # two track events so far = "+totalTwoTrackEvents);
@@ -311,10 +324,11 @@
return f.length() == 0; //return zero also in case file doesn't exist
}
- private void fillTextTuple(MCParticle e, MCParticle p, List<Track> tracks, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, Hep3Vector vtxPosFr, List<HPSEcalCluster> clusters, EventHeader event) {
+ private void fillTextTuple(MCParticle e, MCParticle p, List<CmpTrack> tracks, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, Hep3Vector vtxPosFr, List<HPSEcalCluster> clusters, EventHeader event) {
if(doPrintBranchInfoLine) {
String br_line = "";
br_line+="evtnr/I:";
+ br_line+="ntrks_top/I:ntrks_bot/I:ntrks100_top/I:ntrks100_bot/I:ntrks200_top/I:ntrks200_bot/I:ntrks300_top/I:ntrks300_bot/I:ntrks400_top/I:ntrks400_bot/I:ntrks500_top/I:ntrks500_bot/I:";
br_line+="e_px/F:e_py/F:e_pz/F:p_px/F:p_py/F:p_pz/F:";
br_line+="trk1_d0/F:trk1_phi0/F:trk1_R/F:trk1_z0/F:trk1_slope/F:";
br_line+="trk1_q/I:trk1_chi2/F:trk1_px/F:trk1_py/F:trk1_pz/F:trk1_nhits/I:";
@@ -329,7 +343,7 @@
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk1_strip"+iLayer+"_u/F:";
br_line+="trk1_conv_x/F:trk1_conv_y/F:trk1_conv_z/F:";
br_line+="trk1_fr_conv_x/F:trk1_fr_conv_y/F:trk1_fr_conv_z/F:";
-
+ //
br_line+="trk2_d0/F:trk2_phi0/F:trk2_R/F:trk2_z0/F:trk2_slope/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:";
@@ -343,7 +357,35 @@
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk2_strip"+iLayer+"_u/F:";
br_line+="trk2_conv_x/F:trk2_conv_y/F:trk2_conv_z/F:";
br_line+="trk2_fr_conv_x/F:trk2_fr_conv_y/F:trk2_fr_conv_z/F:";
-
+ //
+ br_line+="trk3_d0/F:trk3_phi0/F:trk3_R/F:trk3_z0/F:trk3_slope/F:";
+ br_line+="trk3_q/I:trk3_chi2/F:trk3_px/F:trk3_py/F:trk3_pz/F:trk3_nhits/I:";
+ for(int iLayer=1;iLayer<=5;++iLayer) br_line+="trk3_hit"+iLayer+"_x/F:"+"trk3_hit"+iLayer+"_y/F:"+"trk3_hit"+iLayer+"_z/F:";
+ for(int iLayer=1;iLayer<=5;++iLayer) {
+ br_line+="trk3_res"+iLayer+"_x/F:"+"trk3_res"+iLayer+"_y/F:"+"trk3_res"+iLayer+"_z/F:";
+ br_line+="trk3_eres"+iLayer+"_x/F:"+"trk3_eres"+iLayer+"_y/F:"+"trk3_eres"+iLayer+"_z/F:";
+ br_line+="trk3_drdphi"+iLayer+"/F:"+"trk3_msdrphi"+iLayer+"/F:";
+ br_line+="trk3_dz"+iLayer+"/F:"+"trk3_msdz"+iLayer+"/F:";
+ //br_line+="trk3_ures"+iLayer+"/F:"+"trk3_ureserr"+iLayer+"/F:";
+ }
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk3_strip"+iLayer+"_u/F:";
+ br_line+="trk3_conv_x/F:trk3_conv_y/F:trk3_conv_z/F:";
+ br_line+="trk3_fr_conv_x/F:trk3_fr_conv_y/F:trk3_fr_conv_z/F:";
+ //
+ br_line+="trk4_d0/F:trk4_phi0/F:trk4_R/F:trk4_z0/F:trk4_slope/F:";
+ br_line+="trk4_q/I:trk4_chi2/F:trk4_px/F:trk4_py/F:trk4_pz/F:trk4_nhits/I:";
+ for(int iLayer=1;iLayer<=5;++iLayer) br_line+="trk4_hit"+iLayer+"_x/F:"+"trk4_hit"+iLayer+"_y/F:"+"trk4_hit"+iLayer+"_z/F:";
+ for(int iLayer=1;iLayer<=5;++iLayer) {
+ br_line+="trk4_res"+iLayer+"_x/F:"+"trk4_res"+iLayer+"_y/F:"+"trk4_res"+iLayer+"_z/F:";
+ br_line+="trk4_eres"+iLayer+"_x/F:"+"trk4_eres"+iLayer+"_y/F:"+"trk4_eres"+iLayer+"_z/F:";
+ br_line+="trk4_drdphi"+iLayer+"/F:"+"trk4_msdrphi"+iLayer+"/F:";
+ br_line+="trk4_dz"+iLayer+"/F:"+"trk4_msdz"+iLayer+"/F:";
+ //br_line+="trk4_ures"+iLayer+"/F:"+"trk4_ureserr"+iLayer+"/F:";
+ }
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk4_strip"+iLayer+"_u/F:";
+ br_line+="trk4_conv_x/F:trk4_conv_y/F:trk4_conv_z/F:";
+ br_line+="trk4_fr_conv_x/F:trk4_fr_conv_y/F:trk4_fr_conv_z/F:";
+ //
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="top_strip"+iLayer+"_n/F:";
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="bot_strip"+iLayer+"_n/F:";
for(int iLayer=1;iLayer<=10;iLayer+=2) br_line+="top_stereo"+iLayer+"_n/F:";
@@ -351,9 +393,9 @@
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+="vtx_fr_x/F:vtx_fr_y/F:vtx_fr_z/F:";
- 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+="cl1_E/F:cl1_ix/I:cl1_iy/I:cl1_x/F:cl1_y/F:";
+ br_line+="cl2_E/F:cl2_ix/I:cl2_iy/I:cl2_x/F:cl2_y/F:";
+ br_line+="cl3_E/F:cl3_ix/I:cl3_iy/I:cl3_x/F:cl3_y/F:";
br_line+="ncl_top/I:ncl_bot/I:";
br_line+="trig_top/I:trig_bot/I";
printWriter.println(br_line);
@@ -362,14 +404,24 @@
//Event info
printWriter.format("%5d ",event.getEventNumber());
+ if(tracks.size()>0) {
+ for(int icut=0;icut<=5;++icut) {
+ int ntrks[] = getNtracks(tracks,icut*0.1);
+ printWriter.format("%5d %5d ",ntrks[0],ntrks[1]);
+ this._ntrks_px.fill(0.1*icut, ntrks[0]+ntrks[1]);
+ }
+ } else {
+ printWriter.format("%5d %5d %5d %5d %5d %5d %5d %5d %5d %5d %5d %5d ", -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999, -9999999 );
+ }
+
//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. );
- for (int itrk=0;itrk<2;itrk++) {
+ for (int itrk=0;itrk<4;itrk++) {
Track trk1 = null;
- if(tracks.size()>itrk) trk1 = tracks.get(itrk);
+ if(tracks.size()>itrk) trk1 = tracks.get(itrk)._track;
if(trk1!=null) {
SeedTrack st1 = (SeedTrack) trk1;
@@ -414,12 +466,13 @@
}
//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
+ Hep3Vector posAtConverter = this.vertexer.extrapolator().extrapolateTrack(SvtTrackExtrapolator.HARP_POSITION);
+ if(posAtConverter!=null) printWriter.format("%5.5f %5.5f %5.5f ", posAtConverter.z(),posAtConverter.x(),posAtConverter.y()); //note rotation from JLab->tracking
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9,-9999999.9,-9999999.9);
HPSTrack hpstrk1 = new HPSTrack(helix1);
Hep3Vector posAtConverterFringe1 = hpstrk1.getPositionAtZMap(100., SvtTrackExtrapolator.HARP_POSITION, 5.0)[0];
- printWriter.format("%5.5f %5.5f %5.5f ", posAtConverterFringe1.z(),posAtConverterFringe1.x(),posAtConverterFringe1.y()); //note rotation from JLab->tracking
-
+ if (posAtConverterFringe1!=null) printWriter.format("%5.5f %5.5f %5.5f ", posAtConverterFringe1.z(),posAtConverterFringe1.x(),posAtConverterFringe1.y()); //note rotation from JLab->tracking
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9,-9999999.9,-9999999.9);
}
else {
@@ -455,8 +508,14 @@
HashMap<Integer,List<SiTrackerHitStrip1D>> allstriphits = this.getAllStripHitsMap(event,true);
//System.out.printf("%s: %d strip hits in event\n",this.getClass().getSimpleName(),allstriphits.size());
for(int iLayer=1;iLayer<=10;++iLayer) {
- if(allstriphits.containsKey(iLayer)) printWriter.format("%5d ", allstriphits.get(iLayer).size());
- else printWriter.format("%5d ", -99999999);
+ if(allstriphits.containsKey(iLayer)) {
+ printWriter.format("%5d ", allstriphits.get(iLayer).size());
+ //System.out.printf("%s: layer %d has %d strip hits\n",this.getClass().getSimpleName(),iLayer,allstriphits.get(iLayer).size());
+ }
+ else {
+ printWriter.format("%5d ", -99999999);
+ //System.out.printf("%s: layer %d has 0 strip hits\n",this.getClass().getSimpleName(),iLayer);
+ }
}
allstriphits = this.getAllStripHitsMap(event,false);
for(int iLayer=1;iLayer<=10;++iLayer) {
@@ -496,16 +555,17 @@
int ncl_t=0; int ncl_b=0;
for(int i=0;i<3;++i) {
if(clusters==null) {
- printWriter.format("%5.5f %5d %5d ",-999999.9,-999999,-999999);
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f ",-999999.9,-999999,-999999,-999999.,-999999.);
}
else if(clusters.size()<=i) {
- printWriter.format("%5.5f %5d %5d ",-999999.9,-999999,-999999);
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f ",-999999.9,-999999,-999999,-999999.,-999999.);
} else {
//for(HPSEcalCluster cl : clusters) {
int iy = clusters.get(i).getSeedHit().getIdentifierFieldValue("iy");
int ix = clusters.get(i).getSeedHit().getIdentifierFieldValue("ix");
+ double pos[] = clusters.get(i).getPosition();
double E = clusters.get(i).getEnergy();
- printWriter.format("%5.5f %5d %5d ",E,ix,iy);
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f ",E,ix,iy,pos[0],pos[1]);
if( iy > 0) ncl_t++;
else ncl_b++;
}
@@ -618,6 +678,37 @@
}
return map;
}
+
+ private int[] getNtracks(List<CmpTrack> tracks) {
+ return this.getNtracks(tracks, 0.);
+ }
+ private int[] getNtracks(List<CmpTrack> tracks, double min_px) {
+ //System.out.printf("%s: getNtracks for %d tracks with min_px=%.3f \n",this.getClass().getSimpleName(),tracks.size(),min_px);
+ int n[] = {0,0};
+ for(CmpTrack track : tracks) {
+ if(track._track.getTrackStates().get(0).getMomentum()[0]<min_px) {
+ continue;
+ }
+ //System.out.printf("%s: track had enough px=%f\n",this.getClass().getSimpleName(),track._track.getTrackStates().get(0).getMomentum()[0]);
+
+ List<TrackerHit> hitsOnTrack = track._track.getTrackerHits();
+ for(TrackerHit hit : hitsOnTrack) {
+ double y = hit.getPosition()[1];
+ if(y>0) {
+ //System.out.printf("%s: this track (chi2=%f) is a top track\n",this.getClass().getSimpleName(),track._track.getChi2());
+ n[0]++;
+ break;
+ } else if(y<0) {
+ //System.out.printf("%s: this track (chi2=%f) is a bot track\n",this.getClass().getSimpleName(),track._track.getChi2());
+ n[1]++;
+ break;
+ }
+ }
+ }
+ //System.out.printf("%s: found %d top and %d bot tracks\n",this.getClass().getSimpleName(),n[0],n[1]);
+
+ return n;
+ }
@@ -628,6 +719,7 @@
_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);
+ _ntrks_px = aida.histogram2D("Track multiplicity vs px cut", 6, 0, 0.6,8,0,8);
_plotterTrackVertex = aida.analysisFactory().createPlotterFactory().create();
_plotterTrackVertex.createRegions(2,2);
@@ -639,20 +731,29 @@
_plotterParticleVertex.region(0).plot(_partvtxpos_x);
_plotterParticleVertex.region(1).plot(_partvtxpos_y);
_plotterParticleVertex.region(2).plot(_partvtxpos_z);
+ _plotterTrackMult = aida.analysisFactory().createPlotterFactory().create();
+ _plotterTrackMult.createRegions(1,1);
+ _plotterTrackMult.region(0).plot(_ntrks_px);
_plotterParticleVertex.setTitle("MC particle Vertex");
_plotterTrackVertex.setTitle("Two Track Vertex");
+ _plotterTrackMult.setTitle("Track multiplicity");
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(i==0) {
+ ((PlotterRegion) _plotterTrackMult.region(i)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) _plotterTrackMult.region(i)).getPlot().setAllowPopupMenus(true);
+ }
}
if(!this.hideFrame) {
this._plotterParticleVertex.show();
this._plotterTrackVertex.show();
+ this._plotterTrackMult.show();
}
}