Skip to content

Commit a7854c2

Browse files
wuctlbymfaggin
wuctlby
and
mfaggin
authored
tunable ratio of ccbar and bbbar for embedding (#2435)
* tunable ccbar and bbbar * Update MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C Co-authored-by: Mattia Faggin <mattia.faggin@cern.ch> * Add test macro --------- Co-authored-by: Mattia Faggin <mattia.faggin@cern.ch>
1 parent a8c42f8 commit a7854c2

4 files changed

Lines changed: 290 additions & 2 deletions

File tree

‎MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ namespace hf_generators
1818
GapTriggeredBeauty, // --> GeneratorPythia8GapTriggeredBeauty: beauty enriched
1919
GapTriggeredCharmAndBeauty, // --> GeneratorPythia8GapTriggeredCharmAndBeauty: charm and beauty enriched (with same ratio)
2020
GapHF, // --> GeneratorPythia8GapHF
21+
GapHFRatio, // --> GeneratorPythia8GapHF, quark list built from b/c ratio
2122
NGenType
2223
};
2324
}
@@ -83,6 +84,12 @@ public:
8384
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
8485
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList, partPdgToReplaceList, freqReplaceList));
8586
break;
87+
88+
case hf_generators::GapHFRatio:
89+
LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GapHFRatio (custom b/c ratio) **********";
90+
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
91+
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList));
92+
break;
8693
default:
8794
LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********";
8895
break;
@@ -387,6 +394,71 @@ private:
387394

388395
};
389396

397+
// Helper: build quarkPdgList from bOverCRatio (= N(beauty)/N(charm), default 1 = 1:1)
398+
// Integer b/c: taken directly; fractional b/c: reduced to nearest nB:nC with nC+nB <= 20
399+
static std::vector<int> BuildQuarkListFromBOverC(float bOverCRatio)
400+
{
401+
const int iterNMax = 19;
402+
if (bOverCRatio <= 0.f) {
403+
LOG(fatal) << "bOverCRatio (b/c) must be > 0";
404+
}
405+
if (bOverCRatio > iterNMax*1.f) {
406+
bOverCRatio = iterNMax*1.f;
407+
LOG(warn) << "bOverCRatio (b/c) too large, using 19:1";
408+
}else if (bOverCRatio < 1.f/iterNMax) {
409+
bOverCRatio = 1.f/iterNMax;
410+
LOG(warn) << "bOverCRatio (b/c) too small, using 1:19";
411+
}
412+
413+
int nC = 1, nB = 1;
414+
float bestErr = 1e9f;
415+
for (int c = 1; c <= iterNMax; ++c) {
416+
int b = static_cast<int>(std::lround(bOverCRatio * c));
417+
if (b < 1 || c + b > 20) {
418+
continue;
419+
}
420+
float err = std::fabs(static_cast<float>(b) / c - bOverCRatio);
421+
if (err < bestErr) {
422+
/// This check here is needed in case bOverCRatio*c is not integer (it can happen with e.g. bOverCRatio=4./9.)
423+
/// In this case, b is its truncation and the desired ratio is not obtained
424+
/// It means that one needs to continue iterating until bOverCRatio*c is integer, i.e. b not truncated
425+
///
426+
/// Possible cases:
427+
/// 1. we want nB = nC*R, with R integer
428+
/// -> we enter here in the first loop iteration
429+
/// 2. we want more charm than beauty by an integer amount, i.e. nB = nC*R with R=1./N, with N integer
430+
/// -> we enter here after N iterations, when c=N
431+
/// 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.)
432+
/// -> In this case, bOverCRatio*c is not integer, namely b is its truncation and the desired ratio is not obtained
433+
/// The code iterates at most until c becomes equal to the denominator of the fraction
434+
/// 4. we want one of the previous cases, but we assign to bOverCRatio a value that is not rational
435+
/// or such as we do not enter here within 19 iterations
436+
/// -> nC and nB, are not touched, therefore we do not have the desired fraction. We'll need to throw a fatal (*)
437+
///
438+
bestErr = err;
439+
nC = c;
440+
nB = b;
441+
442+
/// If we are at this point, we reached already the desired ratio between b and c.
443+
/// Let's break the loop
444+
break;
445+
}
446+
}
447+
448+
// (*) check if we have the desired fraction
449+
bool isRatioUnity = std::fabs(bOverCRatio-1) < 1e-05;
450+
if (!isRatioUnity && nC==1 && nB==1) {
451+
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!";
452+
}
453+
454+
std::vector<int> list;
455+
// Bresenham interleaving
456+
int len = nC + nB;
457+
for (int k = 0; k < len; ++k)
458+
list.push_back((((k + 1) * nC) / len > (k * nC) / len) ? 4 : 5);
459+
return list;
460+
}
461+
390462
// Charm enriched
391463
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 = {})
392464
{
@@ -420,3 +492,18 @@ FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(bool usePtHardBins = false
420492
return myGen;
421493
}
422494

495+
// Charm and beauty enriched with tunable b/c ratio
496+
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 = {})
497+
{
498+
auto myGen = new GeneratorPythia8EmbedHF();
499+
500+
/// build the quark list from the b/c ratio
501+
auto quarkPdgList = BuildQuarkListFromBOverC(bOverCRatio);
502+
503+
/// setup the internal generator for HF events
504+
myGen->setupGeneratorEvHF(hf_generators::GapHFRatio,
505+
usePtHardBins, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax,
506+
quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList);
507+
508+
return myGen;
509+
}

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.ini‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
#NEV_TEST> 10
2+
### The external generator derives from GeneratorPythia8.
3+
[GeneratorExternal]
4+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
5+
funcName=GeneratorPythia8EmbedHFRatio(3.0)
6+
7+
[GeneratorPythia8]
8+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg
9+
includePartonEvent=true

