hps-java/src/main/java/org/lcsim/hps/users/phansson
diff -u -r1.12 -r1.13
--- TwoTrackAnlysis.java 15 Mar 2013 02:52:35 -0000 1.12
+++ TwoTrackAnlysis.java 18 Mar 2013 22:33:44 -0000 1.13
@@ -15,21 +15,22 @@
import java.util.Map;
import java.util.logging.Level;
import java.util.logging.Logger;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.tracker.silicon.SiSensor;
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.HelicalTrackCross;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+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;
import org.lcsim.hps.recon.tracking.HPSTrack;
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;
@@ -62,45 +63,19 @@
private String stereoHitCollectionName = "RotatedHelicalTrackHits";
private String triggerDecisionCollectionName = "TriggerBank";
private String MCParticleCollectionName = "MCParticle";
- private String _stripClusterCollectionName = "HelicalTrackStrips";
+ private String _stripClusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
private boolean _debug;
+ private HitIdentifier _ID = new HitIdentifier();
private TwoTrackVertexer vertexer = new TwoTrackVertexer();
private TwoParticleVertexer particleVertexer = new TwoParticleVertexer();
private IPlotter _plotterParticleVertex;
private IPlotter _plotterTrackVertex;
- private IPlotter _plotterNhits;
- private IPlotter _plotterRes;
- private IPlotter _plotterResAll;
- 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;
- private IHistogram1D _nhits;
- private IHistogram1D _nhits_vtxpospos;
- private IHistogram1D _nhits_vtxposneg;
- private IHistogram1D _layershit;
- private IHistogram1D _layershit_vtxpospos;
- private IHistogram1D _layershit_vtxposneg;
- private IHistogram1D _trk1_res1_z;
- private IHistogram1D _trk2_res1_z;
- private IHistogram1D _trk1_res1_y;
- private IHistogram1D _trk2_res1_y;
- private IHistogram1D _trk1_res2_z;
- private IHistogram1D _trk2_res2_z;
- private IHistogram1D _trk1_res2_y;
- private IHistogram1D _trk2_res2_y;
- private IHistogram1D _trk1_res1_z_all;
- private IHistogram1D _trk2_res1_z_all;
- private IHistogram1D _trk1_res1_y_all;
- private IHistogram1D _trk2_res1_y_all;
- private IHistogram1D _trk1_res2_z_all;
- private IHistogram1D _trk2_res2_z_all;
- private IHistogram1D _trk1_res2_y_all;
- private IHistogram1D _trk2_res2_y_all;
@@ -181,7 +156,7 @@
}
}
- if(tracklist.size()!=2) {
+ if(tracklist.size()==0) {
if(_debug) {
System.out.printf("%s: event %d has only %d tracks \n",this.getClass().getSimpleName(),event.getEventNumber(),tracklist.size());
}
@@ -199,88 +174,54 @@
}
}
+ List<Track> tracks = new ArrayList<Track>();
+ Hep3Vector vtxPos = null;
+ Hep3Vector vtxPosFringe = null;
- //DEBUG
- //if(_debug) {
- {
-
- if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": fill some debug hists for event " + event.getEventNumber());
- }
- for(int i=0;i<tracklist.size();++i) {
- Track trk = tracklist.get(i);
- SeedTrack st1 = (SeedTrack) trk;
- HelicalTrackFit helix = st1.getSeedCandidate().getHelix();
- HashMap<Integer,HelicalTrackHit> hits = this.getHitMap(trk.getTrackerHits(),helix);
- if(i==0) {
- if(hits.get(1) != null) {
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hits.get(1), helix, true);
- _trk1_res1_z_all.fill(res.get("resz"));
- _trk1_res1_y_all.fill(res.get("resy"));
- }
- if(hits.get(3) != null) {
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hits.get(3), helix, true);
- _trk1_res2_z_all.fill(res.get("resz"));
- _trk1_res2_y_all.fill(res.get("resy"));
- }
- }
- if(i==1) {
- if(hits.get(1) != null) {
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hits.get(1), helix, true);
- _trk2_res1_z_all.fill(res.get("resz"));
- _trk2_res1_y_all.fill(res.get("resy"));
- }
- if(hits.get(3) != null) {
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hits.get(3), helix, true);
- _trk2_res2_z_all.fill(res.get("resz"));
- _trk2_res2_y_all.fill(res.get("resy"));
- }
- }
- }
- }
- //}
+ if(tracklist.size()>1) {
+ Track trk1 = tracklist.get(0);
+ Track trk2 = tracklist.get(1);
+ this.vertexer.setTracks(tracklist.get(0), trk2);
+ 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());
+ vtxPos = null;
+ }
-
- 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;
+ vtxPosFringe = this.vertexer.getVertexWithFringe();
+ if(this._debug)
+ System.out.printf("%s: vtxPosFringe=%s\n", this.getClass().getSimpleName(),vtxPosFringe.toString());
+ if(vtxPosFringe.x() != vtxPosFringe.x()) {
+ System.out.printf("%s: vtxPosFringe is NaN -> Skip\n",this.getClass().getSimpleName());
+ vtxPos = null;
+ }
}
- Hep3Vector vtxPosFringe = this.vertexer.getVertexWithFringe();
- if(this._debug)
- System.out.printf("%s: vtxPosFringe=%s\n", this.getClass().getSimpleName(),vtxPosFringe.toString());
- if(vtxPosFringe.x() != vtxPosFringe.x()) {
- System.out.printf("%s: vtxPosFringe is NaN -> Skip\n",this.getClass().getSimpleName());
- return;
+ if(vtxPos!=null) {
+ this._vtxpos_x.fill(vtxPos.x());
+ this._vtxpos_y.fill(vtxPos.y());
+ this._vtxpos_z.fill(vtxPos.z());
}
+ List<HPSEcalCluster> clusters = new ArrayList<HPSEcalCluster>();
if(!event.hasCollection(HPSEcalCluster.class, ecalClusterCollectionName)) {
if(_debug) {
System.out.println(this.getClass().getSimpleName() + ": event doesn't have a ecal cluster collection ");
}
- }
-
+ } else {
+ clusters = event.get(HPSEcalCluster.class, ecalClusterCollectionName);
- List<HPSEcalCluster> clusters = event.get(HPSEcalCluster.class, ecalClusterCollectionName);
-
- if(_debug) {
- System.out.println(this.getClass().getSimpleName() + ": found " + clusters.size() + " ecal clusters " + event.getEventNumber());
+ if(_debug) {
+ System.out.println(this.getClass().getSimpleName() + ": found " + clusters.size() + " ecal clusters " + event.getEventNumber());
+ }
}
-
Hep3Vector vtxPosMC = null;
MCParticle electron=null;
MCParticle positron=null;
@@ -314,7 +255,6 @@
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());
@@ -322,37 +262,11 @@
}
}
-
-
- 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);
- }
-
-
- // Calculate residuals
-
totalTwoTrackEvents++;
- this.fillTextTuple(electron, positron, trk1, trk2, vtxPosMC, vtxPos, vtxPosFringe, clusters, event);
+ this.fillTextTuple(electron, positron, tracklist, vtxPosMC, vtxPos, vtxPosFringe, clusters, event);
if(this._debug) System.out.println(this.getClass().getSimpleName() + ": # two track events so far = "+totalTwoTrackEvents);
@@ -397,7 +311,7 @@
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, Hep3Vector vtxPosFr, List<HPSEcalCluster> clusters, EventHeader event) {
+ private void fillTextTuple(MCParticle e, MCParticle p, List<Track> tracks, Hep3Vector vtxPosParticle, Hep3Vector vtxPos, Hep3Vector vtxPosFr, List<HPSEcalCluster> clusters, EventHeader event) {
if(doPrintBranchInfoLine) {
String br_line = "";
br_line+="evtnr/I:";
@@ -413,8 +327,9 @@
//br_line+="trk1_ures"+iLayer+"/F:"+"trk1_ureserr"+iLayer+"/F:";
}
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk1_strip"+iLayer+"_u/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:";
+ 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:";
@@ -426,15 +341,16 @@
//br_line+="trk2_ures"+iLayer+"/F:"+"trk2_ureserr"+iLayer+"/F:";
}
for(int iLayer=1;iLayer<=10;++iLayer) br_line+="trk2_strip"+iLayer+"_u/F:";
- for(int iLayer=1;iLayer<=10;++iLayer) br_line+="top_stereo"+iLayer+"_n/F:";
- for(int iLayer=1;iLayer<=10;++iLayer) br_line+="bot_stereo"+iLayer+"_n/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:";
+
+ 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+="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+="trk1_fr_conv_x/F:trk1_fr_conv_y/F:trk1_fr_conv_z/F:";
- br_line+="trk2_fr_conv_x/F:trk2_fr_conv_y/F:trk2_fr_conv_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:";
@@ -449,70 +365,95 @@
//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. );
- //Track properties
- SeedTrack st1 = (SeedTrack) trk1;
- HelicalTrackFit helix1 = st1.getSeedCandidate().getHelix();
- printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f ",helix1.dca(),helix1.phi0(),helix1.R(),helix1.z0(),helix1.slope());
- 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,helix1);
- //printWriter.format("\n%s\n","X1");
- // 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);
- }
- //Get the helix for residual calculation
- for(int iLayer=0;iLayer<5;++iLayer) {
- HelicalTrackHit hitOnLayer = hits1.get(iLayer*2+1);// = this.getHitOnLayer(iLayer, hitsOnTrack);
- if (hitOnLayer != null) {
- //printWriter.format("\n%s\n","X11");
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hitOnLayer, helix1, true);
- printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("resy"),res.get("resz"));
- printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("erry"),res.get("errz"));
- printWriter.format("%5.5f %5.5f ", res.get("drphi"),res.get("msdrphi"));
- printWriter.format("%5.5f %5.5f ", res.get("dz_res"),res.get("msdz"));
-
- //DEBUG histos
- if(iLayer==0) {
- _trk1_res1_z.fill(res.get("resz"));
- _trk1_res1_y.fill(res.get("resy"));
- }
- if(iLayer==1) {
- _trk1_res2_z.fill(res.get("resz"));
- _trk1_res2_y.fill(res.get("resy"));
+ for (int itrk=0;itrk<2;itrk++) {
+ Track trk1 = null;
+ if(tracks.size()>itrk) trk1 = tracks.get(itrk);
+
+ if(trk1!=null) {
+ SeedTrack st1 = (SeedTrack) trk1;
+ HelicalTrackFit helix1 = st1.getSeedCandidate().getHelix();
+ List<TrackerHit> hitsOnTrack1 = trk1.getTrackerHits();
+ HashMap<Integer,HelicalTrackHit> hits1 = getHitMap(hitsOnTrack1,helix1);
+
+ printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f ",helix1.dca(),helix1.phi0(),helix1.R(),helix1.z0(),helix1.slope());
+ 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());
+ // 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);
+ }
+
+ //Get the helix for residual calculation
+ for(int iLayer=0;iLayer<5;++iLayer) {
+ HelicalTrackHit hitOnLayer = hits1.get(iLayer*2+1);// = this.getHitOnLayer(iLayer, hitsOnTrack);
+ if (hitOnLayer != null) {
+ //printWriter.format("\n%s\n","X11");
+ Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hitOnLayer, helix1, true);
+ printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("resy"),res.get("resz"));
+ printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("erry"),res.get("errz"));
+ printWriter.format("%5.5f %5.5f ", res.get("drphi"),res.get("msdrphi"));
+ printWriter.format("%5.5f %5.5f ", res.get("dz_res"),res.get("msdz"));
+
+ }
+ else {
+ printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
+ 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;
+ 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);
+ }
+ //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
+ 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
+
+
}
else {
- printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
+
+ printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f ",-9999999.9,-9999999.9,-9999999.9,-9999999.9,-9999999.9);
+ printWriter.format("%5d %5.5f %5.5f %5.5f %5.5f %5d ",-9999999,-9999999.9, -9999999.9,-9999999.9,-9999999.9,-9999999);
+ for(int iLayer=0;iLayer<5;++iLayer) {
+ printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ }
+ for(int iLayer=0;iLayer<5;++iLayer) {
+ printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
+ printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
+ printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
+ }
+ for(int iLayer=1;iLayer<=10;++iLayer) {
+ printWriter.format("%5.5f ", -99999999.9);
+ }
+ 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 %5.5f ",-9999999.9,-9999999.9,-9999999.9); //note rotation from JLab->tracking
+
}
}
-
+ //printWriter.format("\n%s\n","X1");
+
+
+
//printWriter.format("\n%s\n","X2");
+
- List<HelicalTrackHit> stereoHits = new ArrayList<HelicalTrackHit>();
- if(event.hasCollection(HelicalTrackHit.class, stereoHitCollectionName)) {
- stereoHits = event.get(HelicalTrackHit.class, stereoHitCollectionName);
- } else {
- System.out.printf("%s: no collection with stereo hits? \"%s\"\n",this.getClass().getSimpleName(),this.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>> allstriphits = this.getAllStripHitsMap(event,true);
+
+ 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);
@@ -522,66 +463,23 @@
if(allstriphits.containsKey(iLayer)) printWriter.format("%5d ", allstriphits.get(iLayer).size());
else printWriter.format("%5d ", -99999999);
}
- SeedTrack st2 = (SeedTrack) trk2;
- HelicalTrackFit helix2 = st2.getSeedCandidate().getHelix();
- printWriter.format("%5.5f %5.5f %5.5f %5.5f %5.5f ",helix2.dca(),helix2.phi0(),helix2.R(),helix2.z0(),helix2.slope());
- 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,helix2);
- 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);
- }
-
- //printWriter.format("\n%s\n","X3");
-
- for(int iLayer=0;iLayer<5;++iLayer) {
- HelicalTrackHit hitOnLayer = hits2.get(iLayer*2+1);// = this.getHitOnLayer(iLayer, hitsOnTrack);
- if (hitOnLayer != null) {
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hitOnLayer, helix2, true);
- printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("resy"),res.get("resz"));
- printWriter.format("%5.5f %5.5f %5.5f ", 0.0,res.get("erry"),res.get("errz"));
- printWriter.format("%5.5f %5.5f ", res.get("drphi"),res.get("msdrphi"));
- printWriter.format("%5.5f %5.5f ", res.get("dz_res"),res.get("msdz"));
-
-
- if(res.get("erry")>19.7) {
- System.out.printf("%s: y residual error is too large = %f for run %d and event %d \n",this.getClass().getSimpleName(),res.get("erry"),event.getRunNumber(),event.getEventNumber());
- //System.exit(1);
- }
- //DEBUG histos
- if(iLayer==0) {
- _trk2_res1_z.fill(res.get("resz"));
- _trk2_res1_y.fill(res.get("resy"));
- }
- if(iLayer==1) {
- _trk2_res2_z.fill(res.get("resz"));
- _trk2_res2_y.fill(res.get("resy"));
- }
- }
- else {
- printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f %5.5f ", -9999999.9, -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f ", -9999999.9, -9999999.9);
- printWriter.format("%5.5f %5.5f ", -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);
+ List<HelicalTrackHit> stereoHits = new ArrayList<HelicalTrackHit>();
+ if(event.hasCollection(HelicalTrackHit.class, stereoHitCollectionName)) {
+ stereoHits = event.get(HelicalTrackHit.class, stereoHitCollectionName);
+ } else {
+ System.out.printf("%s: no collection with stereo hits? \"%s\"\n",this.getClass().getSimpleName(),this.stereoHitCollectionName);
+
}
+
+
HashMap<Integer,List<HelicalTrackHit>> allstereohits = getAllStereoHitsMap(stereoHits,true);
- for(int iLayer=1;iLayer<=10;++iLayer) {
+ for(int iLayer=1;iLayer<=10;iLayer+=2) {
if(allstereohits.containsKey(iLayer)) printWriter.format("%5d ", allstereohits.get(iLayer).size());
else printWriter.format("%5d ", -99999999);
}
allstereohits = getAllStereoHitsMap(stereoHits,false);
- for(int iLayer=1;iLayer<=10;++iLayer) {
+ for(int iLayer=1;iLayer<=10;iLayer+=2) {
if(allstereohits.containsKey(iLayer)) printWriter.format("%5d ", allstereohits.get(iLayer).size());
else printWriter.format("%5d ", -99999999);
}
@@ -591,21 +489,10 @@
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() );
- printWriter.format("%5.5f %5.5f %5.5f ", vtxPosFr.x(),vtxPosFr.y(),vtxPosFr.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
- 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
- HPSTrack hpstrk2 = new HPSTrack(helix2);
- Hep3Vector posAtConverterFringe2 = hpstrk2.getPositionAtZMap(100., SvtTrackExtrapolator.HARP_POSITION, 5.0)[0];
- printWriter.format("%5.5f %5.5f %5.5f ", posAtConverterFringe2.z(),posAtConverterFringe2.x(),posAtConverterFringe2.y()); //note rotation from JLab->tracking
+ 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() );
+ 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) {
@@ -691,21 +578,24 @@
return map;
}
- private HashMap<Integer,List<HelicalTrackStrip>> getAllStripHitsMap(EventHeader event, boolean top) {
- HashMap<Integer,List<HelicalTrackStrip>> map = new HashMap<Integer,List<HelicalTrackStrip>>();
- if(!event.hasCollection(HelicalTrackStrip.class, this._stripClusterCollectionName)) {
+ private HashMap<Integer,List<SiTrackerHitStrip1D>> getAllStripHitsMap(EventHeader event, boolean top) {
+ HashMap<Integer,List<SiTrackerHitStrip1D>> map = new HashMap<Integer,List<SiTrackerHitStrip1D>>();
+ if(!event.hasCollection(SiTrackerHitStrip1D.class, this._stripClusterCollectionName)) {
return map;
}
- List<HelicalTrackStrip> strips = event.get(HelicalTrackStrip.class, this._stripClusterCollectionName);
- if(this._debug) System.out.printf("%s: asking strips in the %s\n", this.getClass().getSimpleName(),(top?"top":"bottom"));
- for(HelicalTrackStrip strip : strips) {
- if(top && strip.origin().y()<0 ) continue;
- else if(!top && strip.origin().y()>0) continue;
- if(this._debug) System.out.printf("%s: strip at origin %s is selected\n", this.getClass().getSimpleName(),strip.origin().toString());
- if(!map.containsKey(strip.layer())) {
- map.put(strip.layer(), new ArrayList<HelicalTrackStrip>());
+ List<SiTrackerHitStrip1D> strips = event.get(SiTrackerHitStrip1D.class, this._stripClusterCollectionName);
+ if(this._debug) System.out.printf("%s: found %d strips in clollection, asking strips in the %s\n", this.getClass().getSimpleName(),strips.size(),(top?"top":"bottom"));
+ for(SiTrackerHitStrip1D strip : strips) {
+ IDetectorElement de = strip.getSensor();
+ SiSensor sensor = (SiSensor) de;
+ int lyr = _ID.getLayer(de);
+ if(!top && SvtUtils.getInstance().isTopLayer(sensor)) continue;
+ else if (top && !SvtUtils.getInstance().isTopLayer(sensor)) continue;
+ if(this._debug) System.out.printf("%s: strip \"%s\" at %s is selected\n", this.getClass().getSimpleName(),_ID.getName(de),strip.getPositionAsVector().toString());
+ if(!map.containsKey(lyr)) {
+ map.put(lyr, new ArrayList<SiTrackerHitStrip1D>());
}
- map.get(strip.layer()).add(strip);
+ map.get(lyr).add(strip);
}
return map;
@@ -735,129 +625,39 @@
_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);
- _trk1_res1_z = aida.histogram1D("Track1 layer1 residual z", 50,-20,20);
- _trk2_res1_z = aida.histogram1D("Track2 layer1 residual z", 50,-20,20);
- _trk1_res1_y = aida.histogram1D("Track1 layer1 residual y", 50,-20,20);
- _trk2_res1_y = aida.histogram1D("Track2 layer1 residual y", 50,-20,20);
- _trk1_res2_z = aida.histogram1D("Track1 layer2 residual z", 50,-20,20);
- _trk2_res2_z = aida.histogram1D("Track2 layer2 residual z", 50,-20,20);
- _trk1_res2_y = aida.histogram1D("Track1 layer2 residual y", 50,-20,20);
- _trk2_res2_y = aida.histogram1D("Track2 layer2 residual y", 50,-20,20);
-
- _trk1_res1_z_all = aida.histogram1D("Track1 layer1 residual z all", 50,-20,20);
- _trk2_res1_z_all = aida.histogram1D("Track2 layer1 residual z all", 50,-20,20);
- _trk1_res1_y_all = aida.histogram1D("Track1 layer1 residual y all", 50,-20,20);
- _trk2_res1_y_all = aida.histogram1D("Track2 layer1 residual y all", 50,-20,20);
- _trk1_res2_z_all = aida.histogram1D("Track1 layer2 residual z all", 50,-20,20);
- _trk2_res2_z_all = aida.histogram1D("Track2 layer2 residual z all", 50,-20,20);
- _trk1_res2_y_all = aida.histogram1D("Track1 layer2 residual y all", 50,-20,20);
- _trk2_res2_y_all = aida.histogram1D("Track2 layer2 residual y all", 50,-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);
- _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);
- _plotterRes = aida.analysisFactory().createPlotterFactory().create();
- _plotterRes.createRegions(2,2);
- _plotterRes.region(0).plot(this._trk1_res1_z);
- _plotterRes.region(0).plot(this._trk2_res1_z,"mode=overlay");
- _plotterRes.region(1).plot(this._trk1_res2_z);
- _plotterRes.region(1).plot(this._trk2_res2_z,"mode=overlay");
- _plotterRes.region(2).plot(this._trk1_res1_y);
- _plotterRes.region(2).plot(this._trk2_res1_y,"mode=overlay");
- _plotterRes.region(3).plot(this._trk1_res2_y);
- _plotterRes.region(3).plot(this._trk2_res2_y,"mode=overlay");
- _plotterResAll = aida.analysisFactory().createPlotterFactory().create();
- _plotterResAll.createRegions(2,2);
- _plotterResAll.region(0).plot(this._trk1_res1_z_all);
- _plotterResAll.region(0).plot(this._trk2_res1_z_all,"mode=overlay");
- _plotterResAll.region(1).plot(this._trk1_res2_z_all);
- _plotterResAll.region(1).plot(this._trk2_res2_z_all,"mode=overlay");
- _plotterResAll.region(2).plot(this._trk1_res1_y_all);
- _plotterResAll.region(2).plot(this._trk2_res1_y_all,"mode=overlay");
- _plotterResAll.region(3).plot(this._trk1_res2_y_all);
- _plotterResAll.region(3).plot(this._trk2_res2_y_all,"mode=overlay");
_plotterParticleVertex.setTitle("MC particle Vertex");
_plotterTrackVertex.setTitle("Two Track Vertex");
- _plotterNhits.setTitle("Hits on track");
- _plotterRes.setTitle("Track Residuals");
- _plotterRes.setTitle("Track Residuals All");
- for(int i=0;i<4;++i) {
+ 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) _plotterRes.region(i)).getPlot().setAllowUserInteraction(true);
- ((PlotterRegion) _plotterRes.region(i)).getPlot().setAllowPopupMenus(true);
- ((PlotterRegion) _plotterResAll.region(i)).getPlot().setAllowUserInteraction(true);
- ((PlotterRegion) _plotterResAll.region(i)).getPlot().setAllowPopupMenus(true);
}
if(!this.hideFrame) {
this._plotterParticleVertex.show();
this._plotterTrackVertex.show();
- this._plotterNhits.show();
- this._plotterRes.show();
- this._plotterResAll.show();
}
}
-
- void fillResidualPlots(Track track) {
- SeedTrack st2 = (SeedTrack) track;
- HelicalTrackFit helix2 = st2.getSeedCandidate().getHelix();
- for(TrackerHit hit: track.getTrackerHits()) {
- HelicalTrackHit hth = (HelicalTrackHit) hit;
- Map<String,Double> res = TrackUtils.calculateTrackHitResidual(hth, helix2, true);
- //DEBUG histos
- if(hth.Layer()==1) {
- _trk1_res1_z.fill(res.get("resz"));
- _trk1_res1_y.fill(res.get("resy"));
- }
- if(hth.Layer()==3) {
- _trk1_res2_z.fill(res.get("resz"));
- _trk1_res2_y.fill(res.get("resy"));
- }
- }
- }
-
}