diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C index d45e4ad72..1fd07e182 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C @@ -45,12 +45,13 @@ public: /// @brief setup the event generator for HF signals /// \param gentype generator type (only ccbar, only bbbar, both) + /// \param usePtHardBins flag to enable/disable pt-hard bins /// \param yQuarkMin minimum quark rapidity /// \param yQuarkMax maximum quark rapidity /// \param yHadronMin minimum hadron rapidity /// \param yHadronMax maximum hadron rapidity /// \param hadronPdgList list of PDG codes for hadrons to be used in trigger - void setupGeneratorEvHF(int genType, float yQuarkMin, float yQuarkMax, float yHadronMin, float yHadronMax, std::vector hadronPdgList = {}) { + void setupGeneratorEvHF(int genType, bool usePtHardBins, float yQuarkMin, float yQuarkMax, float yHadronMin, float yHadronMax, std::vector hadronPdgList = {}) { mGeneratorEvHF = nullptr; switch (genType) { @@ -80,18 +81,20 @@ public: } // we set pT hard bins - auto seed = dynamic_cast(mGeneratorEvHF)->getUsedSeed(); - float ptHardBins[4] = {2.76, 20., 50., 1000.}; - int iPt{0}; - if (seed % 10 < 7) { - iPt = 0; - } else if (seed % 10 < 9) { - iPt = 1; - } else { - iPt = 2; + if (usePtHardBins) { + auto seed = dynamic_cast(mGeneratorEvHF)->getUsedSeed(); + float ptHardBins[4] = {2.76, 20., 50., 1000.}; + int iPt{0}; + if (seed % 10 < 7) { + iPt = 0; + } else if (seed % 10 < 9) { + iPt = 1; + } else { + iPt = 2; + } + dynamic_cast(mGeneratorEvHF)->readString(Form("PhaseSpace:pTHatMin = %f", ptHardBins[iPt])); + dynamic_cast(mGeneratorEvHF)->readString(Form("PhaseSpace:pTHatMax = %f", ptHardBins[iPt+1])); } - dynamic_cast(mGeneratorEvHF)->readString(Form("PhaseSpace:pTHatMin = %f", ptHardBins[iPt])); - dynamic_cast(mGeneratorEvHF)->readString(Form("PhaseSpace:pTHatMax = %f", ptHardBins[iPt+1])); mGeneratorEvHF->Init(); } @@ -377,34 +380,34 @@ private: }; // Charm enriched -FairGenerator * GeneratorPythia8EmbedHFCharm(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFCharm(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) { auto myGen = new GeneratorPythia8EmbedHF(); /// setup the internal generator for HF events - myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharm, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); + myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharm, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); return myGen; } // Beauty enriched -FairGenerator * GeneratorPythia8EmbedHFBeauty(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) { auto myGen = new GeneratorPythia8EmbedHF(); /// setup the internal generator for HF events - myGen->setupGeneratorEvHF(hf_generators::GapTriggeredBeauty, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); + myGen->setupGeneratorEvHF(hf_generators::GapTriggeredBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); return myGen; } // Charm and beauty enriched (with same ratio) -FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) +FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}) { auto myGen = new GeneratorPythia8EmbedHF(); /// setup the internal generator for HF events - myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharmAndBeauty, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); + myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharmAndBeauty, usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList); return myGen; } diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.ini new file mode 100644 index 000000000..e90890584 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.ini @@ -0,0 +1,8 @@ +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C +funcName=GeneratorPythia8EmbedHFCharmAndBeauty(true) + +[GeneratorPythia8] +config=${O2DPG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg +includePartonEvent=true diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C new file mode 100644 index 000000000..112cda1f4 --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_ptHardBins.C @@ -0,0 +1,134 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + //std::string path{"tf1/sgn_1_Kine.root"}; + + int checkPdgQuarkOne{4}; + int checkPdgQuarkTwo{5}; + float ratioTrigger = 1.; // one event triggered out of 1 + float averagePt = 0.; + + std::vector checkPdgHadron{411, 421, 431, 4122, 4132, 4232, 4332}; + std::map>> checkHadronDecays{ // sorted pdg of daughters + {411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+ + {421, {{-321, 211}, {-321, 111, 211}}}, // D0 + {431, {{211, 333}, {-313, 321}}}, // Ds+ + {4122, {{-313, 2212}, {-321, 2224}, {211, 102134}, {-321, 211, 2212}, {311, 2212}}}, // Lc+ + {4132, {{211, 3312}}}, // Xic0 + {4232, {{-313, 2212}, {-321, 3324}, {211, 211, 3312}, {-321, 211, 2212}}}, // Xic+ + {4332, {{211, 3334}}} // Omegac+ + }; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + o2::dataformats::MCEventHeader *eventHeader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + + int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{}; + int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{}; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + + // check subgenerator information + //if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) { + // bool isValid = false; + // int subGeneratorId = eventHeader->getInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid); + // if (subGeneratorId == 0) { + // nEventsMB++; + // } else if (subGeneratorId == checkPdgQuarkOne) { + // nEventsInjOne++; + // } else if (subGeneratorId == checkPdgQuarkTwo) { + // nEventsInjTwo++; + // } + //} + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + if (std::abs(pdg) == checkPdgQuarkOne) { + nQuarksOne++; + continue; + } + if (std::abs(pdg) == checkPdgQuarkTwo) { + nQuarksTwo++; + continue; + } + if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal + nSignals++; // count signal PDG + averagePt += track.GetPt(); + + std::vector pdgsDecay{}; + std::vector pdgsDecayAntiPart{}; + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + auto pdgDau = tracks->at(j).GetPdgCode(); + pdgsDecay.push_back(pdgDau); + if (pdgDau != 333) { // phi is antiparticle of itself + pdgsDecayAntiPart.push_back(-pdgDau); + } else { + pdgsDecayAntiPart.push_back(pdgDau); + } + } + + std::sort(pdgsDecay.begin(), pdgsDecay.end()); + std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end()); + + for (auto &decay : checkHadronDecays[std::abs(pdg)]) { + if (pdgsDecay == decay || pdgsDecayAntiPart == decay) { + nSignalGoodDecay++; + break; + } + } + } + } + } + + averagePt /= nSignals; + + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + //std::cout << "# MB events: " << nEventsMB << "\n"; + //std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n"; + //std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n"; + std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n"; + std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n"; + std::cout <<"# signal hadrons: " << nSignals << "\n"; + std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n"; + std::cout <<"average pT of signal hadrons: " << averagePt << "\n"; + + //if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small + // std::cerr << "Number of generated MB events different than expected\n"; + // return 1; + //} + //if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) { + // std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n"; + // return 1; + //} + //if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) { + // std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n"; + // return 1; + //} + + if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation + std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n"; + return 1; + } + if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation + std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n"; + return 1; + } + + float fracForcedDecays = float(nSignalGoodDecay) / nSignals; + if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state) + std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n"; + return 1; + } + + return 0; +}