From 5f0a79c5a15970e476f689412300b391ade8b1b0 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Thu, 20 Aug 2026 13:57:15 +0800 Subject: [PATCH 1/3] tunable ccbar and bbbar --- .../generator/generator_pythia8_embed_hf.C | 60 +++++++++++++++++++ .../GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini | 9 +++ ...adronic_with_decays_Mode2_hardQCD_5TeV.cfg | 4 +- 3 files changed, 71 insertions(+), 2 deletions(-) create mode 100644 MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini 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 7a28d631a..c6284fffd 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C @@ -18,6 +18,7 @@ namespace hf_generators GapTriggeredBeauty, // --> GeneratorPythia8GapTriggeredBeauty: beauty enriched GapTriggeredCharmAndBeauty, // --> GeneratorPythia8GapTriggeredCharmAndBeauty: charm and beauty enriched (with same ratio) GapHF, // --> GeneratorPythia8GapHF + GapHFRatio, // --> GeneratorPythia8GapHF, quark list built from b/c ratio NGenType }; } @@ -83,6 +84,12 @@ public: LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs; mGeneratorEvHF = dynamic_cast(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList)); break; + + case hf_generators::GapHFRatio: + LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GapHFRatio (custom b/c ratio) **********"; + LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs; + mGeneratorEvHF = dynamic_cast(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList)); + break; default: LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********"; break; @@ -387,6 +394,44 @@ private: }; +// Helper: build quarkPdgList from bOverCRatio (= N(beauty)/N(charm), default 1 = 1:1) +// Integer b/c: taken directly; fractional b/c: reduced to nearest nB:nC with nC+nB <= 20 +static std::vector BuildQuarkListFromBOverC(float bOverCRatio) +{ + if (bOverCRatio <= 0.f) { + LOG(fatal) << "bOverCRatio (b/c) must be > 0"; + } + if (bOverCRatio > 19.f) { + bOverCRatio = 19.f; + LOG(warn) << "bOverCRatio (b/c) too large, using 19:1"; + }else if (bOverCRatio < 1.f/19.f) { + bOverCRatio = 1.f/19.f; + LOG(warn) << "bOverCRatio (b/c) too small, using 1:19"; + } + + int nC = 1, nB = 1; + float bestErr = 1e9f; + for (int c = 1; c <= 19; ++c) { + int b = (int)std::lround(bOverCRatio * c); + if (b < 1 || c + b > 20) { + continue; + } + float err = std::fabs((float)b / c - bOverCRatio); + if (err < bestErr) { + bestErr = err; + nC = c; + nB = b; + } + } + + std::vector list; + // Bresenham interleaving + int len = nC + nB; + for (int k = 0; k < len; ++k) + list.push_back((((k + 1) * nC) / len > (k * nC) / len) ? 4 : 5); + return list; +} + // Charm enriched FairGenerator * GeneratorPythia8EmbedHFCharm(bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector quarkPdgList = {}, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector freqReplaceList = {}) { @@ -420,3 +465,18 @@ FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false return myGen; } +// Charm and beauty enriched with tunable b/c ratio +FairGenerator * GeneratorPythia8EmbedHFRatio(float bOverCRatio = 1.f, bool usePtHardBins = false, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector hadronPdgList = {}, std::vector> partPdgToReplaceList = {}, std::vector freqReplaceList = {}) +{ + auto myGen = new GeneratorPythia8EmbedHF(); + + /// build the quark list from the b/c ratio + auto quarkPdgList = BuildQuarkListFromBOverC(bOverCRatio); + + /// setup the internal generator for HF events + myGen->setupGeneratorEvHF(hf_generators::GapHFRatio, + usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, + quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList); + + return myGen; +} diff --git a/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini new file mode 100644 index 000000000..418304b17 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini @@ -0,0 +1,9 @@ +#NEV_TEST> 10 +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C +funcName=GeneratorPythia8EmbedHFRatio(3.0) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg +includePartonEvent=true diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg index 81da1a4ba..f407097af 100644 --- a/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg +++ b/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg @@ -107,7 +107,7 @@ BeamRemnants:saturation 5 4122:addChannel = 1 0.04500 100 2224 -321 ### Λc+ -> Delta++ K- 1.08% 4122:addChannel = 1 0.09000 100 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3 ### Λc+ -> p K0S (36%) -4122:addChannel = 1 0.36000 0 2212 311 ### Λc+ -> p K0S 1.59% +4122:addChannel = 1 0.36000 0 2212 310 ### Λc+ -> p K0S 1.59% ### Λc+ -> p K- π+ π0 (small, 3%) 4122:addChannel = 1 0.03000 0 2212 -321 211 111 ### Λc+ -> p K- π+ π0 (non-resonant) 4.6% ### Λc+ -> p π- π+ (12.50%) @@ -256,7 +256,7 @@ BeamRemnants:saturation 5 431:onIfMatch = 221 211 ### Λc -> pK0s -4122:onIfMatch = 2212 311 +4122:onIfMatch = 2212 310 ### Λc -> p K- π+ π0 4122:onIfMatch = 2212 321 211 ### Λc -> p K* From f9f9e7f20ebe1b3ba27105659852d79a77fdd3f5 Mon Sep 17 00:00:00 2001 From: Chuntai <48704924+wuctlby@users.noreply.github.com> Date: Fri, 21 Aug 2026 17:23:17 +0800 Subject: [PATCH 2/3] Update MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C Co-authored-by: Mattia Faggin --- .../generator/generator_pythia8_embed_hf.C | 77 +++++++++++++------ 1 file changed, 52 insertions(+), 25 deletions(-) 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 c6284fffd..7f325eb67 100644 --- a/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C @@ -398,31 +398,58 @@ private: // Integer b/c: taken directly; fractional b/c: reduced to nearest nB:nC with nC+nB <= 20 static std::vector BuildQuarkListFromBOverC(float bOverCRatio) { - if (bOverCRatio <= 0.f) { - LOG(fatal) << "bOverCRatio (b/c) must be > 0"; - } - if (bOverCRatio > 19.f) { - bOverCRatio = 19.f; - LOG(warn) << "bOverCRatio (b/c) too large, using 19:1"; - }else if (bOverCRatio < 1.f/19.f) { - bOverCRatio = 1.f/19.f; - LOG(warn) << "bOverCRatio (b/c) too small, using 1:19"; - } - - int nC = 1, nB = 1; - float bestErr = 1e9f; - for (int c = 1; c <= 19; ++c) { - int b = (int)std::lround(bOverCRatio * c); - if (b < 1 || c + b > 20) { - continue; - } - float err = std::fabs((float)b / c - bOverCRatio); - if (err < bestErr) { - bestErr = err; - nC = c; - nB = b; - } - } +const int iterNMax = 19; + if (bOverCRatio <= 0.f) { + LOG(fatal) << "bOverCRatio (b/c) must be > 0"; + } + if (bOverCRatio > iterNMax*1.f) { + bOverCRatio = iterNMax*1.f; + LOG(warn) << "bOverCRatio (b/c) too large, using 19:1"; + }else if (bOverCRatio < 1.f/iterNMax) { + bOverCRatio = 1.f/iterNMax; + LOG(warn) << "bOverCRatio (b/c) too small, using 1:19"; + } + + int nC = 1, nB = 1; + float bestErr = 1e9f; + for (int c = 1; c <= iterNMax; ++c) { + int b = static_cast(std::lround(bOverCRatio * c)); + if (b < 1 || c + b > 20) { + continue; + } + float err = std::fabs(static_cast(b) / c - bOverCRatio); + if (err < bestErr) { + /// This check here is needed in case bOverCRatio*c is not integer (it can happen with e.g. bOverCRatio=4./9.) + /// In this case, b is its truncation and the desired ratio is not obtained + /// It means that one needs to continue iterating until bOverCRatio*c is integer, i.e. b not truncated + /// + /// Possible cases: + /// 1. we want nB = nC*R, with R integer + /// -> we enter here in the first loop iteration + /// 2. we want more charm than beauty by an integer amount, i.e. nB = nC*R with R=1./N, with N integer + /// -> we enter here after N iterations, when c=N + /// 3. we want either more charm or beauty, but with a factor R that is not integer, as well as its inverse (e.g. bOverCRatio=4./9.) + /// -> In this case, bOverCRatio*c is not integer, namely b is its truncation and the desired ratio is not obtained + /// The code iterates at most until c becomes equal to the denominator of the fraction + /// 4. we want one of the previous cases, but we assign to bOverCRatio a value that is not rational + /// or such as we do not enter here within 19 iterations + /// -> nC and nB, are not touched, therefore we do not have the desired fraction. We'll need to throw a fatal (*) + /// + bestErr = err; + nC = c; + nB = b; + + /// If we are at this point, we reached already the desired ratio between b and c. + /// Let's break the loop + break; + } + } + + // (*) check if we have the desired fraction + bool isRatioUnity = std::fabs(bOverCRatio-1) < 1e-05; + if (!isRatioUnity && nC==1 && nB==1) { + LOG(fatal) << "nC=" << nC << ", nB=" << nB << " but you ask bOverCRatio to be " << bOverCRatio <<", which is different from nB/nC. It means either that bOverCRatio is not rational, or that you need more than " << iterNMax << " iterations. Change it!"; + } std::vector list; // Bresenham interleaving From 1e3f8b2c764f0e774b44270b69fd4266d7fbbd69 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Sun, 23 Aug 2026 14:50:43 +0800 Subject: [PATCH 3/3] Add test macro --- .../GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C | 192 ++++++++++++++++++ 1 file changed, 192 insertions(+) create mode 100644 MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C new file mode 100644 index 000000000..c73693d87 --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C @@ -0,0 +1,192 @@ +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 + + std::vector checkPdgHadron{411, 421, 431, 4122, 4132, 4232, 4332}; + std::map>> checkHadronDecays{ // sorted pdg of daughters + {421, { + {-321, 211}, // D0 -> K-, pi+ + {-321, 211, 111}, // D0 -> K-, pi+, pi0 + {213, -321}, // D0 -> rho(770)+, K- + {-313, 111}, // D0 -> Kbar^*(892)0, pi0 + {-323, 211}, // D0 -> K^*(892)-, pi+ + {-211, 211}, // D0 -> pi-, pi+ + {213, -211}, // D0 -> rho(770)+, pi- + {-211, 211, 111}, // D0 -> pi-, pi+, pi0 + {-321, 321}, // D0 -> K-, K+ + }}, + + {411, { + {-321, 211, 211}, // D+ -> K-, pi+, pi+ + {-10311, 211}, // D+ -> Kbar0^*(1430)0, pi+ + {-313, 211}, // D+ -> Kbar^*(892)0, pi+ + {-321, 211, 211, 111}, // D+ -> K-, pi+, pi+, pi0 + {333, 211}, // D+ -> phi(1020)0, pi+ + {-313, 321}, // D+ -> Kbar^*(892)0, K+ + {-10311, 321}, // D+ -> Kbar0^*(1430)0, K+ + {-321, 321, 211}, // D+ -> K-, K+, pi+ + {113, 211}, // D+ -> rho(770)0, pi+ + {225, 211}, // D+ -> f2(1270)0, pi+ + {-211, 211, 211}, // D+ -> pi-, pi+, pi+ + }}, + + {431, { + {333, 211}, // Ds+ -> phi(1020)0, pi+ + {-313, 321}, // Ds+ -> Kbar^*(892)0, K+ + {333, 213}, // Ds+ -> phi(1020)0, rho(770)+ + {113, 211}, // Ds+ -> rho(770)0, pi+ + {225, 211}, // Ds+ -> f2(1270)0, pi+ + {-211, 211, 211}, // Ds+ -> pi-, pi+, pi+ + {313, 211}, // Ds+ -> K^*(892)0, pi+ + {10221, 321}, // Ds+ -> f0(1370)0, K+ + {113, 321}, // Ds+ -> rho(770)0, K+ + {-211, 321, 211}, // Ds+ -> pi-, K+, pi+ + {221, 211}, // Ds+ -> eta, pi+ + }}, + + {4122, { + {2212, -321, 211}, // Lambdac+ -> p, K-, pi+ + {2212, -313}, // Lambdac+ -> p, Kbar^*(892)0 + {2224, -321}, // Lambdac+ -> Delta(1232)++, K- + {102134, 211}, // Lambdac+ -> 102134, pi+ + {2212, 310}, // Lambdac+ -> p, K0s + {2212, -321, 211, 111}, // Lambdac+ -> p, K-, pi+, pi0 + {2212, -211, 211}, // Lambdac+ -> p, pi-, pi+ + {2212, 333}, // Lambdac+ -> p, phi(1020)0 + }}, + + {4232, { + {2212, -321, 211}, // Xic+ -> p, K-, pi+ + {2212, -313}, // Xic+ -> p, Kbar^*(892)0 + {3312, 211, 211}, // Xic+ -> Xi-, pi+, pi+ + {2212, 333}, // Xic+ -> p, phi(1020)0 + {3222, -211, 211}, // Xic+ -> Sigma+, pi-, pi+ + {3324, 211}, // Xic+ -> Xi(1530)0, pi+ + }}, + + {4132, { + {3312, 211}, // Xic0 -> Xi-, pi+ + }}, + + {4332, { + {3334, 211}, // Omegac0 -> Omega-, pi+ + {3312, 211}, // Omegac0 -> Xi-, pi+ + }}, + }; + + 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 + + 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 && pdgDau != 111 && pdgDau != 221 && pdgDau != 113 && pdgDau != 225) { // 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)]) { + std::sort(decay.begin(), decay.end()); + if (pdgsDecay == decay || pdgsDecayAntiPart == decay) { + nSignalGoodDecay++; + break; + } + } + } + } + } + + 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"; + + //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 = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f; + float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 1.0f; + if (1 - fracForcedDecays > 0.15 + uncFracForcedDecays) { // 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; +}