hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.17 -r1.18
--- TwoTrackAnlysis.java 23 Apr 2013 00:07:37 -0000 1.17
+++ TwoTrackAnlysis.java 26 Jun 2013 22:06:33 -0000 1.18
@@ -4,7 +4,9 @@
import hep.aida.IHistogram2D;
import hep.aida.IPlotter;
import hep.aida.ref.plotter.PlotterRegion;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.SpacePoint;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
@@ -23,12 +25,14 @@
import org.lcsim.hps.evio.TriggerData;
import org.lcsim.hps.recon.ecal.HPSEcalCluster;
import org.lcsim.hps.recon.tracking.*;
+import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.vertexing.TwoParticleVertexer;
import org.lcsim.hps.recon.vertexing.TwoTrackVertexer;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.swim.Helix;
/**
*
@@ -36,8 +40,8 @@
*/
public class TwoTrackAnlysis extends Driver {
- private FileWriter fileWriter;
- private PrintWriter printWriter;
+ private FileWriter fileWriter = null;
+ private PrintWriter printWriter = null;
private String outputNameTextTuple = "twotrackAnlysisTuple.txt";
private String trackCollectionName = "MatchedTracks";
private boolean doPrintBranchInfoLine = true; //firs tline in text file
@@ -60,14 +64,40 @@
private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
private IPlotter _plotterParticleVertex;
private IPlotter _plotterTrackVertex;
+ private IPlotter _plotterTrackVertexNonBend;
private IPlotter _plotterTrackMult;
+ private IPlotter _plotterTrackAtConv;
private IHistogram1D _vtxpos_x;
private IHistogram1D _vtxpos_y;
private IHistogram1D _vtxpos_z;
+ private IHistogram1D _vtxposnonb_x;
+ private IHistogram1D _vtxposnonb_y;
+ private IHistogram1D _vtxposnonb_dy;
+ private IHistogram1D _vtxposnonb_z;
private IHistogram1D _partvtxpos_x;
private IHistogram1D _partvtxpos_y;
private IHistogram1D _partvtxpos_z;
+ private IHistogram1D _vtxposnonb_xAtZ0;
+ private IHistogram1D _vtxposnonb_zAtTarget;
+ private IHistogram1D _vtxposnonb_angle1;
+ private IHistogram1D _vtxposnonb_angle2;
private IHistogram2D _ntrks_px;
+ private IHistogram1D _trk_y_at_conv_top_pos;
+ private IHistogram1D _trk_z_at_conv_top_pos;
+ private IHistogram1D _trk_y_at_conv_top_pos_fr;
+ private IHistogram1D _trk_z_at_conv_top_pos_fr;
+ private IHistogram1D _trk_y_at_conv_bot_pos;
+ private IHistogram1D _trk_z_at_conv_bot_pos;
+ private IHistogram1D _trk_y_at_conv_bot_pos_fr;
+ private IHistogram1D _trk_z_at_conv_bot_pos_fr;
+ private IHistogram1D _trk_y_at_conv_top_neg;
+ private IHistogram1D _trk_z_at_conv_top_neg;
+ private IHistogram1D _trk_y_at_conv_top_neg_fr;
+ private IHistogram1D _trk_z_at_conv_top_neg_fr;
+ private IHistogram1D _trk_y_at_conv_bot_neg;
+ private IHistogram1D _trk_z_at_conv_bot_neg;
+ private IHistogram1D _trk_y_at_conv_bot_neg_fr;
+ private IHistogram1D _trk_z_at_conv_bot_neg_fr;
public class CmpTrack implements Comparable<CmpTrack> {
private Track _track;
@@ -110,14 +140,13 @@
@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);
+ createWriter();
+
+ if(printWriter==null)
+ System.out.println("printWriter is null");
+
+ fillTextTupleBranches();
-
makePlots();
@@ -140,8 +169,8 @@
List<Track> tracklist = null;
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() + " in event " + event.getEventNumber());
+ if(_debug)
+ System.out.println(this.getClass().getSimpleName() + ": Number of Tracks = " + tracklist.size() + " in event " + event.getEventNumber());
} else {
if(_debug) {
System.out.println(this.getClass().getSimpleName() + ": No track collection in event " + event.getEventNumber());
@@ -161,6 +190,7 @@
for(int i=0;i<tracklist.size();++i) {
Track trk = tracklist.get(i);
if(TrackUtils.isGoodTrack(trk, tracklist, EventQuality.Quality.MEDIUM)) {
+ //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(_debug) {
int cuts = TrackUtils.passTrackSelections(trk, tracklist, EventQuality.Quality.MEDIUM);
System.out.printf("%s: track cuts: \n%s\n",this.getClass().getSimpleName(),EventQuality.instance().print(cuts));
@@ -198,6 +228,10 @@
Hep3Vector vtxPos = null;
Hep3Vector vtxPosFringe = null;
+ Hep3Vector vtxPosNonBend = null;
+
+
+
if(tracks.size()>1) {
@@ -220,15 +254,45 @@
System.out.printf("%s: vtxPosFringe is NaN -> Skip\n",this.getClass().getSimpleName());
vtxPos = null;
}
+
+ StraightLineTrack[] slts = this.getSLTs(trk1, trk2, false);
+ double zAtCross = this.getCrossingS(trk1, trk2, false);
+ double[] xyAtZ1 = slts[0].calculateXYAtZ(zAtCross);
+ double[] xyAtZ2 = slts[1].calculateXYAtZ(zAtCross);
+
+ Hep3Vector[] vtxNonBend = {new BasicHep3Vector(xyAtZ1[0],xyAtZ1[1],zAtCross),new BasicHep3Vector(xyAtZ2[0],xyAtZ2[1],zAtCross)};
+ //System.out.printf("%s: vtxNonBend=%s\n", this.getClass().getSimpleName(),vtxNonBend.toString());
+
+ this._vtxposnonb_x.fill(vtxNonBend[0].x());
+ this._vtxposnonb_y.fill(vtxNonBend[0].y());
+ this._vtxposnonb_y.fill(vtxNonBend[1].y());
+ this._vtxposnonb_dy.fill(vtxNonBend[0].y()-vtxNonBend[1].y());
+ this._vtxposnonb_z.fill(vtxNonBend[0].z());
+
+ this._vtxposnonb_xAtZ0.fill(slts[0].calculateXYAtZ(0.)[0]);
+ this._vtxposnonb_zAtTarget.fill(slts[0].getYZAtX(SvtTrackExtrapolator.HARP_POSITION)[1]);
+ this._vtxposnonb_angle1.fill(Math.atan(slts[0].dzdx()));
+ this._vtxposnonb_angle2.fill(Math.atan(slts[1].dzdx()));
+
+ vtxPosNonBend = vtxNonBend[0];
+
+ //slts = this.getSLTs(trk1, trk2, true);
+ //zAtCross = this.getCrossingS(trk1, trk2, true);
+ //xyAtZ1 = slts[0].calculateXYAtZ(zAtCross);
+ //xyAtZ2 = slts[1].calculateXYAtZ(zAtCross);
+ //Hep3Vector[] vtxNonBendFringe = {new BasicHep3Vector(xyAtZ1[0],xyAtZ1[1],zAtCross), new BasicHep3Vector(xyAtZ2[0],xyAtZ2[1],zAtCross)};
+ //System.out.printf("%s: vtxNonBendFridge=%s\n", this.getClass().getSimpleName(),vtxNonBendFringe.toString());
+
}
if(vtxPos!=null) {
this._vtxpos_x.fill(vtxPos.x());
this._vtxpos_y.fill(vtxPos.y());
- this._vtxpos_z.fill(vtxPos.z());
+ this._vtxpos_z.fill(vtxPos.z());
}
+
List<HPSEcalCluster> clusters = new ArrayList<HPSEcalCluster>();
if(!event.hasCollection(HPSEcalCluster.class, ecalClusterCollectionName)) {
@@ -287,13 +351,13 @@
totalTwoTrackEvents++;
try {
- this.fillTextTuple(electron, positron, tracks, vtxPosMC, vtxPos, vtxPosFringe, clusters, event);
+ this.fillTextTuple(electron, positron, tracks, vtxPosMC, vtxPos, vtxPosFringe, vtxPosNonBend, clusters, event);
} catch (IOException ex) {
Logger.getLogger(TwoTrackAnlysis.class.getName()).log(Level.SEVERE, null, ex);
}
if(this._debug) System.out.println(this.getClass().getSimpleName() + ": # two track events so far = "+totalTwoTrackEvents);
-
+ //System.exit(0);
}
@@ -315,8 +379,14 @@
}
}
+ if(printWriter==null) {
+ //Fill branches to text file (this is because detectorChanged() is called only if there are events to process
+ createWriter();
+ fillTextTupleBranches();
+ }
+
-
+ System.out.println("printWriter out: " + printWriter.toString());
printWriter.close();
try {
fileWriter.close();
@@ -324,7 +394,7 @@
Logger.getLogger(TwoTrackAnlysis.class.getName()).log(Level.SEVERE, null, ex);
}
-
+
}
@@ -335,90 +405,70 @@
return f.length() == 0; //return zero also in case file doesn't exist
}
- private void fillTextTuple(MCParticle e, MCParticle p, List<CmpTrack> tracks, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, Hep3Vector vtxPosFr, List<HPSEcalCluster> clusters, EventHeader event) throws IOException {
+
+ private void fillTextTupleBranches() {
+
+ if(!doPrintBranchInfoLine) {
+ throw new RuntimeException("Trying to fill tuple branches with flag being set to false?!");
+ }
+
+ 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:";
+
+ for(int itrk=1;itrk<=4;++itrk) {
+ String trk_str = String.format("trk%d", itrk);
+
+ br_line+=""+trk_str+"_d0/F:"+trk_str+"_phi0/F:"+trk_str+"_R/F:"+trk_str+"_z0/F:"+trk_str+"_slope/F:";
+ br_line+=""+trk_str+"_q/I:"+trk_str+"_chi2/F:"+trk_str+"_px/F:"+trk_str+"_py/F:"+trk_str+"_pz/F:"+trk_str+"_nhits/I:";
+ for(int iLayer=1;iLayer<=5;++iLayer) br_line+=""+trk_str+"_hit"+iLayer+"_x/F:"+""+trk_str+"_hit"+iLayer+"_y/F:"+""+trk_str+"_hit"+iLayer+"_z/F:";
+ for(int iLayer=1;iLayer<=5;++iLayer) {
+ br_line+=""+trk_str+"_res"+iLayer+"_y/F:"+""+trk_str+"_res"+iLayer+"_z/F:";
+ br_line+=""+trk_str+"_eres"+iLayer+"_y/F:"+""+trk_str+"_eres"+iLayer+"_z/F:";
+ br_line+=""+trk_str+"_drdphi"+iLayer+"/F:"+""+trk_str+"_msdrphi"+iLayer+"/F:";
+ br_line+=""+trk_str+"_dz"+iLayer+"/F:"+""+trk_str+"_msdz"+iLayer+"/F:";
+ //br_line+=""+trk_str+"_ures"+iLayer+"/F:"+""+trk_str+"_ureserr"+iLayer+"/F:";
+ }
+ for(int iLayer=1;iLayer<=10;++iLayer) br_line+=""+trk_str+"_strip"+iLayer+"_u/F:"+""+trk_str+"_strip"+iLayer+"_time/F:"+""+trk_str+"_strip"+iLayer+"_E/F:";
+
+ br_line+=""+trk_str+"_conv_y/F:"+trk_str+"_conv_z/F:";
+ br_line+=""+trk_str+"_fr_conv_x/F:"+trk_str+"_fr_conv_y/F:"+trk_str+"_fr_conv_z/F:";
+ br_line+=""+trk_str+"_target_y/F:"+trk_str+"_target_z/F:";
+ br_line+=""+trk_str+"_fr_target_x/F:"+trk_str+"_fr_target_y/F:"+trk_str+"_fr_target_z/F:";
+ br_line+=""+trk_str+"_ecal_y/F:"+trk_str+"_ecal_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:";
+ for(int iLayer=1;iLayer<=10;iLayer+=2) br_line+="bot_stereo"+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+="vtx_fr_x/F:vtx_fr_y/F:vtx_fr_z/F:";
+ br_line+="vtx_nonbend_x/F:vtx_nonbend_y/F:vtx_nonbend_z/F:";
+ br_line+="cl1_E/F:cl1_ix/I:cl1_iy/I:cl1_x/F:cl1_y/F:cl1_n/I:";
+ br_line+="cl2_E/F:cl2_ix/I:cl2_iy/I:cl2_x/F:cl2_y/F:cl2_n/I:";
+ br_line+="cl3_E/F:cl3_ix/I:cl3_iy/I:cl3_x/F:cl3_y/F:cl3_n/I:";
+ br_line+="ncl_top/I:ncl_bot/I:";
+ br_line+="trig_top/I:trig_bot/I";
+
+ if(printWriter==null) {
+ System.out.println("hmm 1");
+ }
+ printWriter.println(br_line);
+
+ doPrintBranchInfoLine = false;
+
+ }
+
+
+
+
+ private void fillTextTuple(MCParticle e, MCParticle p, List<CmpTrack> tracks, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, Hep3Vector vtxPosFr, Hep3Vector vtxPosNonBend, List<HPSEcalCluster> clusters, EventHeader event) throws IOException {
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:";
- for(int itrk=1;itrk<=4;++itrk) {
- String trk_str = String.format("trk%d", itrk);
-
- br_line+=""+trk_str+"_d0/F:"+trk_str+"_phi0/F:"+trk_str+"_R/F:"+trk_str+"_z0/F:"+trk_str+"_slope/F:";
- br_line+=""+trk_str+"_q/I:"+trk_str+"_chi2/F:"+trk_str+"_px/F:"+trk_str+"_py/F:"+trk_str+"_pz/F:"+trk_str+"_nhits/I:";
- for(int iLayer=1;iLayer<=5;++iLayer) br_line+=""+trk_str+"_hit"+iLayer+"_x/F:"+""+trk_str+"_hit"+iLayer+"_y/F:"+""+trk_str+"_hit"+iLayer+"_z/F:";
- for(int iLayer=1;iLayer<=5;++iLayer) {
- br_line+=""+trk_str+"_res"+iLayer+"_y/F:"+""+trk_str+"_res"+iLayer+"_z/F:";
- br_line+=""+trk_str+"_eres"+iLayer+"_y/F:"+""+trk_str+"_eres"+iLayer+"_z/F:";
- br_line+=""+trk_str+"_drdphi"+iLayer+"/F:"+""+trk_str+"_msdrphi"+iLayer+"/F:";
- br_line+=""+trk_str+"_dz"+iLayer+"/F:"+""+trk_str+"_msdz"+iLayer+"/F:";
- //br_line+=""+trk_str+"_ures"+iLayer+"/F:"+""+trk_str+"_ureserr"+iLayer+"/F:";
- }
- for(int iLayer=1;iLayer<=10;++iLayer) br_line+=""+trk_str+"_strip"+iLayer+"_u/F:"+""+trk_str+"_strip"+iLayer+"_time/F:"+""+trk_str+"_strip"+iLayer+"_E/F:";
- br_line+=""+trk_str+"_conv_y/F:"+trk_str+"_conv_z/F:";
- br_line+=""+trk_str+"_fr_conv_x/F:"+trk_str+"_fr_conv_y/F:"+trk_str+"_fr_conv_z/F:";
- br_line+=""+trk_str+"_ecal_y/F:"+trk_str+"_ecal_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:";
-// for(int iLayer=1;iLayer<=5;++iLayer) {
-// br_line+="trk2_res"+iLayer+"_x/F:"+"trk2_res"+iLayer+"_y/F:"+"trk2_res"+iLayer+"_z/F:";
-// br_line+="trk2_eres"+iLayer+"_x/F:"+"trk2_eres"+iLayer+"_y/F:"+"trk2_eres"+iLayer+"_z/F:";
-// br_line+="trk2_drdphi"+iLayer+"/F:"+"trk2_msdrphi"+iLayer+"/F:";
-// br_line+="trk2_dz"+iLayer+"/F:"+"trk2_msdz"+iLayer+"/F:";
-// //br_line+="trk2_ures"+iLayer+"/F:"+"trk2_ureserr"+iLayer+"/F:";
-// }
-// for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk2_strip"+iLayer+"_u/F:"+"trk2_strip"+iLayer+"_time/F:"+"trk2_strip"+iLayer+"_E/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+="trk2_ecal_x/F:trk2_ecal_y/F:trk2_ecal_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+="trk3_ecal_x/F:trk3_ecal_y/F:trk3_ecal_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:";
-// br_line+="trk4_ecal_x/F:trk4_ecal_y/F:trk4_ecal_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:";
- for(int iLayer=1;iLayer<=10;iLayer+=2) br_line+="bot_stereo"+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+="vtx_fr_x/F:vtx_fr_y/F:vtx_fr_z/F:";
- 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);
- doPrintBranchInfoLine = false;
+ throw new RuntimeException("Need to fill tuple branches first!?");
}
//Event info
@@ -476,6 +526,7 @@
printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
}
}
+
HashMap<Integer,List<HelicalTrackStrip>> striphits1 = this.getStripHitsMap(hitsOnTrack1);
for(int iLayer=1;iLayer<=10;++iLayer) {
HelicalTrackStrip strip=null;
@@ -487,21 +538,71 @@
printWriter.format("%5.5f %5.5f %5.5f ", -99999999.9, -99999999.9, -99999999.9);
}
}
+
//Track at converter
this.vertexer.extrapolator().setTrack(trk1);
Hep3Vector posAtConverter = this.vertexer.extrapolator().extrapolateTrack(SvtTrackExtrapolator.HARP_POSITION);
- if(posAtConverter!=null) printWriter.format("%5.5f %5.5f ", posAtConverter.x(),posAtConverter.y()); //note rotation from JLab->tracking
+ if(posAtConverter!=null) printWriter.format("%5.5f %5.5f ", posAtConverter.x(),posAtConverter.y()); //note rotation from JLab->tracking
else printWriter.format("%5.5f %5.5f ", -9999999.9,-9999999.9);
HPSTrack hpstrk1 = new HPSTrack(helix1);
Hep3Vector posAtConverterFringe1 = hpstrk1.getPositionAtZMap(100., SvtTrackExtrapolator.HARP_POSITION, 5.0)[0];
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);
+
+ Hep3Vector posAtNomTarget1 = this.vertexer.extrapolator().extrapolateTrack(0);
+ if(posAtNomTarget1!=null) printWriter.format("%5.5f %5.5f ", posAtNomTarget1.x(),posAtNomTarget1.y()); //note rotation from JLab->tracking
+ else printWriter.format("%5.5f %5.5f ", -9999999.9,-9999999.9);
+
+ Hep3Vector posAtNomTargetFringe1 = hpstrk1.getPositionAtZMap(100., 0.0, 5.0)[0];
+ if (posAtNomTargetFringe1!=null) printWriter.format("%5.5f %5.5f %5.5f ", posAtNomTargetFringe1.z(),posAtNomTargetFringe1.x(),posAtNomTargetFringe1.y()); //note rotation from JLab->tracking
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9,-9999999.9,-9999999.9);
Hep3Vector posAtECal = this.vertexer.extrapolator().extrapolateTrack(SvtTrackExtrapolator.ECAL_FACE);
if(posAtECal!=null && !Double.isNaN(posAtECal.x()) && !Double.isNaN(posAtECal.y())) {
printWriter.format("%5.5f %5.5f ",posAtECal.x(),posAtECal.y()); //note rotation from JLab->tracking
}
else printWriter.format("%5.5f %5.5f ",-9999999.9,-9999999.9);
-
+
+ if(posAtConverter!=null) {
+ if(TrackUtils.isTopTrack(trk1, 4)) {
+ if(trk1.getCharge()>0) {
+ this._trk_y_at_conv_top_pos.fill(posAtConverter.x());
+ this._trk_z_at_conv_top_pos.fill(posAtConverter.y());
+ } else {
+ this._trk_y_at_conv_top_neg.fill(posAtConverter.x());
+ this._trk_z_at_conv_top_neg.fill(posAtConverter.y());
+ }
+ } else {
+ if(trk1.getCharge()>0) {
+ this._trk_y_at_conv_bot_pos.fill(posAtConverter.x());
+ this._trk_z_at_conv_bot_pos.fill(posAtConverter.y());
+ } else {
+ this._trk_y_at_conv_bot_neg.fill(posAtConverter.x());
+ this._trk_z_at_conv_bot_neg.fill(posAtConverter.y());
+ }
+ }
+ }
+ if(posAtConverterFringe1!=null) {
+ if(TrackUtils.isTopTrack(trk1, 4)) {
+ if(trk1.getCharge()>0) {
+ this._trk_y_at_conv_top_pos_fr.fill(posAtConverterFringe1.x());
+ this._trk_z_at_conv_top_pos_fr.fill(posAtConverterFringe1.y());
+ } else {
+ this._trk_y_at_conv_top_neg_fr.fill(posAtConverterFringe1.x());
+ this._trk_z_at_conv_top_neg_fr.fill(posAtConverterFringe1.y());
+ }
+ } else {
+ if(trk1.getCharge()>0) {
+ this._trk_y_at_conv_bot_pos_fr.fill(posAtConverterFringe1.x());
+ this._trk_z_at_conv_bot_pos_fr.fill(posAtConverterFringe1.y());
+ } else {
+ this._trk_y_at_conv_bot_neg_fr.fill(posAtConverterFringe1.x());
+ this._trk_z_at_conv_bot_neg_fr.fill(posAtConverterFringe1.y());
+ }
+ }
+ }
+
+
+
}
else {
@@ -519,10 +620,13 @@
for(int iLayer=1;iLayer<=10;++iLayer) {
printWriter.format("%5.5f %5.5f %5.5f ", -99999999.9, -99999999.9, -99999999.9);
}
+
printWriter.format("%5.5f %5.5f ",-9999999.9,-9999999.9); //note rotation from JLab->tracking
printWriter.format("%5.5f %5.5f %5.5f ",-9999999.9,-9999999.9,-9999999.9); //note rotation from JLab->tracking
printWriter.format("%5.5f %5.5f ",-9999999.9,-9999999.9); //note rotation from JLab->tracking
-
+ printWriter.format("%5.5f %5.5f %5.5f ",-9999999.9,-9999999.9,-9999999.9); //note rotation from JLab->tracking
+ printWriter.format("%5.5f %5.5f ",-9999999.9,-9999999.9); //note rotation from JLab->tracking
+
}
}
//printWriter.format("\n%s\n","X1");
@@ -579,22 +683,25 @@
//Track vtx
if(vtxPos!=null) printWriter.format("%5.5f %5.5f %5.5f ", vtxPos.x(),vtxPos.y(),vtxPos.z() );
else printWriter.format("%5.5f %5.5f %5.5f ", -9999999., -9999999., -9999999. );
- if(vtxPos!=null) printWriter.format("%5.5f %5.5f %5.5f ", vtxPosFr.x(),vtxPosFr.y(),vtxPosFr.z() );
+ if(vtxPosFr!=null) printWriter.format("%5.5f %5.5f %5.5f ", vtxPosFr.x(),vtxPosFr.y(),vtxPosFr.z() );
+ else printWriter.format("%5.5f %5.5f %5.5f ", -9999999., -9999999., -9999999. );
+ if(vtxPosNonBend!=null) printWriter.format("%5.5f %5.5f %5.5f ", vtxPosNonBend.x(),vtxPosNonBend.y(),vtxPosNonBend.z() );
else printWriter.format("%5.5f %5.5f %5.5f ", -9999999., -9999999., -9999999. );
int ncl_t=0; int ncl_b=0;
for(int i=0;i<3;++i) {
if(clusters==null) {
- printWriter.format("%5.5f %5d %5d %5.5f %5.5f ",-999999.9,-999999,-999999,-999999.,-999999.);
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f %5d ",-999999.9,-999999,-999999,-999999.,-999999.,-999999);
}
else if(clusters.size()<=i) {
- printWriter.format("%5.5f %5d %5d %5.5f %5.5f ",-999999.9,-999999,-999999,-999999.,-999999.);
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f %5d ",-999999.9,-999999,-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 %5.5f %5.5f ",E,ix,iy,pos[0],pos[1]);
+ int clsize = clusters.get(i).getSize();
+ printWriter.format("%5.5f %5d %5d %5.5f %5.5f %5d ",E,ix,iy,pos[0],pos[1],clsize);
if( iy > 0) ncl_t++;
else ncl_b++;
}
@@ -744,19 +851,58 @@
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);
- _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);
+ _vtxpos_x = aida.histogram1D("Vertex position X", 100, -800, -500);
+ _vtxpos_y = aida.histogram1D("Vertex position Y", 100, 0, 40);
+ _vtxpos_z = aida.histogram1D("Vertex position Z", 100, -4, 4);
+ _vtxposnonb_x = aida.histogram1D("Vertex position non-bend X", 100, -800, -500);
+ _vtxposnonb_y = aida.histogram1D("Vertex position non-bend Y", 100, 0, 40);
+ _vtxposnonb_dy = aida.histogram1D("Vertex position non-bend Ytrk1-Ytrk2", 100, -20, 20);
+ _vtxposnonb_z = aida.histogram1D("Vertex position non-bend Z", 100, -4, 4);
+ _partvtxpos_x = aida.histogram1D("Particle Vertex position X", 100, -1000, 0);
+ _partvtxpos_y = aida.histogram1D("Particle Vertex position Y", 100, -20, 60);
+ _partvtxpos_z = aida.histogram1D("Particle Vertex position Z", 100, -5, 5);
+ _vtxposnonb_xAtZ0 = aida.histogram1D("SLT position non-bend X at Z=0", 100, -1000, 0);
+ _vtxposnonb_zAtTarget = aida.histogram1D("SLT position non-bend Z at target", 100, -5, 5);
+ _vtxposnonb_angle1 = aida.histogram1D("SLT thetay non-bend (trk1)", 100, -0.05, 0.05);
+ _vtxposnonb_angle2 = aida.histogram1D("SLT thetay non-bend (trk2)", 100, -0.05, 0.05);
+
_ntrks_px = aida.histogram2D("Track multiplicity vs px cut", 6, 0, 0.6,8,0,8);
+ _trk_y_at_conv_top_pos = aida.histogram1D("Track y @ converter top +", 100, -30,60);
+ _trk_z_at_conv_top_pos = aida.histogram1D("Track z @ converter top +", 100, -20,20);
+ _trk_y_at_conv_top_pos_fr = aida.histogram1D("Track y @ converter top + (fr)", 100, -30, 60);
+ _trk_z_at_conv_top_pos_fr = aida.histogram1D("Track z @ converter top + (fr)", 100, -20, 20);
+ _trk_y_at_conv_bot_pos = aida.histogram1D("Track y @ converter bot +", 100, -30,60);
+ _trk_z_at_conv_bot_pos = aida.histogram1D("Track z @ converter bot +", 100, -20,20);
+ _trk_y_at_conv_bot_pos_fr = aida.histogram1D("Track y @ converter bot + (fr)", 100, -30, 30);
+ _trk_z_at_conv_bot_pos_fr = aida.histogram1D("Track z @ converter bot + (fr)", 100, -20, 20);
+
+ _trk_y_at_conv_top_neg = aida.histogram1D("Track y @ converter top -", 100, -30,60);
+ _trk_z_at_conv_top_neg = aida.histogram1D("Track z @ converter top -", 100, -20,20);
+ _trk_y_at_conv_top_neg_fr = aida.histogram1D("Track y @ converter top - (fr)", 100, -30, 60);
+ _trk_z_at_conv_top_neg_fr = aida.histogram1D("Track z @ converter top - (fr)", 100, -20, 20);
+ _trk_y_at_conv_bot_neg = aida.histogram1D("Track y @ converter bot -", 100, -30,60);
+ _trk_z_at_conv_bot_neg = aida.histogram1D("Track z @ converter bot -", 100, -20,20);
+ _trk_y_at_conv_bot_neg_fr = aida.histogram1D("Track y @ converter bot - (fr)", 100, -30, 30);
+ _trk_z_at_conv_bot_neg_fr = aida.histogram1D("Track z @ converter bot - (fr)", 100, -20, 20);
+
+
+
_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);
+ _plotterTrackVertexNonBend = aida.analysisFactory().createPlotterFactory().create();
+ _plotterTrackVertexNonBend.createRegions(2,4);
+ _plotterTrackVertexNonBend.region(0).plot(_vtxposnonb_x);
+ _plotterTrackVertexNonBend.region(1).plot(_vtxposnonb_y);
+ _plotterTrackVertexNonBend.region(2).plot(_vtxposnonb_z);
+ _plotterTrackVertexNonBend.region(3).plot(_vtxposnonb_dy);
+ _plotterTrackVertexNonBend.region(4).plot(_vtxposnonb_xAtZ0);
+ _plotterTrackVertexNonBend.region(5).plot(_vtxposnonb_zAtTarget);
+ _plotterTrackVertexNonBend.region(6).plot(_vtxposnonb_angle1);
+ _plotterTrackVertexNonBend.region(7).plot(_vtxposnonb_angle2);
_plotterParticleVertex = aida.analysisFactory().createPlotterFactory().create();
_plotterParticleVertex.createRegions(2,2);
_plotterParticleVertex.region(0).plot(_partvtxpos_x);
@@ -765,16 +911,39 @@
_plotterTrackMult = aida.analysisFactory().createPlotterFactory().create();
_plotterTrackMult.createRegions(1,1);
_plotterTrackMult.region(0).plot(_ntrks_px);
+ _plotterTrackAtConv = aida.analysisFactory().createPlotterFactory().create();
+ _plotterTrackAtConv.createRegions(4,4);
+ _plotterTrackAtConv.region(0).plot(_trk_y_at_conv_top_pos);
+ _plotterTrackAtConv.region(4).plot(_trk_z_at_conv_top_pos);
+ _plotterTrackAtConv.region(8).plot(_trk_y_at_conv_top_pos_fr);
+ _plotterTrackAtConv.region(12).plot(_trk_z_at_conv_top_pos_fr);
+ _plotterTrackAtConv.region(2).plot(_trk_y_at_conv_bot_pos);
+ _plotterTrackAtConv.region(6).plot(_trk_z_at_conv_bot_pos);
+ _plotterTrackAtConv.region(10).plot(_trk_y_at_conv_bot_pos_fr);
+ _plotterTrackAtConv.region(14).plot(_trk_z_at_conv_bot_pos_fr);
+
+ _plotterTrackAtConv.region(1).plot(_trk_y_at_conv_top_neg);
+ _plotterTrackAtConv.region(5).plot(_trk_z_at_conv_top_neg);
+ _plotterTrackAtConv.region(9).plot(_trk_y_at_conv_top_neg_fr);
+ _plotterTrackAtConv.region(13).plot(_trk_z_at_conv_top_neg_fr);
+ _plotterTrackAtConv.region(3).plot(_trk_y_at_conv_bot_neg);
+ _plotterTrackAtConv.region(7).plot(_trk_z_at_conv_bot_neg);
+ _plotterTrackAtConv.region(11).plot(_trk_y_at_conv_bot_neg_fr);
+ _plotterTrackAtConv.region(15).plot(_trk_z_at_conv_bot_neg_fr);
_plotterParticleVertex.setTitle("MC particle Vertex");
_plotterTrackVertex.setTitle("Two Track Vertex");
+ _plotterTrackVertexNonBend.setTitle("Two Track Vertex Non Bend");
_plotterTrackMult.setTitle("Track multiplicity");
+ _plotterTrackAtConv.setTitle("Track @ converter");
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);
+ ((PlotterRegion) _plotterTrackVertexNonBend.region(i)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) _plotterTrackVertexNonBend.region(i)).getPlot().setAllowPopupMenus(true);
if(i==0) {
((PlotterRegion) _plotterTrackMult.region(i)).getPlot().setAllowUserInteraction(true);
((PlotterRegion) _plotterTrackMult.region(i)).getPlot().setAllowPopupMenus(true);
@@ -784,12 +953,89 @@
if(!this.hideFrame) {
this._plotterParticleVertex.show();
this._plotterTrackVertex.show();
- this._plotterTrackMult.show();
+ this._plotterTrackVertexNonBend.show();
+ //this._plotterTrackMult.show();
+ _plotterTrackAtConv.show();
}
}
-
+ private void createWriter() {
+ try {
+ fileWriter = new FileWriter(outputNameTextTuple);
+ } catch (IOException ex) {
+ Logger.getLogger(TwoTrackAnlysis.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ printWriter = new PrintWriter(fileWriter);
+ }
+
+
+
+ /*
+ * Find the path length where the two helix cross in the non-bend plane
+ */
+
+ public double getCrossingS(Track trk1, Track trk2, boolean useFringe) {
+ double slope_1 = ((SeedTrack)trk1).getSeedCandidate().getHelix().slope();
+ double slope_2 = ((SeedTrack)trk2).getSeedCandidate().getHelix().slope();
+ double z0_1 = ((SeedTrack)trk1).getSeedCandidate().getHelix().z0();
+ double z0_2 = ((SeedTrack)trk2).getSeedCandidate().getHelix().z0();
+ double s = getPathLengthCrossingPoint(slope_1, z0_1, slope_2, z0_2);
+ double zAtCross = z0_1 + s*slope_1;
+ //System.out.printf("%s: s1=%f z1=%f s2=%f z2=%f => path=%f and zAtCross=%f\n",this.getClass().getSimpleName(),slope_1,z0_1,slope_2,z0_2,s,zAtCross);
+ return zAtCross;
+ }
+
+
+ public StraightLineTrack[] getSLTs(Track trk1, Track trk2, boolean useFringe) {
+ // find the point on the x- and y-axis by
+ // 1) go outside the fringe region
+ // 2) assume straight lines
+ // 3) solve for the position where the z-position is at the crossing point
+ double zStart = useFringe==true ? -100. : 0.;
+ StraightLineTrack slt1 = this.findSLTAtZ(trk1,zStart,useFringe);
+ StraightLineTrack slt2 = this.findSLTAtZ(trk2,zStart,useFringe);
+ StraightLineTrack[] vv = {slt1,slt2};
+ return vv;
+ }
+
+
+ public StraightLineTrack findSLTAtZ(Track trk1, double zVal, boolean useFringe) {
+ SeedTrack s1 = (SeedTrack) trk1;
+ HelicalTrackFit htf1 = s1.getSeedCandidate().getHelix();
+ HPSTrack hpstrk1 = new HPSTrack(htf1);
+ Hep3Vector pos1;
+ if(useFringe) {
+ pos1 = hpstrk1.getPositionAtZMap(100.0, zVal, 5.0)[0];
+ } else {
+ this.vertexer.extrapolator().setTrack(trk1);
+ pos1 = this.vertexer.extrapolator().extrapolateTrack(zVal);
+ }
+ //System.out.printf("%s: Position1 at edge of fringe %s\n",this.getClass().getSimpleName(),pos1.toString());
+ Helix traj = (Helix)hpstrk1.getTrajectory();
+ if(traj==null) {
+ SpacePoint r0 = new SpacePoint(HelixUtils.PointOnHelix(htf1,0));
+ traj = new Helix(r0,htf1.R(), htf1.phi0(), Math.atan(htf1.slope()));
+ }
+ StraightLineTrack slt1 = this.vertexer.converter().Convert(traj);
+ //System.out.printf("%s: straight line track: x0=%f,y0=%f,z0=%f dz/dx=%f dydx=%f targetY=%f targetZ=%f \n",this.getClass().getSimpleName(),slt1.x0(),slt1.y0(),slt1.z0(),slt1.dzdx(),slt1.dydx(),slt1.TargetYZ()[0],slt1.TargetYZ()[1]);
+ return slt1;
+ }
+
+
+ /*
+ * Calculate the point where two 1-D lines cross
+ * Use SZ coordinate system nomenclature
+ * Parameters:
+ * slope: slope of track in SZ plane
+ * z0: z-coordinate at which S=0
+ */
+ public double getPathLengthCrossingPoint(double slope_1,double z0_1,double slope_2,double z0_2) {
+ double s; //path length
+ s = (z0_1-z0_2)/(slope_2-slope_1);
+ return s;
+ }
+
}