Skip to content

Commit 76d94ba

Browse files
mfagginMattia Faggin
authored and
Benedikt Volkel
committed
PWGHF: Pb-Pb generator to embed N HF events in a single underlying event. (#1666)
* Pb-Pb generator to embed N HF events in a single underlying event. * Fix HF ev. generator creation. * fix * Fix particle stack filling and mother/daughter indixing. * Add .ini configuration file and tester, and hardQCD.cfg file. * Define mNumSigEvs depending on the collision impact parameter. * Fix o2sim_Kine path in test. * Add description. * Update author list of .cfg file. * Implement suggestions by fgrosa. * Function documentation. * Overriding generate event, to avoid simulating 1 untriggered event. --------- Co-authored-by: Mattia Faggin <mfaggin@alipap1.cern.ch> (cherry picked from commit 36522d4)
1 parent 90f54a1 commit 76d94ba

File tree

4 files changed

+476
-0
lines changed

4 files changed

+476
-0
lines changed

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

Lines changed: 191 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,191 @@
1+
///////////////////////////////////////////////////////////////////////////////
2+
/// ///
3+
/// HF MC generator for Pb-Pb ///
4+
/// Option 1: generate N PYTHIA events triggered on ccbar and/or bbbar ///
5+
/// to be embedded with a underlying Pb-Pb event ///
6+
/// ///
7+
///////////////////////////////////////////////////////////////////////////////
8+
9+
#include "generator_pythia8_gaptriggered_hf.C"
10+
11+
using namespace Pythia8;
12+
13+
namespace hf_generators
14+
{
15+
enum GenType : int {
16+
GapTriggeredCharm = 0, // --> GeneratorPythia8GapTriggeredCharm: charm enriched
17+
GapTriggeredBeauty, // --> GeneratorPythia8GapTriggeredBeauty: beauty enriched
18+
GapTriggeredCharmAndBeauty, // --> GeneratorPythia8GapTriggeredCharmAndBeauty: charm and beauty enriched (with same ratio)
19+
GapHF, // --> GeneratorPythia8GapHF
20+
NGenType
21+
};
22+
}
23+
24+
class GeneratorPythia8EmbedHF : public o2::eventgen::GeneratorPythia8
25+
{
26+
public:
27+
28+
/// default constructor
29+
GeneratorPythia8EmbedHF() = default;
30+
31+
/// constructor
32+
//GeneratorPythia8EmbedHF(int numSigEvs = 1) {
33+
// mNumSigEvs = numSigEvs;
34+
// //mGeneratedEvents = 0;
35+
//}
36+
37+
/// Destructor
38+
~GeneratorPythia8EmbedHF() = default;
39+
40+
/// Init
41+
bool Init() override
42+
{
43+
return o2::eventgen::GeneratorPythia8::Init();
44+
}
45+
46+
/// @brief setup the event generator for HF signals
47+
/// \param gentype generator type (only ccbar, only bbbar, both)
48+
/// \param yQuarkMin minimum quark rapidity
49+
/// \param yQuarkMax maximum quark rapidity
50+
/// \param yHadronMin minimum hadron rapidity
51+
/// \param yHadronMax maximum hadron rapidity
52+
/// \param hadronPdgList list of PDG codes for hadrons to be used in trigger
53+
void setupGeneratorEvHF(int genType, float yQuarkMin, float yQuarkMax, float yHadronMin, float yHadronMax, std::vector<int> hadronPdgList = {}) {
54+
mGeneratorEvHF = nullptr;
55+
switch (genType)
56+
{
57+
case hf_generators::GapTriggeredCharm:
58+
LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredCharm **********";
59+
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
60+
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharm(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
61+
break;
62+
case hf_generators::GapTriggeredBeauty:
63+
LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredBeauty **********";
64+
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
65+
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
66+
break;
67+
case hf_generators::GapTriggeredCharmAndBeauty:
68+
LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapTriggeredCharmAndBeauty **********";
69+
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
70+
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapTriggeredCharmAndBeauty(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
71+
break;
72+
case hf_generators::GapHF:
73+
LOG(info) << "********** [GeneratorPythia8EmbedHF] configuring GeneratorPythia8GapHF **********";
74+
LOG(info) << "********** Default number of HF signal events to be merged (updated by notifyEmbedding): " << mNumSigEvs;
75+
mGeneratorEvHF = dynamic_cast<GeneratorPythia8GapTriggeredHF*>(GeneratorPythia8GapHF(/*no gap trigger*/1, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList));
76+
break;
77+
default:
78+
LOG(fatal) << "********** [GeneratorPythia8EmbedHF] bad configuration, fix it! **********";
79+
break;
80+
}
81+
mGeneratorEvHF->Init();
82+
}
83+
84+
// This function is called by the primary generator
85+
// for each event in case we are in embedding mode.
86+
// We use it to setup the number of signal events
87+
// to be generated and to be embedded on the background.
88+
void notifyEmbedding(const o2::dataformats::MCEventHeader* bkgHeader) override
89+
{
90+
LOG(info) << "[notifyEmbedding] ----- Function called";
91+
92+
/// Impact parameter between the two nuclei
93+
const float x = bkgHeader->GetB();
94+
LOG(info) << "[notifyEmbedding] ----- Collision impact parameter: " << x;
95+
96+
/// number of events to be embedded in a background event
97+
mNumSigEvs = std::max(1.,120.*(x<5.)+80.*(1.-x/20.)*(x>5.)*(x<11.)+240.*(1.-x/13.)*(x>11.));
98+
LOG(info) << "[notifyEmbedding] ----- generating " << mNumSigEvs << " signal events " << std::endl;
99+
};
100+
101+
protected:
102+
103+
/// @brief Main function for event generation
104+
bool generateEvent() override
105+
{
106+
/// Overriding that from GeneratorPythia8, to avoid the simulation of an untriggered event as first
107+
return true;
108+
}
109+
110+
/// @brief Main function to copy the generated particles in mPythia.event into the stack (this.mParticles)
111+
Bool_t importParticles() override
112+
{
113+
/// Import particles from generated event
114+
/// This should not do anything now, since we override generateEvent
115+
GeneratorPythia8::importParticles();
116+
117+
/// Generate mNumSigEvs HF events to be merged in one
118+
int nEvsHF = 0;
119+
while(nEvsHF < mNumSigEvs) {
120+
121+
/// generate the HF event
122+
bool genOk = false;
123+
while(!genOk) {
124+
genOk = (mGeneratorEvHF->generateEvent() && mGeneratorEvHF->importParticles() /*copy particles from mGeneratorEvHF.mPythia.event to mGeneratorEvHF.mParticles*/ );
125+
}
126+
127+
/// copy the particles from the HF event in the particle stack
128+
auto particlesHfEvent = mGeneratorEvHF->getParticles();
129+
int originalSize = mParticles.size(); // stack of this event generator
130+
for(int iPart=0; iPart<particlesHfEvent.size(); iPart++) {
131+
auto particle = particlesHfEvent.at(iPart);
132+
133+
/// adjust the particle mother and daughter indices
134+
if(particle.GetFirstMother() >= 0) particle.SetFirstMother(particle.GetFirstMother() + originalSize);
135+
if(particle.GetFirstDaughter() >= 0) particle.SetFirstDaughter(particle.GetFirstDaughter() + originalSize);
136+
if(particle.GetLastDaughter() >= 0) particle.SetLastDaughter(particle.GetLastDaughter() + originalSize);
137+
138+
/// copy inside this.mParticles from mGeneratorEvHF.mParticles, i.e. the particles generated in mGeneratorEvHF
139+
mParticles.push_back(particle);
140+
}
141+
142+
/// one more event generated, let's update the counter and clear it, to allow the next generation
143+
nEvsHF++;
144+
//mGeneratedEvents++;
145+
mGeneratorEvHF->clearParticles();
146+
}
147+
148+
return true;
149+
}
150+
151+
private:
152+
153+
Generator* mGeneratorEvHF; // to generate HF signal events
154+
155+
int mNumSigEvs{1}; // number of HF signal events to be merged in one Pythia event
156+
//unsigned long long mGeneratedEvents;
157+
158+
};
159+
160+
// Charm enriched
161+
FairGenerator * GeneratorPythia8EmbedHFCharm(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
162+
{
163+
auto myGen = new GeneratorPythia8EmbedHF();
164+
165+
/// setup the internal generator for HF events
166+
myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharm, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
167+
168+
return myGen;
169+
}
170+
171+
// Beauty enriched
172+
FairGenerator * GeneratorPythia8EmbedHFBeauty(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
173+
{
174+
auto myGen = new GeneratorPythia8EmbedHF();
175+
176+
/// setup the internal generator for HF events
177+
myGen->setupGeneratorEvHF(hf_generators::GapTriggeredBeauty, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
178+
179+
return myGen;
180+
}
181+
182+
// Charm and beauty enriched (with same ratio)
183+
FairGenerator * GeneratorPythia8EmbedHFCharmAndBeauty(float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
184+
{
185+
auto myGen = new GeneratorPythia8EmbedHF();
186+
187+
/// setup the internal generator for HF events
188+
myGen->setupGeneratorEvHF(hf_generators::GapTriggeredCharmAndBeauty, yQuarkMin, yQuarkMax, yHadronMin, yHadronMax, hadronPdgList);
189+
190+
return myGen;
191+
}

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb.ini‎

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

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

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
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+
{411, {{-321, 211, 211}, {-313, 211}, {211, 311}, {211, 333}}}, // D+
12+
{421, {{-321, 211}, {-321, 111, 211}}}, // D0
13+
{431, {{211, 333}, {-313, 321}}}, // Ds+
14+
{4122, {{-313, 2212}, {-321, 2224}, {211, 3124}, {-321, 211, 2212}, {311, 2212}}}, // Lc+
15+
{4132, {{211, 3312}}}, // Xic0
16+
{4232, {{-313, 2212}, {-321, 3324}, {211, 211, 3312}, {-321, 211, 2212}}}, // Xic+
17+
{4332, {{211, 3334}}} // Omegac+
18+
};
19+
20+
TFile file(path.c_str(), "READ");
21+
if (file.IsZombie()) {
22+
std::cerr << "Cannot open ROOT file " << path << "\n";
23+
return 1;
24+
}
25+
26+
auto tree = (TTree *)file.Get("o2sim");
27+
std::vector<o2::MCTrack> *tracks{};
28+
tree->SetBranchAddress("MCTrack", &tracks);
29+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
30+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
31+
32+
int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
33+
int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
34+
auto nEvents = tree->GetEntries();
35+
36+
for (int i = 0; i < nEvents; i++) {
37+
tree->GetEntry(i);
38+
39+
// check subgenerator information
40+
//if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
41+
// bool isValid = false;
42+
// int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
43+
// if (subGeneratorId == 0) {
44+
// nEventsMB++;
45+
// } else if (subGeneratorId == checkPdgQuarkOne) {
46+
// nEventsInjOne++;
47+
// } else if (subGeneratorId == checkPdgQuarkTwo) {
48+
// nEventsInjTwo++;
49+
// }
50+
//}
51+
52+
for (auto &track : *tracks) {
53+
auto pdg = track.GetPdgCode();
54+
if (std::abs(pdg) == checkPdgQuarkOne) {
55+
nQuarksOne++;
56+
continue;
57+
}
58+
if (std::abs(pdg) == checkPdgQuarkTwo) {
59+
nQuarksTwo++;
60+
continue;
61+
}
62+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
63+
nSignals++; // count signal PDG
64+
65+
std::vector<int> pdgsDecay{};
66+
std::vector<int> pdgsDecayAntiPart{};
67+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
68+
auto pdgDau = tracks->at(j).GetPdgCode();
69+
pdgsDecay.push_back(pdgDau);
70+
if (pdgDau != 333) { // phi is antiparticle of itself
71+
pdgsDecayAntiPart.push_back(-pdgDau);
72+
} else {
73+
pdgsDecayAntiPart.push_back(pdgDau);
74+
}
75+
}
76+
77+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
78+
std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
79+
80+
for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
81+
if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
82+
nSignalGoodDecay++;
83+
break;
84+
}
85+
}
86+
}
87+
}
88+
}
89+
90+
std::cout << "--------------------------------\n";
91+
std::cout << "# Events: " << nEvents << "\n";
92+
//std::cout << "# MB events: " << nEventsMB << "\n";
93+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
94+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
95+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
96+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
97+
std::cout <<"# signal hadrons: " << nSignals << "\n";
98+
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
99+
100+
//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
101+
// std::cerr << "Number of generated MB events different than expected\n";
102+
// return 1;
103+
//}
104+
//if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
105+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
106+
// return 1;
107+
//}
108+
//if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
109+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
110+
// return 1;
111+
//}
112+
113+
if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
114+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
115+
return 1;
116+
}
117+
if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
118+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
119+
return 1;
120+
}
121+
122+
float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
123+
if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
124+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
125+
return 1;
126+
}
127+
128+
return 0;
129+
}

0 commit comments

Comments
 (0)