From 5f0a79c5a15970e476f689412300b391ade8b1b0 Mon Sep 17 00:00:00 2001 From: wulbyu Date: Thu, 20 Aug 2026 13:57:15 +0800 Subject: [PATCH 1/2] 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/2] 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