Commit in hps-java/src/main/java/org/lcsim/hps/users/phansson on MAIN | |||
TwoTrackAnlysis.java | +160 | -360 | 1.12 -> 1.13 |
Added 1 track flexibility
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")); - } - } - } -
}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1