Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
60 changes: 60 additions & 0 deletions MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
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,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<int> 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;
}
}
Comment thread
wuctlby marked this conversation as resolved.
Outdated

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 +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<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;
}
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
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%

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

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