File indexing completed on 2025-08-05 08:12:43
0001 #include "HFFastSim.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <phool/getClass.h>
0006
0007 #include <phool/PHCompositeNode.h>
0008
0009 #include <phhepmc/PHGenIntegral.h>
0010
0011 #include <TDatabasePDG.h>
0012 #include <TFile.h>
0013 #include <TH1.h>
0014 #include <TH2.h>
0015 #include <TH3.h>
0016 #include <TLorentzVector.h>
0017 #include <TString.h>
0018 #include <TTree.h>
0019 #include <g4jets/Jet.h>
0020 #include <g4jets/JetMap.h>
0021
0022 #include <HepMC/GenEvent.h>
0023 #include <HepMC/GenRanges.h>
0024 #include <HepMC/GenVertex.h>
0025 #include <phhepmc/PHHepMCGenEvent.h>
0026 #include <phhepmc/PHHepMCGenEventMap.h>
0027
0028 #include <gsl/gsl_const.h>
0029
0030 #include <cassert>
0031 #include <cmath>
0032 #include <cstddef>
0033 #include <iostream>
0034
0035 using namespace std;
0036
0037 HFFastSim::HFFastSim(std::string filename, int flavor, int maxevent)
0038 : SubsysReco("HFFastSim")
0039 , _verbose(0)
0040 , _ievent(0)
0041 , _total_pass(0)
0042 , _f(nullptr)
0043 , _h2(nullptr)
0044 , _h2all(nullptr)
0045 , _h2_b(nullptr)
0046 , _h2_c(nullptr)
0047 , m_hNorm(nullptr)
0048 , m_DRapidity(nullptr)
0049 , m_tSingleD(nullptr)
0050 , m_singleD(nullptr)
0051 , m_tSingleLc(nullptr)
0052 , m_singleLc(nullptr)
0053 , _embedding_id(1)
0054 {
0055 _foutname = filename;
0056
0057 _flavor = flavor;
0058
0059 _maxevent = maxevent;
0060
0061 _pt_min = 0;
0062 _pt_max = 100;
0063
0064 _eta_min = -1.1;
0065 _eta_max = +1.1;
0066
0067 _rejection_action = Fun4AllReturnCodes::DISCARDEVENT;
0068 }
0069
0070 int HFFastSim::Init(PHCompositeNode* topNode)
0071 {
0072 _verbose = true;
0073
0074 _ievent = 0;
0075
0076 _total_pass = 0;
0077
0078 _f = new TFile(_foutname.c_str(), "RECREATE");
0079
0080 m_hNorm = new TH1D("hNormalization",
0081 "Normalization;Items;Summed quantity", 10, .5, 10.5);
0082 int i = 1;
0083 m_hNorm->GetXaxis()->SetBinLabel(i++, "IntegratedLumi");
0084 m_hNorm->GetXaxis()->SetBinLabel(i++, "Event");
0085 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0");
0086 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK");
0087 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0-Pair");
0088 m_hNorm->GetXaxis()->SetBinLabel(i++, "D0->PiK-Pair");
0089 m_hNorm->GetXaxis()->SetBinLabel(i++, "Lc");
0090 m_hNorm->GetXaxis()->SetBinLabel(i++, "Lc->pPiK");
0091 m_hNorm->GetXaxis()->SetBinLabel(i++, "Accepted");
0092
0093 m_hNorm->GetXaxis()->LabelsOption("v");
0094
0095 m_DRapidity = new TH2F("hDRapidity",
0096 "hDRapidity;Rapidity of D0 meson;Accepted", 1000, -5, 5, 2, -.5, 1.5);
0097
0098 _h2 = new TH2D("h2", "", 100, 0, 100.0, 40, -2, +2);
0099 _h2_b = new TH2D("h2_b", "", 100, 0, 100.0, 40, -2, +2);
0100 _h2_c = new TH2D("h2_c", "", 100, 0, 100.0, 40, -2, +2);
0101 _h2all = new TH2D("h2all", "", 100, 0, 100.0, 40, -2, +2);
0102
0103 m_tSingleD = new TTree("TSingleD", "TSingleD");
0104
0105 m_singleD = new D02PiK;
0106 m_tSingleD->Branch("D02PiK", "HFFastSim::D02PiK", &m_singleD);
0107
0108 m_tSingleLc = new TTree("TSingleLc", "TSingleLc");
0109
0110 m_singleLc = new Lc2pPiK;
0111 m_tSingleLc->Branch("Lc2pPiK", "HFFastSim::Lc2pPiK", &m_singleLc);
0112
0113 return 0;
0114 }
0115
0116 std::pair<int, TString> HFFastSim::quark_trace(const HepMC::GenParticle* p, HepMC::GenEvent* theEvent)
0117 {
0118 int max_abs_q_pid = 0;
0119 TString chain;
0120
0121 while (p)
0122 {
0123 int pidabs = 0;
0124
0125 TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle((p)->pdg_id());
0126 if (pdg_p)
0127 {
0128 chain = TString(" -> ") + TString(pdg_p->GetName()) + chain;
0129
0130 if (TString(pdg_p->ParticleClass()).BeginsWith("B-"))
0131 {
0132 pidabs = TDatabasePDG::Instance()->GetParticle("b")->PdgCode();
0133 }
0134 else if (TString(pdg_p->ParticleClass()).BeginsWith("Charmed"))
0135 {
0136 pidabs = TDatabasePDG::Instance()->GetParticle("c")->PdgCode();
0137 }
0138 }
0139 else
0140 {
0141 chain = TString(" -> ") + Form("%d", (p)->pdg_id()) + chain;
0142 }
0143
0144 if (pidabs > max_abs_q_pid) max_abs_q_pid = pidabs;
0145
0146 const HepMC::GenVertex* vtx = p->production_vertex();
0147 if (vtx)
0148 {
0149 if (vtx->particles_in_size() == 1)
0150 {
0151
0152 p = *(vtx->particles_in_const_begin());
0153 }
0154 else
0155 {
0156
0157
0158
0159 break;
0160 }
0161 }
0162 }
0163
0164 return make_pair(max_abs_q_pid, chain);
0165 }
0166
0167 bool HFFastSim::process_D02PiK(HepMC::GenEvent* theEvent)
0168 {
0169 assert(theEvent);
0170
0171 const double mom_factor = HepMC::Units::conversion_factor(theEvent->momentum_unit(), HepMC::Units::GEV);
0172 const double length_factor = HepMC::Units::conversion_factor(theEvent->length_unit(), HepMC::Units::CM);
0173
0174
0175 TDatabasePDG* pdg = TDatabasePDG::Instance();
0176
0177 int targetPID = std::abs(pdg->GetParticle("D0")->PdgCode());
0178 int daughter1PID = std::abs(pdg->GetParticle("pi+")->PdgCode());
0179 int daughter2PID = std::abs(pdg->GetParticle("K+")->PdgCode());
0180
0181 unsigned int nD0(0);
0182 unsigned int nD0PiK(0);
0183
0184 auto range = theEvent->particle_range();
0185 for (HepMC::GenEvent::particle_const_iterator piter = range.begin(); piter != range.end(); ++piter)
0186 {
0187 const HepMC::GenParticle* p = *piter;
0188 assert(p);
0189
0190 if (std::abs(p->pdg_id()) == targetPID)
0191 {
0192 if (Verbosity() >= VERBOSITY_MORE)
0193 {
0194 cout << "process_D02PiK - Accept signal particle : ";
0195 p->print();
0196 cout << endl;
0197 }
0198
0199 m_hNorm->Fill("D0", 1);
0200 ++nD0;
0201
0202 assert(m_DRapidity);
0203 const double rapidity = 0.5 * log((p->momentum().e() + p->momentum().z()) /
0204 (p->momentum().e() - p->momentum().z()));
0205
0206 m_DRapidity->Fill(rapidity, 0);
0207
0208 const HepMC::GenVertex* decayVertex = p->end_vertex();
0209
0210 int hasDecay1(0);
0211 int hasDecay2(0);
0212 int hasDecayOther(0);
0213
0214 if (decayVertex)
0215 {
0216 for (auto diter = decayVertex->particles_out_const_begin();
0217 diter != decayVertex->particles_out_const_end();
0218 ++diter)
0219
0220 {
0221 const HepMC::GenParticle* pd = *diter;
0222 assert(pd);
0223
0224 if (Verbosity() >= VERBOSITY_MORE)
0225 {
0226 cout << "process_D02PiK - Testing daughter particle: ";
0227 pd->print();
0228 cout << endl;
0229 }
0230
0231 if (std::abs(pd->pdg_id()) == daughter1PID)
0232 {
0233 const double eta = pd->momentum().eta();
0234
0235 if (eta > _eta_min and eta < _eta_max)
0236 ++hasDecay1;
0237 }
0238 else if (std::abs(pd->pdg_id()) == daughter2PID)
0239 {
0240 const double eta = pd->momentum().eta();
0241
0242 if (eta > _eta_min and eta < _eta_max)
0243 ++hasDecay2;
0244 }
0245 else
0246 ++hasDecayOther;
0247 }
0248
0249
0250 if (hasDecay1 == 1 and hasDecay2 == 1 and hasDecayOther == 0)
0251 {
0252 m_hNorm->Fill("D0->PiK", 1);
0253 ++nD0PiK;
0254
0255 assert(m_singleD);
0256 m_singleD->pid = p->pdg_id();
0257 m_singleD->y = rapidity;
0258 m_singleD->pt = p->momentum().perp() * mom_factor;
0259
0260 m_singleD->decayL = decayVertex->point3d().r() * length_factor;
0261 m_singleD->prodL = p->production_vertex()->point3d().r() * length_factor;
0262
0263 auto trace = quark_trace(p, theEvent);
0264 m_singleD->hadron_origion_qpid = trace.first;
0265 m_singleD->decayChain = trace.second;
0266
0267 if (Verbosity())
0268 {
0269 cout << "process_D02PiK - Found D0->piK:";
0270 cout << "m_singleD->pid = " << m_singleD->pid << ", ";
0271 cout << "m_singleD->y = " << m_singleD->y << ", ";
0272 cout << "m_singleD->pt = " << m_singleD->pt;
0273 cout << " origin " << m_singleD->decayChain << endl;
0274 }
0275
0276 assert(m_tSingleD);
0277 m_tSingleD->Fill();
0278 m_singleD->Clear();
0279
0280 m_DRapidity->Fill(rapidity, 1);
0281 }
0282
0283 }
0284 else
0285 {
0286 cout << "process_D02PiK - Warning - target particle did not have decay vertex : ";
0287 p->print();
0288 cout << endl;
0289 }
0290
0291 }
0292 }
0293
0294 if (nD0 >= 2)
0295 {
0296 if (Verbosity() >= VERBOSITY_MORE)
0297 cout << "process_D02PiK - D0-Pair with nD0 = " << nD0 << endl;
0298 m_hNorm->Fill("D0-Pair", nD0 * (nD0 - 1) / 2);
0299 }
0300 if (nD0PiK >= 2)
0301 {
0302 m_hNorm->Fill("D0->PiK-Pair", nD0PiK * (nD0PiK - 1) / 2);
0303 }
0304
0305 return nD0PiK >= 1;
0306 }
0307
0308 bool HFFastSim::process_Lc2pPiK(HepMC::GenEvent* theEvent)
0309 {
0310 assert(theEvent);
0311
0312 const double mom_factor = HepMC::Units::conversion_factor(theEvent->momentum_unit(), HepMC::Units::GEV);
0313 const double length_factor = HepMC::Units::conversion_factor(theEvent->length_unit(), HepMC::Units::CM);
0314
0315
0316 TDatabasePDG* pdg = TDatabasePDG::Instance();
0317
0318 const int targetPID = std::abs(pdg->GetParticle("Lambda_c+")->PdgCode());
0319 assert(targetPID == 4122);
0320 const int daughter1PID = std::abs(pdg->GetParticle("pi+")->PdgCode());
0321 assert(daughter1PID == 211);
0322 const int daughter2PID = std::abs(pdg->GetParticle("K+")->PdgCode());
0323 assert(daughter2PID == 321);
0324 const int daughter3PID = std::abs(pdg->GetParticle("proton")->PdgCode());
0325 assert(daughter3PID == 2212);
0326
0327 unsigned int nLc(0);
0328 unsigned int nLc2pPiK(0);
0329
0330 auto range = theEvent->particle_range();
0331 for (HepMC::GenEvent::particle_const_iterator piter = range.begin(); piter != range.end(); ++piter)
0332 {
0333 const HepMC::GenParticle* p = *piter;
0334 assert(p);
0335
0336 if (std::abs(p->pdg_id()) == targetPID)
0337 {
0338 if (Verbosity() >= VERBOSITY_MORE)
0339 {
0340 cout << "HFFastSim::process_Lc2pPiK - Accept signal particle : ";
0341 p->print();
0342 cout << endl;
0343 }
0344
0345 m_hNorm->Fill("Lc", 1);
0346 ++nLc;
0347
0348 const double rapidity = 0.5 * log((p->momentum().e() + p->momentum().z()) /
0349 (p->momentum().e() - p->momentum().z()));
0350
0351 const HepMC::GenVertex* decayVertex = p->end_vertex();
0352
0353 int hasDecay1(0);
0354 int hasDecay2(0);
0355 int hasDecay3(0);
0356 int hasDecayOther(0);
0357
0358 if (decayVertex)
0359 {
0360 for (auto diter = decayVertex->particles_out_const_begin();
0361 diter != decayVertex->particles_out_const_end();
0362 ++diter)
0363
0364 {
0365 const HepMC::GenParticle* pd = *diter;
0366 assert(pd);
0367
0368 if (Verbosity() >= VERBOSITY_MORE)
0369 {
0370 cout << "HFFastSim::process_Lc2pPiK - Testing daughter particle: ";
0371 pd->print();
0372 cout << endl;
0373 }
0374
0375 if (std::abs(pd->pdg_id()) == daughter1PID)
0376 {
0377 const double eta = pd->momentum().eta();
0378
0379 if (eta > _eta_min and eta < _eta_max)
0380 ++hasDecay1;
0381 }
0382 else if (std::abs(pd->pdg_id()) == daughter2PID)
0383 {
0384 const double eta = pd->momentum().eta();
0385
0386 if (eta > _eta_min and eta < _eta_max)
0387 ++hasDecay2;
0388 }
0389 else if (std::abs(pd->pdg_id()) == daughter3PID)
0390 {
0391 const double eta = pd->momentum().eta();
0392
0393 if (eta > _eta_min and eta < _eta_max)
0394 ++hasDecay3;
0395 }
0396 else
0397 ++hasDecayOther;
0398 }
0399
0400
0401 if (hasDecay1 == 1 and hasDecay2 == 1 and hasDecay3 == 1 and hasDecayOther == 0)
0402 {
0403 m_hNorm->Fill("Lc->pPiK", 1);
0404 ++nLc2pPiK;
0405
0406 assert(m_singleLc);
0407 m_singleLc->pid = p->pdg_id();
0408 m_singleLc->y = rapidity;
0409 m_singleLc->pt = p->momentum().perp() * mom_factor;
0410
0411 m_singleLc->decayL = decayVertex->point3d().r() * length_factor;
0412 m_singleLc->prodL = p->production_vertex()->point3d().r() * length_factor;
0413
0414 auto trace = quark_trace(p, theEvent);
0415 m_singleLc->hadron_origion_qpid = trace.first;
0416 m_singleLc->decayChain = trace.second;
0417
0418 if (Verbosity())
0419 {
0420 cout << "HFFastSim::process_Lc2pPiK - Found D0->piK:";
0421 cout << "m_singleLc->pid = " << m_singleLc->pid << ", ";
0422 cout << "m_singleLc->y = " << m_singleLc->y << ", ";
0423 cout << "m_singleLc->pt = " << m_singleLc->pt;
0424 cout << " origin " << m_singleLc->decayChain << endl;
0425 }
0426
0427 assert(m_tSingleLc);
0428 m_tSingleLc->Fill();
0429 m_singleLc->Clear();
0430
0431 m_DRapidity->Fill(rapidity, 1);
0432 }
0433
0434 }
0435 else
0436 {
0437 cout << "HFFastSim::process_Lc2pPiK - Warning - target particle did not have decay vertex : ";
0438 p->print();
0439 cout << endl;
0440 }
0441
0442 }
0443 }
0444
0445 return nLc2pPiK >= 1;
0446 }
0447
0448 int HFFastSim::process_event(PHCompositeNode* topNode)
0449 {
0450 assert(m_hNorm);
0451 m_hNorm->Fill("Event", 1);
0452
0453 PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0454 if (!geneventmap)
0455 {
0456 std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
0457 return Fun4AllReturnCodes::ABORTRUN;
0458 }
0459
0460 PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
0461 if (!genevt)
0462 {
0463 std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
0464 std::cout << ". Print PHHepMCGenEventMap:";
0465 geneventmap->identify();
0466 return Fun4AllReturnCodes::ABORTRUN;
0467 }
0468
0469 HepMC::GenEvent* theEvent = genevt->getEvent();
0470
0471 assert(theEvent);
0472
0473 if (Verbosity() >= VERBOSITY_MORE)
0474 {
0475 cout << "HFFastSim::process_event - process HepMC::GenEvent with signal_process_id = "
0476 << theEvent->signal_process_id();
0477 if (theEvent->signal_process_vertex())
0478 {
0479 cout << " and signal_process_vertex : ";
0480 theEvent->signal_process_vertex()->print();
0481 }
0482 cout << " - Event record:" << endl;
0483 theEvent->print();
0484 }
0485
0486
0487 for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
0488 p != theEvent->particles_end(); ++p)
0489 {
0490 int pidabs = abs((*p)->pdg_id());
0491 if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode()
0492 or pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode()
0493 or pidabs == TDatabasePDG::Instance()->GetParticle("t")->PdgCode())
0494 {
0495 if (pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
0496 {
0497 if (Verbosity() >= VERBOSITY_MORE)
0498 std::cout << __PRETTY_FUNCTION__
0499 << " --BOTTOM--> pt / eta / phi = "
0500 << ((*p)->momentum().perp()) << " / "
0501 << (*p)->momentum().pseudoRapidity() << " / "
0502 << (*p)->momentum().phi() << std::endl;
0503
0504 _h2_b->Fill((*p)->momentum().perp(), (*p)->momentum().pseudoRapidity());
0505 }
0506 else if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
0507 {
0508 if (Verbosity() >= VERBOSITY_MORE)
0509 std::cout << __PRETTY_FUNCTION__
0510 << " --CHARM --> pt / eta / phi = "
0511 << (*p)->momentum().perp() << " / "
0512 << (*p)->momentum().pseudoRapidity() << " / "
0513 << (*p)->momentum().phi() << std::endl;
0514
0515 _h2_c->Fill((*p)->momentum().perp(), (*p)->momentum().pseudoRapidity());
0516 }
0517 }
0518
0519 }
0520
0521 bool pass_event = process_D02PiK(theEvent);
0522 process_Lc2pPiK(theEvent);
0523
0524 if (pass_event && _total_pass >= _maxevent)
0525 {
0526 if (Verbosity() >= HFFastSim::VERBOSITY_MORE)
0527 std::cout << __PRETTY_FUNCTION__ << " --> FAIL due to max events = "
0528 << _total_pass << std::endl;
0529 _ievent++;
0530 return _rejection_action;
0531 }
0532 else if (pass_event)
0533 {
0534 _total_pass++;
0535 if (Verbosity() >= HFFastSim::VERBOSITY_MORE)
0536 std::cout << __PRETTY_FUNCTION__ << " --> PASS, total = " << _total_pass
0537 << std::endl;
0538 _ievent++;
0539 m_hNorm->Fill("Accepted", 1);
0540 return Fun4AllReturnCodes::EVENT_OK;
0541 }
0542 else
0543 {
0544 if (Verbosity() >= HFFastSim::VERBOSITY_MORE)
0545 std::cout << __PRETTY_FUNCTION__ << " --> FAIL " << std::endl;
0546 _ievent++;
0547 return _rejection_action;
0548 }
0549 }
0550
0551 int HFFastSim::End(PHCompositeNode* topNode)
0552 {
0553 if (Verbosity() >= HFFastSim::VERBOSITY_SOME)
0554 std::cout << __PRETTY_FUNCTION__ << " PASSED " << _total_pass
0555 << " events" << std::endl;
0556
0557 PHGenIntegral* integral_node = findNode::getClass<PHGenIntegral>(topNode, "PHGenIntegral");
0558 if (integral_node)
0559 {
0560 integral_node->identify();
0561 assert(m_hNorm);
0562 m_hNorm->Fill("IntegratedLumi", integral_node->get_Integrated_Lumi());
0563 }
0564
0565 _f->cd();
0566 _f->Write();
0567
0568
0569 return 0;
0570 }
