Skip to content

Commit 4580f77

Browse files
JinjooSeoJinjooSeo
committed
DQ PbPb MC
1 parent 41467f3 commit 4580f77

File tree

3 files changed

+419
-0
lines changed

3 files changed

+419
-0
lines changed

‎MC/config/PWGDQ/external/generator_pythia8_HadronTriggered_PbPb.C‎

Lines changed: 390 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,390 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/Pythia.h"
4+
#include "TRandom.h"
5+
#include "TParticle.h"
6+
7+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen)
8+
#include "GeneratorEvtGen.C"
9+
10+
#include <string>
11+
12+
using namespace o2::eventgen;
13+
14+
namespace o2
15+
{
16+
namespace eventgen
17+
{
18+
19+
class GeneratorPythia8HadronTriggeredPbPb : public o2::eventgen::GeneratorPythia8 {
20+
public:
21+
22+
/// constructor
23+
GeneratorPythia8HadronTriggeredPbPb(int inputTriggerRatio = 5) {
24+
25+
mGeneratedEvents = 0;
26+
mInverseTriggerRatio = inputTriggerRatio;
27+
// define minimum bias event generator
28+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
29+
TString pathconfigMB = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_PbPb_5TeV.cfg");
30+
pythiaMBgen.readFile(pathconfigMB.Data());
31+
pythiaMBgen.readString("Random:setSeed on");
32+
pythiaMBgen.readString("Random:seed " + std::to_string(seed));
33+
mConfigMBdecays = "";
34+
mRapidityMin = -1.;
35+
mRapidityMax = 1.;
36+
mVerbose = false;
37+
}
38+
39+
/// Destructor
40+
~GeneratorPythia8HadronTriggeredPbPb() = default;
41+
42+
void addHadronPDGs(int pdg) { mHadronsPDGs.push_back(pdg); };
43+
44+
void setRapidityRange(double valMin, double valMax)
45+
{
46+
mRapidityMin = valMin;
47+
mRapidityMax = valMax;
48+
};
49+
50+
void setTriggerGap(int triggerGap) {mInverseTriggerRatio = triggerGap;}
51+
52+
void setConfigMBdecays(TString val){mConfigMBdecays = val;}
53+
54+
void setVerbose(bool val) { mVerbose = val; };
55+
56+
protected:
57+
58+
bool generateEvent() override {
59+
return true;
60+
}
61+
62+
bool Init() override {
63+
if(mConfigMBdecays.Contains("cfg")) {
64+
pythiaMBgen.readFile(mConfigMBdecays.Data());
65+
}
66+
GeneratorPythia8::Init();
67+
pythiaMBgen.init();
68+
return true;
69+
}
70+
71+
std::vector<int> findAllCharmonia(const Pythia8::Event& event) {
72+
std::vector<int> out; out.reserve(4);
73+
74+
for (int ipa = 0; ipa < event.size(); ++ipa) {
75+
76+
auto daughterList = event[ipa].daughterList();
77+
78+
for (auto ida : daughterList) {
79+
for (int pdg : mHadronsPDGs) { // check that at least one of the pdg code is found in the event
80+
if (event[ida].id() == pdg) {
81+
if ((event[ida].y() > mRapidityMin) && (event[ida].y() < mRapidityMax)) {
82+
cout << "============= Found jpsi y,pt " << event[ida].y() << ", " << event[ida].pT() << endl;
83+
out.push_back(ida);
84+
}
85+
}
86+
}
87+
}
88+
}
89+
90+
return out;
91+
}
92+
93+
void collectAncestors(const Pythia8::Event& event, int idx, std::vector<int>& decayChains, std::vector<char>& visited) {
94+
if (idx < 0 || idx >= event.size()) return;
95+
if (!visited[idx]) {
96+
visited[idx] = 1;
97+
decayChains.push_back(idx);
98+
}
99+
100+
const int idabs = std::abs(event[idx].id());
101+
if (idabs == 4 || idabs == 5 || idabs == 21) return;
102+
103+
int mother1 = event[idx].mother1();
104+
int mother2 = event[idx].mother2();
105+
if (mother1 < 0) return;
106+
if (mother2 < mother1) mother2 = mother1;
107+
for (int m = mother1; m <= mother2; ++m) {
108+
if (m == idx) continue;
109+
collectAncestors(event, m, decayChains, visited);
110+
}
111+
}
112+
113+
void collectDaughters(const Pythia8::Event& event, int idx, std::vector<int>& decayChains, std::vector<char>& visited) {
114+
if (idx < 0 || idx >= event.size()) return;
115+
if (visited[idx] == 0) {
116+
decayChains.push_back(idx);
117+
}
118+
119+
if (visited[idx] == 2) return;
120+
visited[idx] = 2;
121+
122+
int daughter1 = event[idx].daughter1();
123+
int daughter2 = event[idx].daughter2();
124+
if (daughter1 < 0) return;
125+
if (daughter2 < daughter1) daughter2 = daughter1;
126+
for (int d = daughter1; d <= daughter2; ++d) {
127+
if (d == idx) continue;
128+
collectDaughters(event, d, decayChains, visited);
129+
}
130+
}
131+
132+
TParticle makeTParticleTemp(const Pythia8::Event& event, int idx) {
133+
const auto& q = event[idx];
134+
int status = q.status();
135+
if (status < 0) {
136+
return TParticle(0, 0, -1, -1, -1, -1,
137+
0.,0.,0.,0., 0.,0.,0.,0.);
138+
}
139+
140+
int m1 = q.mother1();
141+
int m2 = q.mother2();
142+
int d1 = q.daughter1();
143+
int d2 = q.daughter2();
144+
return TParticle(q.id(), status, m1, m2, d1, d2,
145+
q.px(), q.py(), q.pz(), q.e(),
146+
q.xProd(), q.yProd(), q.zProd(), q.tProd());
147+
}
148+
149+
Bool_t importParticles() override
150+
{
151+
//LOG(info) << "";
152+
//LOG(info) << "*************************************************************";
153+
//LOG(info) << "************** New signal event considered **************";
154+
//LOG(info) << "*************************************************************";
155+
//LOG(info) << "";
156+
157+
const int nSig = std::max(1, (int)std::lround(mNumSigEvs));
158+
for (int isig=0; isig<nSig; ++isig) {
159+
160+
bool genOk=false;
161+
std::vector<int> charmonia;
162+
while (! (genOk && !charmonia.empty())) {
163+
/// reset event
164+
mPythia.event.reset();
165+
genOk = GeneratorPythia8::generateEvent();
166+
if (!genOk) continue;
167+
charmonia = findAllCharmonia(mPythia.event);
168+
}
169+
170+
std::vector<int> decayChains;
171+
std::vector<char> visited(mPythia.event.size(), 0);
172+
decayChains.reserve(256);
173+
174+
// find all ancestors of the charmonia
175+
for (size_t ic = 0; ic < charmonia.size(); ++ic) {
176+
int cidx = charmonia[ic];
177+
collectAncestors(mPythia.event, cidx, decayChains, visited);
178+
}
179+
180+
// find all daughters of the charmonia
181+
for (size_t ic = 0; ic < charmonia.size(); ++ic) {
182+
int cidx = charmonia[ic];
183+
collectDaughters(mPythia.event, cidx, decayChains, visited);
184+
}
185+
186+
std::vector<int> idxMap(mPythia.event.size(), -1);
187+
mParticles.reserve(mParticles.size() + (int)decayChains.size());
188+
189+
for (int i = 0; i < (int)decayChains.size(); ++i) {
190+
const int srcIdx = decayChains[i];
191+
if (srcIdx < 0 || srcIdx >= mPythia.event.size()) continue;
192+
193+
TParticle part = makeTParticleTemp(mPythia.event, srcIdx);
194+
if(part.GetPdgCode() == 0) continue;
195+
196+
int newIdx = (int)mParticles.size();
197+
mParticles.push_back(part);
198+
idxMap[srcIdx] = newIdx;
199+
}
200+
201+
for (int iLoc = 0; iLoc < (int) decayChains.size(); ++iLoc) {
202+
const int srcIdx = decayChains[iLoc];
203+
if (srcIdx < 0 || srcIdx >= (int)idxMap.size()) continue;
204+
const int outIdx = idxMap[srcIdx];
205+
if (outIdx < 0) continue;
206+
207+
const auto& src = mPythia.event[srcIdx];
208+
209+
const int mother1 = (src.mother1() >= 0 ? idxMap[src.mother1()] : -1);
210+
const int mother2 = (src.mother2() >= 0 ? idxMap[src.mother2()] : -1);
211+
const int daughter1 = (src.daughter1()>= 0 ? idxMap[src.daughter1()] : -1);
212+
const int daughter2 = (src.daughter2()>= 0 ? idxMap[src.daughter2()] : -1);
213+
214+
// update TParticle
215+
TParticle& particle = mParticles[outIdx];
216+
particle.SetFirstMother(mother1);
217+
particle.SetLastMother(mother2);
218+
particle.SetFirstDaughter(daughter1);
219+
particle.SetLastDaughter(daughter2);
220+
}
221+
LOG(info) << "-----------------------------------------------";
222+
LOG(info) << "============ After event " << isig << " (size " << decayChains.size() << ")";
223+
LOG(info) << "Full stack (size " << mParticles.size() << "):";
224+
LOG(info) << "-----------------------------------------------";
225+
// printParticleVector(mParticles);
226+
}
227+
228+
if (mVerbose) mOutputEvent.list();
229+
230+
return kTRUE;
231+
}
232+
233+
void notifyEmbedding(const o2::dataformats::MCEventHeader* bkgHeader) override {
234+
LOG(info) << "[notifyEmbedding] ----- Function called";
235+
236+
/// Impact parameter between the two nuclei
237+
const float x = bkgHeader->GetB();
238+
LOG(info) << "[notifyEmbedding] ----- Collision impact parameter: " << x;
239+
240+
/// number of events to be embedded in a background event
241+
mNumSigEvs = 5 + 0.886202881*std::pow(std::max(0.0f, 17.5f - x),1.7);
242+
LOG(info) << "[notifyEmbedding] ----- generating " << mNumSigEvs << " signal events " << std::endl;
243+
};
244+
245+
private:
246+
// Interface to override import particles
247+
Pythia8::Event mOutputEvent;
248+
249+
// Control gap-triggering
250+
unsigned long long mGeneratedEvents;
251+
int mInverseTriggerRatio;
252+
Pythia8::Pythia pythiaMBgen; // minimum bias event
253+
TString mConfigMBdecays;
254+
std::vector<int> mHadronsPDGs;
255+
double mRapidityMin;
256+
double mRapidityMax;
257+
bool mVerbose;
258+
259+
// number of signal events to be embedded in a background event
260+
int mNumSigEvs{1};
261+
};
262+
263+
}
264+
265+
}
266+
267+
// Predefined generators:
268+
FairGenerator*
269+
GeneratorPromptJpsi_EvtGenMidY(int triggerGap, double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, bool embedding = false)
270+
{
271+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8HadronTriggeredPbPb>();
272+
gen->setTriggerGap(triggerGap);
273+
gen->setRapidityRange(rapidityMin, rapidityMax);
274+
gen->addHadronPDGs(443);
275+
gen->setVerbose(verbose);
276+
277+
TString pathO2table = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffJpsi.cfg");
278+
gen->readFile(pathO2table.Data());
279+
gen->setConfigMBdecays(pathO2table);
280+
gen->PrintDebug(true);
281+
282+
gen->SetSizePdg(1);
283+
gen->AddPdg(443, 0);
284+
285+
gen->SetForceDecay(kEvtDiElectron);
286+
287+
// set random seed
288+
gen->readString("Random:setSeed on");
289+
uint random_seed;
290+
unsigned long long int random_value = 0;
291+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
292+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
293+
gen->readString(Form("Random:seed = %llu", random_value % 900000001));
294+
295+
// print debug
296+
// gen->PrintDebug();
297+
298+
return gen;
299+
}
300+
301+
FairGenerator*
302+
GeneratorPromptJpsiPsi2S_EvtGenMidY(int triggerGap, double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, bool embedding = false)
303+
{
304+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8HadronTriggeredPbPb>();
305+
gen->setTriggerGap(triggerGap);
306+
gen->setRapidityRange(rapidityMin, rapidityMax);
307+
gen->addHadronPDGs(443);
308+
gen->addHadronPDGs(100443);
309+
gen->setVerbose(verbose);
310+
TString pathO2table = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffJpsi.cfg");
311+
gen->readFile(pathO2table.Data());
312+
gen->setConfigMBdecays(pathO2table);
313+
gen->PrintDebug(true);
314+
gen->SetSizePdg(2);
315+
gen->AddPdg(443, 0);
316+
gen->AddPdg(100443, 1);
317+
gen->SetForceDecay(kEvtDiElectron);
318+
// set random seed
319+
gen->readString("Random:setSeed on");
320+
uint random_seed;
321+
unsigned long long int random_value = 0;
322+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
323+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
324+
gen->readString(Form("Random:seed = %llu", random_value % 900000001));
325+
// print debug
326+
// gen->PrintDebug();
327+
return gen;
328+
}
329+
330+
FairGenerator*
331+
GeneratorPromptJpsi_EvtGenFwdy(int triggerGap, double rapidityMin = -4.3, double rapidityMax = -2.3, bool verbose = false, bool embedding = false)
332+
{
333+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8HadronTriggeredPbPb>();
334+
gen->setTriggerGap(triggerGap);
335+
gen->setRapidityRange(rapidityMin, rapidityMax);
336+
gen->addHadronPDGs(443);
337+
gen->setVerbose(verbose);
338+
339+
TString pathO2table = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffJpsi.cfg");
340+
gen->readFile(pathO2table.Data());
341+
gen->setConfigMBdecays(pathO2table);
342+
gen->PrintDebug(true);
343+
344+
gen->SetSizePdg(1);
345+
gen->AddPdg(443, 0);
346+
347+
gen->SetForceDecay(kEvtDiMuon);
348+
349+
// set random seed
350+
gen->readString("Random:setSeed on");
351+
uint random_seed;
352+
unsigned long long int random_value = 0;
353+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
354+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
355+
gen->readString(Form("Random:seed = %llu", random_value % 900000001));
356+
357+
// print debug
358+
// gen->PrintDebug();
359+
360+
return gen;
361+
}
362+
363+
FairGenerator*
364+
GeneratorPromptJpsiPsi2S_EvtGenFwdY(int triggerGap, double rapidityMin = -4.3, double rapidityMax = -2.3, bool verbose = false, bool embedding = false)
365+
{
366+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8HadronTriggeredPbPb>();
367+
gen->setTriggerGap(triggerGap);
368+
gen->setRapidityRange(rapidityMin, rapidityMax);
369+
gen->addHadronPDGs(443);
370+
gen->addHadronPDGs(100443);
371+
gen->setVerbose(verbose);
372+
TString pathO2table = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffJpsi.cfg");
373+
gen->readFile(pathO2table.Data());
374+
gen->setConfigMBdecays(pathO2table);
375+
gen->PrintDebug(true);
376+
gen->SetSizePdg(2);
377+
gen->AddPdg(443, 0);
378+
gen->AddPdg(100443, 1);
379+
gen->SetForceDecay(kEvtDiMuon);
380+
// set random seed
381+
gen->readString("Random:setSeed on");
382+
uint random_seed;
383+
unsigned long long int random_value = 0;
384+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
385+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
386+
gen->readString(Form("Random:seed = %llu", random_value % 900000001));
387+
// print debug
388+
// gen->PrintDebug();
389+
return gen;
390+
}

‎MC/config/PWGDQ/ini/Generator_InjectedCharmoniaMidy_PbPb5TeV.ini‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
### The external generator derives from GeneratorPythia8.
2+
[GeneratorExternal]
3+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_HadronTriggered_PbPb.C
4+
funcName=GeneratorPromptJpsiPsi2S_EvtGenMidY(1,-1,1,false)
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_PbPb_5TeV.cfg

0 commit comments

Comments
 (0)