Skip to content

Commit fa2d14d

Browse files
rashiguptrashigupt
authored
Create generator_pythia8_gaptriggered_nonhfe.C
Add generator_pythia8_gaptriggered_nonhfe.C file for enance pion and eta
1 parent 41fabe0 commit fa2d14d

File tree

1 file changed

+218
-0
lines changed

1 file changed

+218
-0
lines changed

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

Lines changed: 218 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
#include "FairGenerator.h"
2+
3+
#include "TDatabasePDG.h"
4+
#include "TRandom.h"
5+
6+
#include "Generators/GeneratorPythia8.h"
7+
#include "Generators/GeneratorPythia8Param.h"
8+
#include "Pythia8/Pythia.h"
9+
#include <fairlogger/Logger.h>
10+
11+
#include <string>
12+
#include <vector>
13+
14+
using namespace Pythia8;
15+
16+
class GeneratorPythia8GapTriggeredPionEta : public o2::eventgen::GeneratorPythia8
17+
{
18+
public:
19+
/// default constructor
20+
GeneratorPythia8GapTriggeredPionEta() = default;
21+
GeneratorPythia8GapTriggeredPionEta(int inputTriggerRatio = 5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {})
22+
{
23+
mGeneratedEvents = 0;
24+
mInverseTriggerRatio = inputTriggerRatio;
25+
mQuarkRapidityMin = -1.5;
26+
mQuarkRapidityMax = 1.5;
27+
mHadRapidityMin = -1.5;
28+
mHadRapidityMax = 1.5;
29+
30+
mQuarkPdg = 0;
31+
mHadronPdg = 0;
32+
mQuarkPdgList = quarkPdgList;
33+
mHadronPdgList = hadronPdgList;
34+
Print();
35+
}
36+
37+
/// Destructor
38+
GeneratorPythia8GapTriggeredPionEta() = default;
39+
/// Print the input
40+
void Print()
41+
{
42+
LOG(info) << "********** GeneratorPythia8GapTriggeredHF configuration dump **********";
43+
LOG(info) << Form("* Trigger ratio: %d", mInverseTriggerRatio);
44+
LOG(info) << Form("* Quark pdg: %d", mQuarkPdg);
45+
LOG(info) << Form("* Quark rapidity: %f - %f", mQuarkRapidityMin, mQuarkRapidityMax);
46+
LOG(info) << Form("* Hadron pdg: %d", mHadronPdg);
47+
LOG(info) << Form("* Hadron rapidity: %f - %f", mHadRapidityMin, mHadRapidityMax);
48+
LOG(info) << Form("* Quark pdg list: ");
49+
for (auto pdg : mQuarkPdgList) {
50+
LOG(info) << Form("* %d ", pdg);
51+
}
52+
LOG(info) << Form("* Hadron pdg list: ");
53+
for (auto pdg : mHadronPdgList) {
54+
LOG(info) << Form("* %d ", pdg);
55+
}
56+
LOG(info) << "***********************************************************************";
57+
}
58+
59+
bool Init() override
60+
{
61+
addSubGenerator(0, "Minimum bias");
62+
addSubGenerator(1, "Down injected");
63+
addSubGenerator(2, "Up injected");
64+
addSubGenerator(3, "Strange injected");
65+
66+
return o2::eventgen::GeneratorPythia8::Init();
67+
}
68+
void setQuarkRapidity(float yMin, float yMax)
69+
{
70+
mQuarkRapidityMin = yMin;
71+
mQuarkRapidityMax = yMax;
72+
};
73+
void setHadronRapidity(float yMin, float yMax)
74+
{
75+
mHadRapidityMin = yMin;
76+
mHadRapidityMax = yMax;
77+
};
78+
void setUsedSeed(unsigned int seed)
79+
{
80+
mUsedSeed = seed;
81+
};
82+
unsigned int getUsedSeed() const
83+
{
84+
return mUsedSeed;
85+
};
86+
87+
protected:
88+
//__________________________________________________________________
89+
bool generateEvent() override
90+
{
91+
92+
// Simple straightforward check to alternate generators
93+
if (mGeneratedEvents % mInverseTriggerRatio == 0) {
94+
int nInjectedEvents = mGeneratedEvents / mInverseTriggerRatio;
95+
// Alternate quarks if enabled (with the same ratio)
96+
if (mQuarkPdgList.size() >= 1) {
97+
int iQuark = nInjectedEvents % mQuarkPdgList.size();
98+
mQuarkPdg = mQuarkPdgList[iQuark];
99+
LOG(debug) << "SELECTED quark: " << mQuarkPdgList[iQuark];
100+
}
101+
// Alternate hadrons if enabled (with the same ratio)
102+
if (mHadronPdgList.size() >= 1) {
103+
int iHadron = (nInjectedEvents / std::max(mQuarkPdgList.size(), 1ul)) % mHadronPdgList.size();
104+
mHadronPdg = mHadronPdgList[iHadron];
105+
LOG(debug) << "SELECTED hadron: " << mHadronPdgList[iHadron];
106+
}
107+
108+
// Generate event of interest
109+
bool genOk = false;
110+
while (!genOk) {
111+
if (GeneratorPythia8::generateEvent()) {
112+
genOk = selectEvent();
113+
}
114+
}
115+
notifySubGenerator(mQuarkPdg);
116+
} else {
117+
// Generate minimum-bias event
118+
bool genOk = false;
119+
while (!genOk) {
120+
genOk = GeneratorPythia8::generateEvent();
121+
}
122+
notifySubGenerator(0);
123+
}
124+
125+
mGeneratedEvents++;
126+
127+
return true;
128+
}
129+
bool selectEvent()
130+
{
131+
132+
bool isGoodAtPartonLevel{mQuarkPdgList.size() == 0};
133+
bool isGoodAtHadronLevel{mHadronPdgList.size() == 0};
134+
135+
for (auto iPart{0}; iPart < mPythia.event.size(); ++iPart) {
136+
// search for Q-Qbar mother with at least one Q in rapidity window
137+
if (!isGoodAtPartonLevel) {
138+
auto daughterList = mPythia.event[iPart].daughterList();
139+
bool hasQ = false, hasQbar = false, atSelectedY = false;
140+
for (auto iDau : daughterList) {
141+
if (mPythia.event[iDau].id() == mQuarkPdg) {
142+
hasQ = true;
143+
if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax)) {
144+
atSelectedY = true;
145+
}
146+
}
147+
if (mPythia.event[iDau].id() == -mQuarkPdg) {
148+
hasQbar = true;
149+
if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax)) {
150+
atSelectedY = true;
151+
}
152+
}
153+
}
154+
if (hasQ && hasQbar && atSelectedY) {
155+
isGoodAtPartonLevel = true;
156+
}
157+
}
158+
// search for hadron in rapidity window
159+
if (!isGoodAtHadronLevel) {
160+
int id = std::abs(mPythia.event[iPart].id());
161+
float rap = mPythia.event[iPart].y();
162+
if (id == mHadronPdg && rap > mHadRapidityMin && rap < mHadRapidityMax) {
163+
isGoodAtHadronLevel = true;
164+
}
165+
}
166+
// we send the trigger immediately (if there are no particles to replace, that can be different from the trigger ones)
167+
if (isGoodAtPartonLevel && isGoodAtHadronLevel) {
168+
LOG(debug) << "EVENT SELECTED: Found particle " << mPythia.event[iPart].id() << " at rapidity " << mPythia.event[iPart].y() << "\n";
169+
return true;
170+
}
171+
}
172+
// we send the trigger
173+
if (isGoodAtPartonLevel && isGoodAtHadronLevel) {
174+
return true;
175+
}
176+
177+
return false;
178+
};
179+
180+
private:
181+
// Interface to override import particles
182+
Pythia8::Event mOutputEvent;
183+
184+
// Properties of selection
185+
int mQuarkPdg;
186+
float mQuarkRapidityMin;
187+
float mQuarkRapidityMax;
188+
int mHadronPdg;
189+
float mHadRapidityMin;
190+
float mHadRapidityMax;
191+
unsigned int mUsedSeed;
192+
193+
// Control gap-triggering
194+
unsigned long long mGeneratedEvents;
195+
int mInverseTriggerRatio;
196+
197+
// Control alternate trigger on charm and beauty quarks
198+
std::vector<int> mQuarkPdgList = {};
199+
200+
// Control alternate trigger on different hadrons
201+
std::vector<int> mHadronPdgList = {};
202+
};
203+
// Predefined generators:
204+
// Predefined generators:
205+
// Charm-enriched
206+
FairGenerator* GeneratorPythia8GapTriggeredPionAndEta(int inputTriggerRatio, 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 = {})
207+
{
208+
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{1, 2, 3}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
209+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
210+
myGen->setUsedSeed(seed);
211+
myGen->readString("Random:setSeed on");
212+
myGen->readString("Random:seed " + std::to_string(seed));
213+
myGen->setQuarkRapidity(yQuarkMin, yQuarkMax);
214+
if (hadronPdgList.size() != 0) {
215+
myGen->setHadronRapidity(yHadronMin, yHadronMax);
216+
}
217+
return myGen;
218+
}

0 commit comments

Comments
 (0)