22 #include <TDatabasePDG.h>
26 #include <TParticlePDG.h>
41 , _svtxEvalStack(nullptr)
43 , _truthContainer(nullptr)
53 std::cout <<
"QAG4SimulationUpsilon::InitRun - Fatal Error - "
54 <<
"unable to find DST node "
55 <<
"G4TruthInfo" << std::endl;
78 ";Reco p_{T}/Truth p_{T}", 500, 0, 2);
82 ";Truth p_{T} [GeV/c];Reco p_{T}/Truth p_{T}", 200, 0, 50, 500, 0, 2);
88 "Reco tracks at truth p_{T};Truth p_{T} [GeV/c]", 200, 0.1, 50.5);
93 ";Truth p_{T} [GeV/c];Track count / bin", 200, 0.1, 50.5);
99 "Reco tracks at truth #eta;Truth #eta", 200, -2, 2);
104 ";Truth #eta;Track count / bin", 200, -2, 2);
109 ";Truth Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
113 ";Reco Invariant Mass [GeV/c^2];Pair count / bin", 450, 0, 15);
121 h->GetXaxis()->SetBinLabel(i++,
"Event");
122 h->GetXaxis()->SetBinLabel(i++,
"Truth Track+");
123 h->GetXaxis()->SetBinLabel(i++,
"Truth Track-");
124 h->GetXaxis()->SetBinLabel(i++,
"Reco Track");
125 h->GetXaxis()->SetBinLabel(i++,
"Truth Upsilon");
126 h->GetXaxis()->SetBinLabel(i++,
"Truth Upsilon in Acc.");
127 h->GetXaxis()->SetBinLabel(i++,
"Reco Upsilon");
128 h->GetXaxis()->LabelsOption(
"v");
142 std::cout <<
"QAG4SimulationUpsilon::process_event() entered" << std::endl;
158 assert(h_pTRecoGenRatio);
162 assert(h_pTRecoGenRatio);
166 assert(h_nReco_pTGen);
170 assert(h_nGen_pTGen);
174 assert(h_nReco_etaGen);
178 assert(h_nGen_etaGen);
182 assert(h_nGen_etaGen);
185 assert(h_nGen_etaGen);
190 h_norm->Fill(
"Event", 1);
193 typedef std::set<std::pair<PHG4Particle *, SvtxTrack *>> truth_reco_set_t;
194 truth_reco_set_t truth_reco_set_pos;
195 truth_reco_set_t truth_reco_set_neg;
200 std::cout <<
"QAG4SimulationUpsilon::process_event - fatal error - missing _truthContainer! ";
212 std::cout <<
"QAG4SimulationUpsilon::process_event - processing ";
219 int candidate_embedding_id = trutheval->
get_embed(g4particle);
220 if (candidate_embedding_id < 0) candidate_embedding_id = -1;
226 double gpx = g4particle->
get_px();
227 double gpy = g4particle->
get_py();
228 double gpz = g4particle->
get_pz();
232 if (gpx != 0 && gpy != 0)
234 TVector3 gv(gpx, gpy, gpz);
243 std::cout <<
"QAG4SimulationUpsilon::process_event - accept particle eta = " << geta << std::endl;
249 std::cout <<
"QAG4SimulationUpsilon::process_event - ignore particle eta = " << geta << std::endl;
258 std::cout <<
"QAG4SimulationUpsilon::process_event - ignore particle PID = " << pid <<
" as m_daughterAbsPID = " <<
m_daughterAbsPID << std::endl;
262 TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(pid);
265 std::cout <<
"QAG4SimulationUpsilon::process_event - Error - invalid particle ID = " << pid << std::endl;
269 const double gcharge = pdg_p->Charge() / 3;
272 h_norm->Fill(
"Truth Track+", 1);
274 else if (gcharge < 0)
276 h_norm->Fill(
"Truth Track-", 1);
281 std::cout <<
"QAG4SimulationUpsilon::process_event - invalid neutral decay particle ID = " << pid << std::endl;
286 h_nGen_pTGen->Fill(gpt);
287 h_nGen_etaGen->Fill(geta);
293 h_nReco_etaGen->Fill(geta);
294 h_nReco_pTGen->Fill(gpt);
296 double px = track->
get_px();
297 double py = track->
get_py();
298 double pz = track->
get_pz();
300 TVector3
v(px, py, pz);
305 float pt_ratio = (gpt != 0) ? pt / gpt : 0;
306 h_pTRecoGenRatio->Fill(pt_ratio);
307 h_pTRecoGenRatio_pTGen->Fill(gpt, pt_ratio);
308 h_norm->Fill(
"Reco Track", 1);
313 truth_reco_set_pos.insert(std::make_pair(g4particle, track));
315 else if (gcharge < 0)
317 truth_reco_set_neg.insert(std::make_pair(g4particle, track));
322 TParticlePDG *pdg_p = TDatabasePDG::Instance()->GetParticle(
m_quarkoniaPID);
325 std::cout <<
"QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_quarkoniaPID = " <<
m_quarkoniaPID << std::endl;
329 const double quarkonium_mass = pdg_p->Mass();
330 TParticlePDG *pdg_d = TDatabasePDG::Instance()->GetParticle(
m_daughterAbsPID);
333 std::cout <<
"QAG4SimulationUpsilon::process_event - Fatal Error - invalid particle ID m_daughterAbsPID = " <<
m_daughterAbsPID << std::endl;
337 const double daughter_mass = pdg_d->Mass();
339 for (
const auto &pair_pos : truth_reco_set_pos)
340 for (
const auto &pair_neg : truth_reco_set_neg)
342 assert(pair_pos.first);
343 assert(pair_neg.first);
346 pair_pos.first->get_px(),
347 pair_pos.first->get_py(),
348 pair_pos.first->get_pz(),
349 pair_pos.first->get_e());
352 pair_neg.first->get_px(),
353 pair_neg.first->get_py(),
354 pair_neg.first->get_pz(),
355 pair_neg.first->get_e());
359 if (fabs(quarkonium_mass - gv_quakonium.m()) > 1
e-3)
363 std::cout <<
"QAG4SimulationUpsilon::process_event - invalid pair with in compativle mass with " << quarkonium_mass <<
"GeV for PID = " <<
m_quarkoniaPID <<
": " << std::endl;
364 pair_pos.first->identify();
365 pair_neg.first->identify();
370 h_nGen_Pair_InvMassGen->Fill(gv_quakonium.m());
371 h_norm->Fill(
"Truth Upsilon in Acc.", 1);
373 if (pair_pos.second and pair_neg.second)
380 pair_pos.second->get_px(),
381 pair_pos.second->get_py(),
382 pair_pos.second->get_pz()),
386 pair_neg.second->get_px(),
387 pair_neg.second->get_py(),
388 pair_neg.second->get_pz()),
393 h_nReco_Pair_InvMassReco->Fill(v_quakonium.
m());
394 h_norm->Fill(
"Reco Upsilon", 1);
406 return std::string(
"h_") +
Name() + std::string(
"_");