‎MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar3_PbPb.C‎

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
//std::string path{"tf1/sgn_1_Kine.root"};
4+
5+
int checkPdgQuarkOne{4};
6+
int checkPdgQuarkTwo{5};
7+
float ratioTrigger = 1.; // one event triggered out of 1
8+
9+
std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4132, 4232, 4332};
10+
std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
11+
{421, {
12+
{-321, 211}, // D0 -> K-, pi+
13+
{-321, 211, 111}, // D0 -> K-, pi+, pi0
14+
{213, -321}, // D0 -> rho(770)+, K-
15+
{-313, 111}, // D0 -> Kbar^*(892)0, pi0
16+
{-323, 211}, // D0 -> K^*(892)-, pi+
17+
{-211, 211}, // D0 -> pi-, pi+
18+
{213, -211}, // D0 -> rho(770)+, pi-
19+
{-211, 211, 111}, // D0 -> pi-, pi+, pi0
20+
{-321, 321}, // D0 -> K-, K+
21+
}},
22+
23+
{411, {
24+
{-321, 211, 211}, // D+ -> K-, pi+, pi+
25+
{-10311, 211}, // D+ -> Kbar0^*(1430)0, pi+
26+
{-313, 211}, // D+ -> Kbar^*(892)0, pi+
27+
{-321, 211, 211, 111}, // D+ -> K-, pi+, pi+, pi0
28+
{333, 211}, // D+ -> phi(1020)0, pi+
29+
{-313, 321}, // D+ -> Kbar^*(892)0, K+
30+
{-10311, 321}, // D+ -> Kbar0^*(1430)0, K+
31+
{-321, 321, 211}, // D+ -> K-, K+, pi+
32+
{113, 211}, // D+ -> rho(770)0, pi+
33+
{225, 211}, // D+ -> f2(1270)0, pi+
34+
{-211, 211, 211}, // D+ -> pi-, pi+, pi+
35+
}},
36+
37+
{431, {
38+
{333, 211}, // Ds+ -> phi(1020)0, pi+
39+
{-313, 321}, // Ds+ -> Kbar^*(892)0, K+
40+
{333, 213}, // Ds+ -> phi(1020)0, rho(770)+
41+
{113, 211}, // Ds+ -> rho(770)0, pi+
42+
{225, 211}, // Ds+ -> f2(1270)0, pi+
43+
{-211, 211, 211}, // Ds+ -> pi-, pi+, pi+
44+
{313, 211}, // Ds+ -> K^*(892)0, pi+
45+
{10221, 321}, // Ds+ -> f0(1370)0, K+
46+
{113, 321}, // Ds+ -> rho(770)0, K+
47+
{-211, 321, 211}, // Ds+ -> pi-, K+, pi+
48+
{221, 211}, // Ds+ -> eta, pi+
49+
}},
50+
51+
{4122, {
52+
{2212, -321, 211}, // Lambdac+ -> p, K-, pi+
53+
{2212, -313}, // Lambdac+ -> p, Kbar^*(892)0
54+
{2224, -321}, // Lambdac+ -> Delta(1232)++, K-
55+
{102134, 211}, // Lambdac+ -> 102134, pi+
56+
{2212, 310}, // Lambdac+ -> p, K0s
57+
{2212, -321, 211, 111}, // Lambdac+ -> p, K-, pi+, pi0
58+
{2212, -211, 211}, // Lambdac+ -> p, pi-, pi+
59+
{2212, 333}, // Lambdac+ -> p, phi(1020)0
60+
}},
61+
62+
{4232, {
63+
{2212, -321, 211}, // Xic+ -> p, K-, pi+
64+
{2212, -313}, // Xic+ -> p, Kbar^*(892)0
65+
{3312, 211, 211}, // Xic+ -> Xi-, pi+, pi+
66+
{2212, 333}, // Xic+ -> p, phi(1020)0
67+
{3222, -211, 211}, // Xic+ -> Sigma+, pi-, pi+
68+
{3324, 211}, // Xic+ -> Xi(1530)0, pi+
69+
}},
70+
71+
{4132, {
72+
{3312, 211}, // Xic0 -> Xi-, pi+
73+
}},
74+
75+
{4332, {
76+
{3334, 211}, // Omegac0 -> Omega-, pi+
77+
{3312, 211}, // Omegac0 -> Xi-, pi+
78+
}},
79+
};
80+
81+
TFile file(path.c_str(), "READ");
82+
if (file.IsZombie()) {
83+
std::cerr << "Cannot open ROOT file " << path << "\n";
84+
return 1;
85+
}
86+
87+
auto tree = (TTree *)file.Get("o2sim");
88+
std::vector<o2::MCTrack> *tracks{};
89+
tree->SetBranchAddress("MCTrack", &tracks);
90+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
91+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
92+
93+
int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
94+
int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
95+
auto nEvents = tree->GetEntries();
96+
97+
for (int i = 0; i < nEvents; i++) {
98+
tree->GetEntry(i);
99+
100+
// check subgenerator information
101+
//if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
102+
// bool isValid = false;
103+
// int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
104+
// if (subGeneratorId == 0) {
105+
// nEventsMB++;
106+
// } else if (subGeneratorId == checkPdgQuarkOne) {
107+
// nEventsInjOne++;
108+
// } else if (subGeneratorId == checkPdgQuarkTwo) {
109+
// nEventsInjTwo++;
110+
// }
111+
//}
112+
113+
for (auto &track : *tracks) {
114+
auto pdg = track.GetPdgCode();
115+
if (std::abs(pdg) == checkPdgQuarkOne) {
116+
nQuarksOne++;
117+
continue;
118+
}
119+
if (std::abs(pdg) == checkPdgQuarkTwo) {
120+
nQuarksTwo++;
121+
continue;
122+
}
123+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
124+
nSignals++; // count signal PDG
125+
126+
std::vector<int> pdgsDecay{};
127+
std::vector<int> pdgsDecayAntiPart{};
128+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
129+
auto pdgDau = tracks->at(j).GetPdgCode();
130+
pdgsDecay.push_back(pdgDau);
131+
if (pdgDau != 333 && pdgDau != 111 && pdgDau != 221 && pdgDau != 113 && pdgDau != 225) { // phi is antiparticle of itself
132+
pdgsDecayAntiPart.push_back(-pdgDau);
133+
} else {
134+
pdgsDecayAntiPart.push_back(pdgDau);
135+
}
136+
}
137+
138+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
139+
std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
140+
141+
for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
142+
std::sort(decay.begin(), decay.end());
143+
if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
144+
nSignalGoodDecay++;
145+
break;
146+
}
147+
}
148+
}
149+
}
150+
}
151+
152+
std::cout << "--------------------------------\n";
153+
std::cout << "# Events: " << nEvents << "\n";
154+
//std::cout << "# MB events: " << nEventsMB << "\n";
155+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
156+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
157+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
158+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
159+
std::cout <<"# signal hadrons: " << nSignals << "\n";
160+
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
161+
162+
//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
163+
// std::cerr << "Number of generated MB events different than expected\n";
164+
// return 1;
165+
//}
166+
//if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
167+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
168+
// return 1;
169+
//}
170+
//if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
171+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
172+
// return 1;
173+
//}
174+
175+
if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
176+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
177+
return 1;
178+
}
179+
if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
180+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
181+
return 1;
182+
}
183+
184+
float fracForcedDecays = nSignals ? float(nSignalGoodDecay) / nSignals : 0.0f;
185+
float uncFracForcedDecays = nSignals ? std::sqrt(fracForcedDecays * (1 - fracForcedDecays) / nSignals) : 1.0f;
186+
if (1 - fracForcedDecays > 0.15 + uncFracForcedDecays) { // we put some tolerance (e.g. due to oscillations which might change the final state)
187+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
188+
return 1;
189+
}
190+
191+
return 0;
192+
}

