11 #include <tracking/modules/trackingPerformanceEvaluation/TrackingPerformanceEvaluationModule.h>
13 #include <framework/datastore/StoreArray.h>
14 #include <framework/datastore/RelationVector.h>
16 #include <framework/geometry/BFieldManager.h>
18 #include <vxd/geometry/GeoCache.h>
20 #include <mdst/dataobjects/MCParticle.h>
21 #include <mdst/dataobjects/Track.h>
22 #include <mdst/dataobjects/HitPatternCDC.h>
23 #include <mdst/dataobjects/HitPatternVXD.h>
25 #include <pxd/reconstruction/PXDRecoHit.h>
26 #include <svd/reconstruction/SVDRecoHit.h>
27 #include <svd/reconstruction/SVDRecoHit2D.h>
28 #include <cdc/dataobjects/CDCRecoHit.h>
30 #include <tracking/dataobjects/RecoTrack.h>
32 #include <pxd/dataobjects/PXDTrueHit.h>
33 #include <pxd/dataobjects/PXDCluster.h>
34 #include <svd/dataobjects/SVDCluster.h>
35 #include <cdc/dataobjects/CDCHit.h>
37 #include <genfit/KalmanFitterInfo.h>
39 #include <root/TObject.h>
41 #include <boost/foreach.hpp>
50 TrackingPerformanceEvaluationModule::TrackingPerformanceEvaluationModule() :
54 setDescription(
"This module evaluates the tracking package performance");
56 addParam(
"outputFileName",
m_rootFileName,
"Name of output root file.",
57 std::string(
"TrackingPerformanceEvaluation_output.root"));
58 addParam(
"MCParticlesName", m_MCParticlesName,
"Name of MC Particle collection.", std::string(
""));
59 addParam(
"TracksName", m_TracksName,
"Name of Track collection.", std::string(
""));
60 addParam(
"RecoTracksName", m_RecoTracksName,
"Name of RecoTrack collection.", std::string(
"RecoTracks"));
61 addParam(
"MCRecoTracksName", m_MCRecoTracksName,
"Name of MCRecoTrack collection.", std::string(
"MCRecoTracks"));
62 addParam(
"ParticleHypothesis", m_ParticleHypothesis,
"Particle Hypothesis used in the track fit.",
int(211));
66 TrackingPerformanceEvaluationModule::~TrackingPerformanceEvaluationModule()
179 2000, -10, 10,
"x (cm)",
180 2000, -10, 10,
"y (cm)",
184 2000, -30, 40,
"z (cm)",
185 2000, 0, 15,
"r_{t} (cm)",
193 100, 0, 3,
"p_{t} (GeV/c)",
194 1000, 0, 0.2,
"#sigma_{#omega}/#omega",
199 100, 0, 3,
"p_{t} (GeV/c)",
200 100, 0, 0.1,
"#sigma_{z0} (cm)",
215 100, 0, 3,
"p_{t} (GeV/c)",
216 100, 0, 0.1,
"#sigma_{d0} (cm)",
229 50, 0, 2.5,
"p_{t} (GeV/c)",
252 2000, -15, 15,
"x (cm)",
253 2000, -15, 15,
"y (cm)",
257 2000, -30, 40,
"z (cm)",
258 2000, 0, 15,
"r_{t} (cm)",
264 Double_t bins_pt[10 + 1] = {0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 1, 1.5, 2, 3.5};
265 Double_t bins_theta[10 + 1] = {0, 0.25, 0.5, 0.75, 0.75 + 0.32, 0.75 + 2 * 0.32, 0.75 + 3 * 0.32, 0.75 + 4 * 0.32, 0.75 + 5 * 0.32, 2.65, TMath::Pi()};
266 Double_t bins_phi[14 + 1];
267 Double_t width_phi = 2 * TMath::Pi() / 14;
268 for (
int bin = 0; bin < 14 + 1; bin++)
269 bins_phi[bin] = - TMath::Pi() + bin * width_phi;
273 10, bins_pt,
"p_{t} (GeV/c)",
274 10, bins_theta,
"#theta",
275 14, bins_phi,
"#phi" );
278 "entry per Track connected to a MCParticle",
282 "entry per Track with PXD hits connected to a MCParticle",
286 "entry per RecoTrack with PXD hits connected to a MCParticle",
290 "entry per RecoTrack with PXD hits connected to a MCParticle with PXD hits",
294 "entry per MCParticle with PXD hits",
298 "entry per MCRecoTrack connected to the MCParticle",
302 "entry per Track connected to an MCRecoTrack",
309 "entry per Track connected to a positive MCParticle",
313 "entry per MCRecoTrack connected to the positive MCParticle",
317 "entry per Track connected to a positive MCRecoTrack",
326 "entry per Track connected to a negative MCParticle",
330 "entry per MCRecoTrack connected to the negative MCParticle",
334 "entry per Track connected to a negative MCRecoTrack",
350 "entry per MCParticle connected to a Track",
366 bool hasTrack =
false;
367 B2DEBUG(99,
"+++++ 1. loop on MCParticles");
368 BOOST_FOREACH(
MCParticle & mcParticle, mcParticles) {
373 int pdgCode = mcParticle.
getPDG();
374 B2DEBUG(99,
"MCParticle has PDG code " << pdgCode);
376 int nFittedTracksMCRT = 0;
377 int nFittedTracks = 0;
378 int nFittedTrackswPXDHits = 0;
402 if (MCRecoTracks_fromMCParticle.
size() > 0)
403 if (MCRecoTracks_fromMCParticle[0]->hasPXDHits()) {
415 if (MCRecoTracks_fromMCParticle.
size() > 0) {
430 B2DEBUG(99, Tracks_fromMCParticle.
size() <<
" Tracks related to this MCParticle");
432 for (
int trk = 0; trk < (int)Tracks_fromMCParticle.
size(); trk++) {
437 B2WARNING(
" the TrackFitResult is not found!");
449 nFittedTrackswPXDHits++;
460 if (MCRecoTracks_fromMCParticle.
size() > 0) {
481 if (MCRecoTracks_fromMCParticle.
size() > 0)
487 B2DEBUG(99,
"+++++ 2. loop on Tracks");
492 BOOST_FOREACH(
Track & track, tracks) {
494 int nMCParticles = 0;
505 m_h3_Tracks->Fill(momentum.Pt(), momentum.Theta(), momentum.Phi());
511 for (
int layer = 0; layer < 56; layer++) {
515 for (
int layer = 1; layer <= 2; layer++) {
519 for (
int layer = 3; layer <= 6; layer++) {
524 for (
int i = 0; i < N; i++)
531 for (
int mcp = 0; mcp < (int)MCParticles_fromTrack.
size(); mcp++)
543 B2DEBUG(99,
"+++++ 3. loop on MCRecoTracks");
552 BOOST_FOREACH(
RecoTrack & mcRecoTrack, mcRecoTracks) {
555 bool hasRecoTrack =
false;
559 B2DEBUG(99,
"~ " << RecoTracks_fromMCRecoTrack.
size() <<
" RecoTracks related to this MCRecoTrack");
565 B2DEBUG(99,
"~~~ " << MCParticles_fromMCRecoTrack.
size() <<
" MCParticles related to this MCRecoTrack");
566 for (
int mcp = 0; mcp < (int)MCParticles_fromMCRecoTrack.
size(); mcp++) {
570 (MCParticles_fromMCRecoTrack[mcp]);
572 B2DEBUG(99,
"~~~~~ " << RecoTracks_fromMCParticle.
size() <<
" RecoTracks related to this MCParticle");
573 for (
int tc = 0; tc < (int)RecoTracks_fromMCParticle.
size(); tc++)
584 B2DEBUG(99,
"+++++ 4. loop on RecoTracks");
589 BOOST_FOREACH(
RecoTrack & recoTrack, RecoTracks) {
631 double efficiency = num / den ;
632 double efficiencyErr = sqrt(efficiency * (1 - efficiency)) / sqrt(den);
634 double nFittedTracksMCRT = 0;
640 double nRecoTrack = 0;
646 double nMCRecoTrack = 0;
652 double nMCParticles = 0;
659 B2INFO(
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
660 B2INFO(
"~ Tracking Performance Evaluation ~ SHORT SUMMARY ~");
662 B2INFO(
" + overall:");
663 B2INFO(
" efficiency = (" << efficiency * 100 <<
" +/- " << efficiencyErr * 100 <<
")% ");
664 B2INFO(
" purity = " << purity * 100 <<
" +/- " << purityErr * 100 <<
")% ");
666 B2INFO(
" + factorizing geometrical acceptance:");
667 B2INFO(
" efficiency = " << efficiencyMCRT * 100 <<
" +/- " << efficiencyMCRTErr * 100 <<
")% ");
669 B2INFO(
" + pattern recognition:");
670 B2INFO(
" efficiency = " << efficiencyPR * 100 <<
" +/- " << efficiencyPRErr * 100 <<
")% ");
671 B2INFO(
" purity = " << purityPR * 100 <<
" +/- " << purityPRErr * 100 <<
")% ");
673 B2INFO(
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
688 TDirectory* oldDir = gDirectory;
690 TDirectory* dir_multiplicity = oldDir->mkdir(
"multiplicity");
691 dir_multiplicity->cd();
694 while ((obj = nextH_multiplicity()))
697 TDirectory* dir_efficiency = oldDir->mkdir(
"efficiency");
698 dir_efficiency->cd();
700 while ((obj = nextH_efficiency()))
703 TDirectory* dir_trkQuality = oldDir->mkdir(
"trkQuality");
704 dir_trkQuality->cd();
706 while ((obj = nextH_trkQuality()))
709 TDirectory* dir_firstHit = oldDir->mkdir(
"firstHit");
712 while ((obj = nextH_firstHit()))
736 double d0_res = fitResult->
getD0() - mcParticleInfo.
getD0();
737 double phi_res = TMath::ASin(TMath::Sin(fitResult->
getPhi() - mcParticleInfo.
getPhi()));
739 double z0_res = fitResult->
getZ0() - mcParticleInfo.
getZ0();
746 double p_res = (fitResult->
getMomentum().Mag() - mcParticleInfo.
getP()) / mcParticleInfo.
getP();
753 double r_res = fitResult->
getPosition().Perp() - sqrt(mcParticleInfo.
getX() * mcParticleInfo.
getX() + mcParticleInfo.
getY() *
754 mcParticleInfo.
getY());
755 double rtot_res = fitResult->
getPosition().Mag() - sqrt(mcParticleInfo.
getX() * mcParticleInfo.
getX() + mcParticleInfo.
getY() *
756 mcParticleInfo.
getY() + mcParticleInfo.
getZ() * mcParticleInfo.
getZ());
807 double px = momentum.Px();
808 double py = momentum.Py();
809 double pz = momentum.Pz();
810 double pt = momentum.Pt();
811 double p = momentum.Mag();
813 double beta = p / sqrt(p * p + mass * mass);
814 double sinTheta = TMath::Sin(momentum.Theta());
817 d0_err / phi_err * py / pt);
819 z0_err / cotTheta_err * py / pt);
827 m_h2_d0errMSVSpt->Fill(pt, d0_err * beta * p * pow(sinTheta, 3 / 2) / 0.0136);
840 bool hasPXDhit =
false;
841 bool isTrueHit =
false;
843 double d0_err = -999;
844 double z0_err = -999;
848 if ((fitResult != NULL)
855 bool hasCDChit[56] = {
false };
859 for (
int tc = 0; tc < (int)RecoTracks_fromTrack.
size(); tc++) {
861 const std::vector< genfit::TrackPoint* >& tp_vector = RecoTracks_fromTrack[tc]->getHitPointsWithMeasurement();
862 for (
int i = 0; i < (int) tp_vector.size(); i++) {
865 int nMea = tp->getNumRawMeasurements();
866 for (
int mea = 0; mea < nMea; mea++) {
871 std::vector<double> weights;
876 B2WARNING(
" No KalmanFitterInfo associated to the TrackPoint!");
879 TVector3 globalHit(-999, -999, -999);
891 weight = weights.at(mea);
894 double uCoor = pxdHit->
getU();
895 double vCoor = pxdHit->
getV();
904 globalHit = aSensorInfo.
pointToGlobal(TVector3(uCoor, vCoor, 0),
true);
910 if ((
int)pxdth_fromcl.
size() != 0) {
916 for (
int mcp = 0; mcp < (int)MCParticles_fromTrack.
size(); mcp++) {
918 for (
int th = 0; th < (int)trueHit_fromMCParticles.
size(); th++) {
919 if (trueHit_fromMCParticles[th]->getArrayIndex() == trueHitIndex)
926 }
else if (svdHit2D) {
929 weight = weights.at(mea);
932 double uCoor = svdHit2D->
getU();
933 double vCoor = svdHit2D->
getV();
943 globalHit = aSensorInfo.
pointToGlobal(TVector3(uCoor, vCoor, 0),
true);
948 weight = weights.at(mea);
965 globalHit = aSensorInfo.
pointToGlobal(TVector3(uCoor, vCoor, 0),
true);
969 weight = weights.at(mea);
1019 histoList->Add(h_ineff_pt);
1021 TH1F* h_ineff_theta =
createHistogramsRatio(
"hinefftheta",
"inefficiency VS #theta, normalized to MCParticles",
1023 histoList->Add(h_ineff_theta);
1027 histoList->Add(h_ineff_phi);
1030 TH1F* h_ineffMCRT_pt =
createHistogramsRatio(
"hineffMCRTpt",
"inefficiency VS pt, normalized to MCRecoTrack",
1032 histoList->Add(h_ineffMCRT_pt);
1034 TH1F* h_ineffMCRT_theta =
createHistogramsRatio(
"hineffMCRTtheta",
"inefficiency VS #theta, normalized to MCRecoTrack",
1036 histoList->Add(h_ineffMCRT_theta);
1038 TH1F* h_ineffMCRT_phi =
createHistogramsRatio(
"hineffMCRTphi",
"inefficiency VS #phi, normalized to MCRecoTrack",
1040 histoList->Add(h_ineffMCRT_phi);
1048 TH1F* h_MCPwPXDhits_pt =
createHistogramsRatio(
"hMCPwPXDhits",
"fraction of MCParticles with PXD hits VS pt",
1051 histoList->Add(h_MCPwPXDhits_pt);
1054 "fraction of MCParticles with PXD Hits with RecoTracks with PXD hits VS pt",
1057 histoList->Add(h_RTwPXDhitsMCPwPXDHits_pt);
1059 TH1F* h_wPXDhits_pt =
createHistogramsRatio(
"hTrkswPXDhits",
"fraction of tracks with PXD hits VS pt",
1062 histoList->Add(h_wPXDhits_pt);
1067 histoList->Add(h_eff_pt);
1072 histoList->Add(h_eff_theta);
1076 histoList->Add(h_eff_phi);
1081 histoList->Add(h_effMCRT_pt);
1083 TH1F* h_effMCRT_theta =
createHistogramsRatio(
"heffMCRTtheta",
"efficiency VS #theta, normalized to MCRecoTrack",
1085 histoList->Add(h_effMCRT_theta);
1087 TH1F* h_effMCRT_phi =
createHistogramsRatio(
"heffMCRTphi",
"efficiency VS #phi, normalized to MCRecoTrack",
1089 histoList->Add(h_effMCRT_phi);
1094 TH1F* h_eff_pt_plus =
createHistogramsRatio(
"heffpt_plus",
"efficiency VS pt, normalized to positive MCParticles",
1096 histoList->Add(h_eff_pt_plus);
1099 TH1F* h_eff_theta_plus =
createHistogramsRatio(
"hefftheta_plus",
"efficiency VS #theta, normalized to positive MCParticles",
1101 histoList->Add(h_eff_theta_plus);
1103 TH1F* h_eff_phi_plus =
createHistogramsRatio(
"heffphi_plus",
"efficiency VS #phi, normalized to positive MCParticles",
1105 histoList->Add(h_eff_phi_plus);
1108 TH1F* h_effMCRT_pt_plus =
createHistogramsRatio(
"heffMCRTpt_plus",
"efficiency VS pt, normalized to positive MCRecoTrack",
1110 histoList->Add(h_effMCRT_pt_plus);
1112 TH1F* h_effMCRT_theta_plus =
createHistogramsRatio(
"heffMCRTtheta_plus",
"efficiency VS #theta, normalized to positive MCRecoTrack",
1114 histoList->Add(h_effMCRT_theta_plus);
1116 TH1F* h_effMCRT_phi_plus =
createHistogramsRatio(
"heffMCRTphi_plus",
"efficiency VS #phi, normalized to positive MCRecoTrack",
1118 histoList->Add(h_effMCRT_phi_plus);
1123 TH1F* h_eff_pt_minus =
createHistogramsRatio(
"heffpt_minus",
"efficiency VS pt, normalized to positive MCParticles",
1125 histoList->Add(h_eff_pt_minus);
1128 TH1F* h_eff_theta_minus =
createHistogramsRatio(
"hefftheta_minus",
"efficiency VS #theta, normalized to positive MCParticles",
1130 histoList->Add(h_eff_theta_minus);
1132 TH1F* h_eff_phi_minus =
createHistogramsRatio(
"heffphi_minus",
"efficiency VS #phi, normalized to positive MCParticles",
1134 histoList->Add(h_eff_phi_minus);
1137 TH1F* h_effMCRT_pt_minus =
createHistogramsRatio(
"heffMCRTpt_minus",
"efficiency VS pt, normalized to positive MCRecoTrack",
1139 histoList->Add(h_effMCRT_pt_minus);
1143 histoList->Add(h_effMCRT_theta_minus);
1145 TH1F* h_effMCRT_phi_minus =
createHistogramsRatio(
"heffMCRTphi_minus",
"efficiency VS #phi, normalized to positive MCRecoTrack",
1147 histoList->Add(h_effMCRT_phi_minus);
1150 TH1F* h_effPR =
createHistogramsRatio(
"heffPR",
"PR efficiency VS VXD Layer, normalized to MCRecoTrack",
1152 histoList->Add(h_effPR);
1157 histoList->Add(h_effVXDHitFit);
1161 histoList->Add(h_effCDCHitFit);
1175 return (isPrimary && isChargedStable);