11 #include <tracking/modules/trackingPerformanceEvaluation/V0findingPerformanceEvaluationModule.h>
13 #include <tracking/dataobjects/MCParticleInfo.h>
14 #include <tracking/dataobjects/V0ValidationVertex.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/datastore/RelationVector.h>
19 #include <mdst/dataobjects/Track.h>
21 #include <framework/geometry/BFieldManager.h>
23 #include <root/TAxis.h>
24 #include <root/TObject.h>
26 #include <boost/foreach.hpp>
36 V0findingPerformanceEvaluationModule::V0findingPerformanceEvaluationModule() :
40 setDescription(
"This module evaluates the V0 finding package performance");
43 std::string(
"V0findingPerformanceEvaluation_output.root"));
44 addParam(
"V0sName", m_V0sName,
"Name of V0 collection.", std::string(
"V0ValidationVertexs"));
49 V0findingPerformanceEvaluationModule::~V0findingPerformanceEvaluationModule()
60 v0ValidationVertices.isRequired(
m_V0sName);
88 m_h2_vtxTvsR_err =
createHistogram2D(
"h2vtxTerrVsR",
"vtxT error vs R", 100, 0, 100,
"R (cm)", 100, 0, 0.3,
"#sigma_{vtxT} (cm)",
113 m_h2_mom =
createHistogram2D(
"h2mom",
"reco VS true momentum", 100, 0, 3,
"V0 mom (GeV/c)", 100, 0, 3,
"MC mom (GeV/c)",
115 m_h2_mass =
createHistogram2D(
"h2mass",
"reco VS true mass", 100, 0, 1.5,
"V0 mass (GeV/c2)", 100, 0, 1.5,
"MC mass (GeV/c)",
119 Double_t bins_pt[9 + 1] = {0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5, 1, 2, 3.5};
120 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()};
121 Double_t bins_phi[14 + 1];
122 Double_t width_phi = 2 * TMath::Pi() / 14;
123 for (
int bin = 0; bin < 14 + 1; bin++)
124 bins_phi[bin] = - TMath::Pi() + bin * width_phi;
133 9, bins_pt,
"p_{t} (GeV/c)",
134 10, bins_theta,
"#theta",
135 14, bins_phi,
"#phi" );
138 "entry per V0 connected to a MCParticle",
146 "entry per MCParticle connected to a V0",
162 B2DEBUG(99,
"+++++ 1. loop on MCParticles");
163 BOOST_FOREACH(
MCParticle & mcParticle, mcParticles) {
165 if (!
isV0(mcParticle))
172 if (nMatchedDau != 2)
175 int pdgCode = mcParticle.
getPDG();
176 B2DEBUG(99,
"MCParticle has PDG code " << pdgCode);
181 TVector3 MC_prodvtx = mcParticle.
getVertex();
184 float MC_mass = mcParticle.
getMass();
185 TVector3 MC_FL = MC_vtx - MC_prodvtx;
186 float flightR = sqrt(MC_FL.X() * MC_FL.X() + MC_FL.Y() * MC_FL.Y());
195 DataStore::getRelationsWithObj<V0ValidationVertex>(&mcParticle,
m_V0sName);
199 if (V0s_toMCParticle.
size() > 0)
202 for (
int v0 = 0; v0 < (int)V0s_toMCParticle.
size(); v0++) {
204 TVector3 V0_vtx = V0s_toMCParticle[v0]->getVertexPosition();
205 float V0_mom = V0s_toMCParticle[v0]->getFittedMomentum();
206 float V0_chi2 = V0s_toMCParticle[v0]->getVertexChi2();
207 float V0_mass = V0s_toMCParticle[v0]->getFittedInvariantMass();
208 TMatrixDSym V0_cov = V0s_toMCParticle[v0]->getVertexPositionCovariance();
226 m_h1_vtxX_pll->Fill((V0_vtx.X() - MC_vtx.X()) / sqrt(V0_cov[0][0]));
227 m_h1_vtxY_pll->Fill((V0_vtx.Y() - MC_vtx.Y()) / sqrt(V0_cov[1][1]));
228 m_h1_vtxZ_pll->Fill((V0_vtx.Z() - MC_vtx.Z()) / sqrt(V0_cov[2][2]));
238 B2DEBUG(99,
"+++++ 2. loop on V0s");
246 int nMCParticles = 0;
253 nMCParticles = MCParticles_fromV0.
size();
255 if (nMCParticles == 0)
275 double efficiency = num / den ;
276 double efficiencyErr = sqrt(efficiency * (1 - efficiency)) / sqrt(den);
278 double nMCParticles = 0;
285 B2INFO(
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
286 B2INFO(
"~ V0 Finding Performance Evaluation ~ SHORT SUMMARY ~");
288 B2INFO(
" + overall:");
289 B2INFO(
" efficiency = (" << efficiency * 100 <<
" +/- " << efficiencyErr * 100 <<
")% ");
290 B2INFO(
" purity = (" << purity * 100 <<
" +/- " << purityErr * 100 <<
")% ");
292 B2INFO(
"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
300 for (
int bin = 0; bin < h_eff_R->GetXaxis()->GetNbins(); bin++) {
307 eff = (double)num / den;
308 err = sqrt(eff * (1 - eff)) / sqrt(den);
311 h_eff_R->SetBinContent(bin + 1, eff);
312 h_eff_R->SetBinError(bin + 1, err);
324 TDirectory* oldDir = gDirectory;
326 TDirectory* dir_multiplicity = oldDir->mkdir(
"multiplicity");
327 dir_multiplicity->cd();
330 while ((obj = nextH_multiplicity()))
333 TDirectory* dir_efficiency = oldDir->mkdir(
"efficiency");
334 dir_efficiency->cd();
336 while ((obj = nextH_efficiency()))
339 TDirectory* dir_purity = oldDir->mkdir(
"purity");
342 while ((obj = nextH_purity()))
345 TDirectory* dir_trkQuality = oldDir->mkdir(
"trkQuality");
346 dir_trkQuality->cd();
348 while ((obj = nextH_trkQuality()))
352 while ((obj = nextH()))
365 bool isGamma =
false;
366 if (abs(the_mcParticle.
getPDG()) == 22)
370 if (abs(the_mcParticle.
getPDG()) == 310)
374 if (abs(the_mcParticle.
getPDG()) == 311)
377 bool isLambda =
false;
378 if (abs(the_mcParticle.
getPDG()) == 3122)
381 bool twoProngs =
false;
382 bool twoChargedProngs =
false;
389 twoChargedProngs =
true;
391 return ((isGamma || isK_S0 || isK_0 || isLambda) && twoChargedProngs);
400 std::vector< MCParticle* > MCPart_dau = the_mcParticle.
getDaughters();
405 RelationVector<Track> Tracks_fromMCParticle_0 = DataStore::getRelationsWithObj<Track>(MCPart_dau[0]);
406 if (Tracks_fromMCParticle_0.
size() > 0)
409 RelationVector<Track> Tracks_fromMCParticle_1 = DataStore::getRelationsWithObj<Track>(MCPart_dau[1]);
410 if (Tracks_fromMCParticle_1.
size() > 0)