‎MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV.cfg‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -107,7 +107,7 @@ BeamRemnants:saturation 5
107107
4122:addChannel = 1 0.04500 100 2224 -321 ### Λc+ -> Delta++ K- 1.08%
108108
4122:addChannel = 1 0.09000 100 102134 211 ### Λc+ -> Lambda(1520) K- 2.20e-3
109109
### Λc+ -> p K0S (36%)
110-
4122:addChannel = 1 0.36000 0 2212 311 ### Λc+ -> p K0S 1.59%
110+
4122:addChannel = 1 0.36000 0 2212 310 ### Λc+ -> p K0S 1.59%
111111
### Λc+ -> p K- π+ π0 (small, 3%)
112112
4122:addChannel = 1 0.03000 0 2212 -321 211 111 ### Λc+ -> p K- π+ π0 (non-resonant) 4.6%
113113
### Λc+ -> p π- π+ (12.50%)
@@ -256,7 +256,7 @@ BeamRemnants:saturation 5
256256
431:onIfMatch = 221 211
257257

258258
### Λc -> pK0s
259-
4122:onIfMatch = 2212 311
259+
4122:onIfMatch = 2212 310
260260
### Λc -> p K- π+ π0
261261
4122:onIfMatch = 2212 321 211
262262
### Λc -> p K*

0 commit comments

Comments
 (0)