Skip to content

Commit 9813fd4

Browse files
ffiondaBenedikt Volkel
authored andcommitted
update prompt / non-prompt charmonia simulation for introducing trigger gap (#1284)
1 parent 5f65a32 commit 9813fd4

13 files changed

Lines changed: 1229 additions & 0 deletions

MC/config/PWGDQ/external/generator/GeneratorPromptCharmonia.C

Lines changed: 460 additions & 0 deletions
Large diffs are not rendered by default.
Lines changed: 276 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,276 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/Pythia.h"
4+
#include "TRandom.h"
5+
6+
R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen)
7+
#include "GeneratorEvtGen.C"
8+
9+
#include <string>
10+
11+
using namespace o2::eventgen;
12+
using namespace Pythia8;
13+
14+
namespace o2
15+
{
16+
namespace eventgen
17+
{
18+
19+
class GeneratorPythia8NonPromptInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 {
20+
public:
21+
22+
/// constructor
23+
GeneratorPythia8NonPromptInjectedGapTriggeredDQ(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_ROOT/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg");
30+
pythiaMBgen.readFile(pathconfigMB.Data());
31+
pythiaMBgen.readString("Random:setSeed on");
32+
pythiaMBgen.readString("Random:seed " + std::to_string(seed));
33+
mConfigMBdecays = "";
34+
mPDG = 5;
35+
mRapidityMin = -1.;
36+
mRapidityMax = 1.;
37+
mHadronMultiplicity = -1;
38+
mHadronRapidityMin = -1.;
39+
mHadronRapidityMax = 1.;
40+
mVerbose = false;
41+
}
42+
43+
/// Destructor
44+
~GeneratorPythia8NonPromptInjectedGapTriggeredDQ() = default;
45+
46+
void setPDG(int val) { mPDG = val; };
47+
void addHadronPDGs(int pdg) { mHadronsPDGs.push_back(pdg); };
48+
void setHadronMultiplicity(int val) { mHadronMultiplicity = val; };
49+
void setRapidity(double valMin, double valMax)
50+
{
51+
mRapidityMin = valMin;
52+
mRapidityMax = valMax;
53+
};
54+
55+
void setRapidityHadron(double valMin, double valMax)
56+
{
57+
mHadronRapidityMin = valMin;
58+
mHadronRapidityMax = valMax;
59+
};
60+
61+
void setConfigMBdecays(TString val){mConfigMBdecays = val;}
62+
63+
void setVerbose(bool val) { mVerbose = val; };
64+
65+
protected:
66+
67+
Bool_t generateEvent() override
68+
{
69+
// reset event
70+
bool genOk = false;
71+
if (mGeneratedEvents % mInverseTriggerRatio == 0){
72+
bool ancestor = false;
73+
while (! (genOk && ancestor) ){
74+
/// reset event
75+
mPythia.event.reset();
76+
genOk = GeneratorPythia8::generateEvent();
77+
// find the q-qbar ancestor
78+
ancestor = findHeavyQuarkPair(mPythia.event);
79+
}
80+
} else {
81+
/// reset event
82+
pythiaMBgen.event.reset();
83+
while (!genOk) {
84+
genOk = pythiaMBgen.next();
85+
}
86+
mPythia.event = pythiaMBgen.event;
87+
}
88+
mGeneratedEvents++;
89+
if (mVerbose) mOutputEvent.list();
90+
return true;
91+
}
92+
93+
Bool_t Init() override
94+
{
95+
if(mConfigMBdecays.Contains("cfg")) pythiaMBgen.readFile(mConfigMBdecays.Data());
96+
GeneratorPythia8::Init();
97+
pythiaMBgen.init();
98+
return true;
99+
}
100+
101+
// search for q-qbar mother with at least one q in a selected rapidity window
102+
bool findHeavyQuarkPair(Pythia8::Event& event)
103+
{
104+
int countH[mHadronsPDGs.size()]; for(int ipdg=0; ipdg < mHadronsPDGs.size(); ipdg++) countH[ipdg]=0;
105+
bool hasq = false, hasqbar = false, atSelectedY = false, isOkAtPartonicLevel = false;
106+
for (int ipa = 0; ipa < event.size(); ++ipa) {
107+
108+
if(!isOkAtPartonicLevel){
109+
auto daughterList = event[ipa].daughterList();
110+
hasq = false; hasqbar = false; atSelectedY = false;
111+
for (auto ida : daughterList) {
112+
if (event[ida].id() == mPDG)
113+
hasq = true;
114+
if (event[ida].id() == -mPDG)
115+
hasqbar = true;
116+
if ((event[ida].y() > mRapidityMin) && (event[ida].y() < mRapidityMax))
117+
atSelectedY = true;
118+
}
119+
if (hasq && hasqbar && atSelectedY) isOkAtPartonicLevel = true;
120+
}
121+
122+
if( (mHadronMultiplicity <= 0) && isOkAtPartonicLevel) return true; // no selection at hadron level
123+
124+
/// check at hadron level if needed
125+
int ipdg=0;
126+
for (auto& pdgVal : mHadronsPDGs){
127+
if ( (TMath::Abs(event[ipa].id()) == pdgVal) && (event[ipa].y() > mHadronRapidityMin) && (event[ipa].y() < mHadronRapidityMax) ) { countH[ipdg]++; printf("ndaughters %d \n",event[ipa].daughterList().size()); }
128+
if(isOkAtPartonicLevel && countH[ipdg] >= mHadronMultiplicity) return true;
129+
ipdg++;
130+
}
131+
}
132+
return false;
133+
};
134+
135+
136+
private:
137+
// Interface to override import particles
138+
Pythia8::Event mOutputEvent;
139+
140+
// Control gap-triggering
141+
unsigned long long mGeneratedEvents;
142+
int mInverseTriggerRatio;
143+
Pythia pythiaMBgen; // minimum bias event
144+
TString mConfigMBdecays;
145+
int mPDG;
146+
std::vector<int> mHadronsPDGs;
147+
int mHadronMultiplicity;
148+
double mRapidityMin;
149+
double mRapidityMax;
150+
double mHadronRapidityMin;
151+
double mHadronRapidityMax;
152+
bool mVerbose;
153+
};
154+
155+
}
156+
157+
}
158+
159+
// Predefined generators:
160+
FairGenerator*
161+
GeneratorBeautyToJpsi_EvtGenMidY(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "511;521;531;541;5112;5122;5232;5132;5332")
162+
{
163+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8NonPromptInjectedGapTriggeredDQ>();
164+
gen->setRapidity(rapidityMin, rapidityMax);
165+
gen->setPDG(5);
166+
gen->setRapidityHadron(-1.5,1.5);
167+
gen->setHadronMultiplicity(1);
168+
TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBhadrons.cfg");
169+
gen->readFile(pathO2table.Data());
170+
gen->setConfigMBdecays(pathO2table);
171+
gen->setVerbose(verbose);
172+
173+
std::string spdg;
174+
TObjArray* obj = pdgs.Tokenize(";");
175+
gen->SetSizePdg(obj->GetEntriesFast());
176+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
177+
spdg = obj->At(i)->GetName();
178+
gen->AddPdg(std::stoi(spdg), i);
179+
gen->addHadronPDGs(std::stoi(spdg));
180+
printf("PDG %d \n", std::stoi(spdg));
181+
}
182+
gen->SetForceDecay(kEvtBJpsiDiElectron);
183+
184+
// set random seed
185+
gen->readString("Random:setSeed on");
186+
uint random_seed;
187+
unsigned long long int random_value = 0;
188+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
189+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
190+
gen->readString(Form("Random:seed = %d", random_value % 900000001));
191+
192+
// print debug
193+
// gen->PrintDebug();
194+
195+
return gen;
196+
}
197+
198+
FairGenerator*
199+
GeneratorBeautyToPsiAndJpsi_EvtGenMidY(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "511;521;531;541;5112;5122;5232;5132;5332")
200+
{
201+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8NonPromptInjectedGapTriggeredDQ>();
202+
gen->setRapidity(rapidityMin, rapidityMax);
203+
gen->setPDG(5);
204+
gen->setRapidityHadron(rapidityMin,rapidityMax);
205+
gen->setHadronMultiplicity(1);
206+
TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBhadrons.cfg");
207+
gen->readFile(pathO2table.Data());
208+
gen->setConfigMBdecays(pathO2table);
209+
gen->setVerbose(verbose);
210+
211+
std::string spdg;
212+
TObjArray* obj = pdgs.Tokenize(";");
213+
gen->SetSizePdg(obj->GetEntriesFast());
214+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
215+
spdg = obj->At(i)->GetName();
216+
gen->AddPdg(std::stoi(spdg), i);
217+
printf("PDG %d \n", std::stoi(spdg));
218+
gen->addHadronPDGs(std::stoi(spdg));
219+
}
220+
gen->SetForceDecay(kEvtBPsiAndJpsiDiElectron);
221+
222+
// set random seed
223+
gen->readString("Random:setSeed on");
224+
uint random_seed;
225+
unsigned long long int random_value = 0;
226+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
227+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
228+
gen->readString(Form("Random:seed = %d", random_value % 900000001));
229+
230+
// print debug
231+
// gen->PrintDebug();
232+
233+
return gen;
234+
}
235+
236+
FairGenerator*
237+
GeneratorBplusToJpsiKaon_EvtGen(double rapidityMin = -1.5, double rapidityMax = 1.5, bool verbose = false, TString pdgs = "521")
238+
{
239+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorPythia8NonPromptInjectedGapTriggeredDQ>();
240+
gen->setRapidity(rapidityMin, rapidityMax);
241+
gen->setPDG(5);
242+
gen->addHadronPDGs(521);
243+
gen->setRapidityHadron(rapidityMin,rapidityMax);
244+
gen->setHadronMultiplicity(2);
245+
TString pathO2table = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffBplus.cfg");
246+
gen->readFile(pathO2table.Data());
247+
gen->setConfigMBdecays(pathO2table);
248+
gen->setVerbose(verbose);
249+
250+
std::string spdg;
251+
TObjArray* obj = pdgs.Tokenize(";");
252+
gen->SetSizePdg(obj->GetEntriesFast());
253+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
254+
spdg = obj->At(i)->GetName();
255+
gen->AddPdg(std::stoi(spdg), i);
256+
gen->addHadronPDGs(std::stoi(spdg));
257+
printf("PDG %d \n", std::stoi(spdg));
258+
}
259+
260+
TString pathO2 = gSystem->ExpandPathName("$O2DPG_ROOT/MC/config/PWGDQ/EvtGen/DecayTablesEvtgen");
261+
gen->SetDecayTable(Form("%s/BPLUSTOKAONJPSITOELE.DEC", pathO2.Data()));
262+
// print debug
263+
// gen->PrintDebug();
264+
// set random seed
265+
gen->readString("Random:setSeed on");
266+
uint random_seed;
267+
unsigned long long int random_value = 0;
268+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
269+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
270+
gen->readString(Form("Random:seed = %d", random_value % 900000001));
271+
272+
return gen;
273+
}
274+
275+
276+
Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
#include "FairGenerator.h"
2+
#include "Generators/GeneratorPythia8.h"
3+
#include "Pythia8/Pythia.h"
4+
#include "TRandom.h"
5+
#include "GeneratorPromptCharmonia.C"
6+
7+
#include <string>
8+
9+
using namespace o2::eventgen;
10+
using namespace Pythia8;
11+
12+
class GeneratorPythia8PromptInjectedGapTriggeredDQ : public o2::eventgen::GeneratorPythia8 {
13+
public:
14+
15+
/// default constructor
16+
GeneratorPythia8PromptInjectedGapTriggeredDQ() = default;
17+
18+
/// constructor
19+
GeneratorPythia8PromptInjectedGapTriggeredDQ(int inputTriggerRatio = 5, int gentype = 0) {
20+
21+
mGeneratedEvents = 0;
22+
mGeneratorParam = 0x0;
23+
mInverseTriggerRatio = inputTriggerRatio;
24+
switch (gentype) {
25+
case 0: // generate prompt jpsi + psi2s cocktail at midrapidity
26+
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToElectronEvtGen_pp13TeV();
27+
break;
28+
case 1: // generate prompt jpsi at midrapidity
29+
mGeneratorParam = (Generator*)GeneratorParamPromptJpsiToElectronEvtGen_pp13TeV("443");
30+
break;
31+
case 2: // generate prompt psi2S at midrapidity
32+
mGeneratorParam = (Generator*)GeneratorParamPromptPsiToElectronEvtGen_pp13TeV("100443");
33+
break;
34+
case 3: // generate prompt jpsi + psi2s cocktail at forward rapidity
35+
mGeneratorParam = (Generator*)GeneratorCocktailPromptCharmoniaToMuonEvtGen_pp13TeV();
36+
break;
37+
case 4: // generate prompt jpsi at forward rapidity
38+
mGeneratorParam = (Generator*)GeneratorParamPromptJpsiToMuonEvtGen_pp13TeV("443");
39+
break;
40+
case 5: // generate prompt psi2S at forward rapidity
41+
mGeneratorParam = (Generator*)GeneratorParamPromptPsiToMuonEvtGen_pp13TeV("100443");
42+
break;
43+
}
44+
mGeneratorParam->Init();
45+
}
46+
47+
/// Destructor
48+
~GeneratorPythia8PromptInjectedGapTriggeredDQ() = default;
49+
50+
protected:
51+
52+
Bool_t importParticles() override
53+
{
54+
GeneratorPythia8::importParticles();
55+
bool genOk = false;
56+
if (mGeneratedEvents % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
57+
bool genOk = false;
58+
while (!genOk){
59+
genOk = (mGeneratorParam->generateEvent() && mGeneratorParam->importParticles()) ? true : false ;
60+
}
61+
int originalSize = mParticles.size();
62+
for(int ipart=0; ipart < mGeneratorParam->getParticles().size(); ipart++){
63+
TParticle part = TParticle(mGeneratorParam->getParticles().at(ipart));
64+
if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
65+
if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
66+
if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
67+
mParticles.push_back(part);
68+
// encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C
69+
}
70+
mGeneratorParam->clearParticles();
71+
}
72+
73+
mGeneratedEvents++;
74+
return true;
75+
}
76+
77+
78+
private:
79+
Generator* mGeneratorParam;
80+
// Control gap-triggering
81+
unsigned long long mGeneratedEvents;
82+
int mInverseTriggerRatio;
83+
};
84+
85+
// Predefined generators:
86+
FairGenerator *GeneratorPythia8InjectedPromptCharmoniaGapTriggered(int inputTriggerRatio, int gentype) {
87+
auto myGen = new GeneratorPythia8PromptInjectedGapTriggeredDQ(inputTriggerRatio,gentype);
88+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
89+
myGen->readString("Random:setSeed on");
90+
myGen->readString("Random:seed " + std::to_string(seed));
91+
return myGen;
92+
}
Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
### The setup uses an external event generator
2+
### This part sets the path of the file and the function call to retrieve it
3+
4+
[GeneratorExternal]
5+
fileName = ${O2DPG_ROOT}/MC/config/PWGDQ/external/generator/generator_pythia8_NonPromptSignals_gaptriggered_dq.C
6+
funcName = GeneratorBplusToJpsiKaon_EvtGen()
7+
8+
### The external generator derives from GeneratorPythia8.
9+
### This part configures the bits of the interface: configuration and user hooks
10+
11+
[GeneratorPythia8]
12+
config = ${O2DPG_ROOT}/MC/config/common/pythia8/generator/pythia8_hf.cfg
13+
hooksFileName = ${O2DPG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
14+
hooksFuncName = pythia8_userhooks_bbbar(-1.5,1.5)
15+
16+
### The setup inhibits transport of primary particles which are produce at forward rapidity.
17+
### The settings below only transports particles in the barrel, which is currently defined as |eta| < 2
18+
19+
[Stack]
20+
transportPrimary = barrel

0 commit comments

Comments
 (0)