Skip to content

tunable ratio of ccbar and bbbar for embedding #2435

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wuctlby wants to merge 3 commits into
base: master
Choose a base branch
from
Open

tunable ratio of ccbar and bbbar for embedding #2435

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5f0a79c
tunable ccbar and bbbar
wuctlby Aug 20, 2026
f9f9e7f
Update MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
wuctlby Aug 21, 2026
1e3f8b2
Add test macro
wuctlby Aug 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 87 additions & 0 deletions MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
};
}
Expand Down Expand Up @@ -83,6 +84,12 @@ public:
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(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<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList));
break;
default:
LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********";
break;
Expand Down Expand Up @@ -387,6 +394,71 @@ 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<int> BuildQuarkListFromBOverC(float bOverCRatio)
{
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<int>(std::lround(bOverCRatio * c));
if (b < 1 || c + b > 20) {
continue;
}
float err = std::fabs(static_cast<float>(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<int> 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<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
{
Expand Down Expand Up @@ -420,3 +492,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<int> hadronPdgList = {}, std::vector<std::array<int,2>> partPdgToReplaceList = {}, std::vector<float> 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;
}
9 changes: 9 additions & 0 deletions MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
192 changes: 192 additions & 0 deletions MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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<int> checkPdgHadron{411, 421, 431, 4122, 4132, 4232, 4332};
std::map<int, std::vector<std::vector<int>>> 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<o2::MCTrack> *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<int>(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<int> pdgsDecay{};
std::vector<int> 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;
}
4 changes: 2 additions & 2 deletions MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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%

@mfaggin mfaggin Aug 20, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In principle this should work in pythia, but did you check that this works?
The physics wants that only 4122:addChannel = 1 0.36000 0 2212 311 is correct, since K0, K0bar are the only states with a well-defined strangeness content, that materialize as K0s or K0L only in decay.
Did you check that the correct and expected decays are present, i.e. the fact that we have Lc decaying into p and K0s with K0s decaying in 2 pions?

@wuctlby wuctlby Aug 21, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I did a private test, and it works.
I checked the decay channels produced in one of the kinematic files, and they look good.
ScreenShot_2026-08-21_184841_592
Then I further checked the decay of produced K0s, as follows:
Found K(S)0 daughters for Lambdac+ decay: 111 (#pi^{0}) 111 (#pi^{0})
Found K(S)0 daughters for Lambdac+ decay: 211 (#pi^{+}) -211 (#pi^{-})

### Λ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%)
Expand Down Expand Up @@ -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*
Expand Down
Loading