Skip to content

Commit 418fe24

Browse files
rbailhacrbailhac
authored
Generator for HF enhanced PbPb MC for dielectron analysis (#2304)
* HF pp embedded in Pb--Pb simulated events * HF pp events embedded in Pb--Pb events for dielectron analysis * Fix the test of the generators
1 parent b5e99d7 commit 418fe24

7 files changed

+894
-0
lines changed

‎MC/config/PWGEM/external/generator/Generator_pythia8_HFLepton_pp.C‎

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
1+
#include "Pythia8/Pythia.h"
2+
#include "Pythia8/HeavyIons.h"
3+
#include "FairGenerator.h"
4+
#include "FairPrimaryGenerator.h"
5+
#include "Generators/GeneratorPythia8.h"
6+
#include "TRandom3.h"
7+
#include "TParticlePDG.h"
8+
#include "TDatabasePDG.h"
9+
10+
#include <map>
11+
#include <unordered_set>
12+
13+
using namespace Pythia8;
14+
15+
class GeneratorPythia8HFLeptonpp : public o2::eventgen::GeneratorPythia8
16+
{
17+
public:
18+
/// default constructor
19+
GeneratorPythia8HFLeptonpp() = default;
20+
21+
/// constructor
22+
GeneratorPythia8HFLeptonpp(TString configsignal, int quarkPdg = 4, int lInputExternalID = 0)
23+
{
24+
25+
lGeneratedEvents = 0;
26+
lExternalID = lInputExternalID;
27+
mQuarkPdg = quarkPdg;
28+
29+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
30+
31+
int offset = (int)(gRandom->Uniform(1)); // create offset to mitigate edge effects due to small number of events per job
32+
lGeneratedEvents += offset;
33+
34+
cout << "Initalizing PYTHIA object used to generate signal events..." << endl;
35+
TString pathconfigSignal = gSystem->ExpandPathName(configsignal.Data());
36+
pythiaObjectSignal.readFile(pathconfigSignal.Data());
37+
pythiaObjectSignal.readString("Random:setSeed on");
38+
pythiaObjectSignal.readString("Random:seed " + std::to_string(seed));
39+
pythiaObjectSignal.readString("Beams:eCM = 5360.0");
40+
pythiaObjectSignal.init();
41+
cout << "Initalization of signal event is complete" << endl;
42+
43+
// flag the generators using type
44+
// addCocktailConstituent(type, "interesting");
45+
// addCocktailConstitent(0, "minbias");
46+
// Add Sub generators
47+
addSubGenerator(1, "charm lepton");
48+
addSubGenerator(2, "beauty forced decay");
49+
addSubGenerator(3, "beauty no foced decay");
50+
}
51+
52+
/// Destructor
53+
~GeneratorPythia8HFLeptonpp() = default;
54+
55+
void addTriggerOnDaughter(int nb, int pdg)
56+
{
57+
mNbDaughter = nb;
58+
mPdgDaughter = pdg;
59+
};
60+
void setQuarkRapidity(float yMin, float yMax)
61+
{
62+
mQuarkRapidityMin = yMin;
63+
mQuarkRapidityMax = yMax;
64+
};
65+
void setDaughterRapidity(float yMin, float yMax)
66+
{
67+
mDaughterRapidityMin = yMin;
68+
mDaughterRapidityMax = yMax;
69+
};
70+
71+
protected:
72+
//__________________________________________________________________
73+
Bool_t generateEvent() override
74+
{
75+
/// reset event
76+
mPythia.event.reset();
77+
78+
// Generate event of interest
79+
Bool_t lGenerationOK = kFALSE;
80+
while (!lGenerationOK) {
81+
if (pythiaObjectSignal.next()) {
82+
lGenerationOK = selectEvent(pythiaObjectSignal.event);
83+
}
84+
}
85+
mPythia.event = pythiaObjectSignal.event;
86+
notifySubGenerator(lExternalID);
87+
88+
lGeneratedEvents++;
89+
// mPythia.next();
90+
91+
return true;
92+
}
93+
94+
bool selectEvent(const Pythia8::Event& event)
95+
{
96+
bool isGoodAtPartonLevel = false, isGoodAtDaughterLevel = (mPdgDaughter != 0) ? false : true;
97+
int nbDaughter = 0;
98+
for (auto iPart{0}; iPart < event.size(); ++iPart) {
99+
// search for Q-Qbar mother with at least one Q in rapidity window
100+
if (!isGoodAtPartonLevel) {
101+
auto daughterList = event[iPart].daughterList();
102+
bool hasQ = false, hasQbar = false, atSelectedY = false;
103+
for (auto iDau : daughterList) {
104+
if (event[iDau].id() == mQuarkPdg) {
105+
hasQ = true;
106+
}
107+
if (event[iDau].id() == -mQuarkPdg) {
108+
hasQbar = true;
109+
}
110+
if ((std::abs(event[iDau].id()) == mQuarkPdg) && (event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax))
111+
atSelectedY = true;
112+
}
113+
if (hasQ && hasQbar && atSelectedY) {
114+
isGoodAtPartonLevel = true;
115+
}
116+
}
117+
// search for mNbDaughter daughters of type mPdgDaughter in rapidity window
118+
if (!isGoodAtDaughterLevel) {
119+
int id = std::abs(event[iPart].id());
120+
float rap = event[iPart].y();
121+
if (id == mPdgDaughter) {
122+
int motherindexa = event[iPart].mother1();
123+
if (motherindexa > 0) {
124+
int idmother = std::abs(event[motherindexa].id());
125+
if (int(std::abs(idmother) / 100.) == 4 || int(std::abs(idmother) / 1000.) == 4 || int(std::abs(idmother) / 100.) == 5 || int(std::abs(idmother) / 1000.) == 5) {
126+
if (rap > mDaughterRapidityMin && rap < mDaughterRapidityMax) {
127+
nbDaughter++;
128+
if (nbDaughter >= mNbDaughter) isGoodAtDaughterLevel = true;
129+
}
130+
}
131+
}
132+
}
133+
}
134+
// we send the trigger
135+
if (isGoodAtPartonLevel && isGoodAtDaughterLevel) {
136+
return true;
137+
}
138+
}
139+
return false;
140+
};
141+
142+
private:
143+
// Interface to override import particles
144+
Pythia8::Event mOutputEvent;
145+
146+
// Properties of selection
147+
int mQuarkPdg;
148+
float mQuarkRapidityMin;
149+
float mQuarkRapidityMax;
150+
int mPdgDaughter;
151+
int mNbDaughter;
152+
float mDaughterRapidityMin;
153+
float mDaughterRapidityMax;
154+
155+
Long64_t lGeneratedEvents;
156+
// ID for different generators
157+
int lExternalID;
158+
159+
// Base event generators
160+
Pythia8::Pythia pythiaObjectSignal; ///Signal collision generator
161+
};
162+
163+
// Predefined generators:
164+
165+
// Charm-enriched forced decay
166+
FairGenerator* GeneratorPythia8CharmLepton(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
167+
{
168+
auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_pp_cr2_forceddecayscharm.cfg", 4, inputExternalID);
169+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
170+
myGen->readString("Random:setSeed on");
171+
myGen->readString("Random:seed " + std::to_string(seed));
172+
myGen->setQuarkRapidity(yMinQ, yMaxQ);
173+
myGen->addTriggerOnDaughter(2, pdgLepton);
174+
myGen->setDaughterRapidity(yMinL, yMaxL);
175+
return myGen;
176+
}
177+
178+
// Beauty-enriched forced decay
179+
FairGenerator* GeneratorPythia8BeautyForcedDecays(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
180+
{
181+
auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar_forceddecayscharmbeauty.cfg", 5, inputExternalID);
182+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
183+
myGen->readString("Random:setSeed on");
184+
myGen->readString("Random:seed " + std::to_string(seed));
185+
myGen->setQuarkRapidity(yMinQ, yMaxQ);
186+
myGen->addTriggerOnDaughter(2, pdgLepton);
187+
myGen->setDaughterRapidity(yMinL, yMaxL);
188+
return myGen;
189+
}
190+
191+
// Beauty-enriched no forced decay
192+
FairGenerator* GeneratorPythia8BeautyNoForcedDecays(int inputExternalID, int pdgLepton, float yMinQ = -1.5, float yMaxQ = 1.5, float yMinL = -1, float yMaxL = 1)
193+
{
194+
auto myGen = new GeneratorPythia8HFLeptonpp("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar.cfg", 5, inputExternalID);
195+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
196+
myGen->readString("Random:setSeed on");
197+
myGen->readString("Random:seed " + std::to_string(seed));
198+
myGen->setQuarkRapidity(yMinQ, yMaxQ);
199+
myGen->addTriggerOnDaughter(2, pdgLepton);
200+
myGen->setDaughterRapidity(yMinL, yMaxL);
201+
return myGen;
202+
}

0 commit comments

Comments
 (0)