Skip to content

Commit 54b415f

Browse files
committed
Files required for an O2DPG production of prompt-photon MC, with a gap trigger
1 parent d59baad commit 54b415f

File tree

4 files changed

+314
-0
lines changed

4 files changed

+314
-0
lines changed
Lines changed: 279 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,279 @@
1+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
2+
///#include "FairGenerator.h"
3+
//#include "Generators/GeneratorPythia8.h"
4+
#include "Pythia8/Pythia.h"
5+
#include "MC/config/PWGGAJE/hooks/prompt_gamma_hook.C"
6+
//#include "TRandom.h"
7+
//#include <fairlogger/Logger.h>
8+
//
9+
//#include <string>
10+
//#include <vector>
11+
12+
// Prompt-photon custom event generator
13+
// that alternates between 2 gun generators.
14+
// set up to inject MB events alongside prompt-photon events
15+
// in 'MB-gap' mode.
16+
// The number of MB events injected, and the PYTHIA config
17+
// for each event type is defined by the user in the .ini
18+
// generator file (e.g. GeneratorJE_gapgen5_hook.ini)
19+
//
20+
// Author: Adrian Fereydon Nassirpour (adrian.fereydon.nassirpour@cern.ch), based on code from Jaime Norman (jaime.norman@cern.ch)
21+
22+
namespace o2
23+
{
24+
namespace eventgen
25+
{
26+
27+
using namespace Pythia8;
28+
29+
30+
/// A very simple gap generator alternating between 2 different particle guns
31+
class GeneratorPythia8GapGenJEPhoton : public o2::eventgen::GeneratorPythia8
32+
{
33+
public:
34+
/// default constructor
35+
GeneratorPythia8GapGenJEPhoton(int inputTriggerRatio = 5,std::string pathMB = "",std::string pathSignal = "") {
36+
37+
mGeneratedEvents = 0;
38+
mInverseTriggerRatio = inputTriggerRatio;
39+
40+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
41+
42+
cout << "Initalizing extra PYTHIA object used to generate min-bias events..." << endl;
43+
TString pathconfigMB = gSystem->ExpandPathName(TString(pathMB));
44+
pythiaObjectMinimumBias.readFile(pathconfigMB.Data());
45+
pythiaObjectMinimumBias.readString("Random:setSeed on");
46+
pythiaObjectMinimumBias.readString("Random:seed " + std::to_string(seed));
47+
pythiaObjectMinimumBias.init();
48+
cout << "Initalization complete" << endl;
49+
cout << "Initalizing extra PYTHIA object used to generate signal events..." << endl;
50+
TString pathconfigSignal = gSystem->ExpandPathName(TString(pathSignal));
51+
pythiaObjectSignal.readFile(pathconfigSignal.Data());
52+
pythiaObjectSignal.readString("Random:setSeed on");
53+
pythiaObjectSignal.readString("Random:seed " + std::to_string(seed));
54+
// load jet hook to ensure at least one prompt photon is within detector acceptance
55+
Pythia8::UserHooks *hook = pythia8_userhooks_promptgamma();
56+
pythiaObjectSignal.setUserHooksPtr(std::shared_ptr<Pythia8::UserHooks>(hook));
57+
pythiaObjectSignal.init();
58+
cout << "Initalization complete" << endl;
59+
// Add Sub generators
60+
addSubGenerator(0, "MB generator");
61+
addSubGenerator(1, "jet-jet generator");
62+
}
63+
64+
65+
/// Destructor
66+
~GeneratorPythia8GapGenJEPhoton() = default;
67+
68+
void setUsedSeed(unsigned int seed)
69+
{
70+
mUsedSeed = seed;
71+
};
72+
unsigned int getUsedSeed() const
73+
{
74+
return mUsedSeed;
75+
};
76+
77+
bool generateEvent() override
78+
{
79+
80+
// Simple straightforward check to alternate generators
81+
mPythia.event.reset();
82+
83+
if (mGeneratedEvents % mInverseTriggerRatio == 0) {
84+
LOG(info) << "Event " << mGeneratedEvents << ", generate signal event";
85+
// Generate event of interest
86+
Bool_t mGenerationOK = kFALSE;
87+
while (!mGenerationOK) {
88+
mGenerationOK = pythiaObjectSignal.next();
89+
}
90+
mPythia.event = pythiaObjectSignal.event;
91+
setEventHeaderProperties(pythiaObjectSignal);
92+
LOG(info) << "--- Print info properties custom...";
93+
printEventHeaderProperties(pythiaObjectSignal);
94+
notifySubGenerator(1);
95+
}
96+
else {
97+
LOG(info) << "Event " << mGeneratedEvents << ", generate mb event";
98+
// Generate minimum-bias event
99+
Bool_t mGenerationOK = kFALSE;
100+
while (!mGenerationOK) {
101+
mGenerationOK = pythiaObjectMinimumBias.next();
102+
}
103+
mPythia.event = pythiaObjectMinimumBias.event;
104+
setEventHeaderProperties(pythiaObjectMinimumBias);
105+
LOG(info) << "--- Print info properties custom...";
106+
printEventHeaderProperties(pythiaObjectMinimumBias);
107+
notifySubGenerator(0);
108+
}
109+
mGeneratedEvents++;
110+
return true;
111+
}
112+
113+
// for testing
114+
void printEventHeaderProperties (Pythia8::Pythia &pythiaObject) {
115+
LOG(info) << "Info name = " << pythiaObject.info.name();
116+
LOG(info) << "Info code = " << pythiaObject.info.code();
117+
LOG(info) << "Info weight = " << pythiaObject.info.weight();
118+
LOG(info) << "Info id1pdf = " << pythiaObject.info.id1pdf();
119+
LOG(info) << "Info id2pdf = " << pythiaObject.info.id2pdf();
120+
121+
LOG(info) << "Info x1pdf = " << pythiaObject.info.x1pdf();
122+
LOG(info) << "Info x2pdf = " << pythiaObject.info.x2pdf();
123+
LOG(info) << "Info QFac = " << pythiaObject.info.QFac();
124+
LOG(info) << "Info pdf1 = " << pythiaObject.info.pdf1();
125+
LOG(info) << "Info pdf2 = " << pythiaObject.info.pdf2();
126+
127+
// Set cross section
128+
LOG(info) << "Info sigmaGen = " << pythiaObject.info.sigmaGen();
129+
LOG(info) << "Info sigmaErr = " << pythiaObject.info.sigmaErr();
130+
131+
// Set event scale and nMPI
132+
LOG(info) << "Info QRen = " << pythiaObject.info.QRen();
133+
LOG(info) << "Info nMPI = " << pythiaObject.info.nMPI();
134+
135+
// Set accepted and attempted values
136+
LOG(info) << "Info accepted = " << pythiaObject.info.nAccepted();
137+
LOG(info) << "Info attempted = " << pythiaObject.info.nTried();
138+
139+
// Set weights (overrides cross-section for each weight)
140+
size_t iw = 0;
141+
auto xsecErr = pythiaObject.info.weightContainerPtr->getTotalXsecErr();
142+
for (auto w : pythiaObject.info.weightContainerPtr->getTotalXsec()) {
143+
std::string post = (iw == 0 ? "" : "_" + std::to_string(iw));
144+
LOG(info) << "Info weight by index = " << pythiaObject.info.weightValueByIndex(iw);
145+
iw++;
146+
}
147+
148+
}
149+
150+
// in order to save the event weight we need to override the following function
151+
// from the inherited o2::eventgen::GeneratorPythia8 class. The event header properties
152+
// are created as members of this class, and are set using the active event generator
153+
// (MB or jet-jet), then propagated to the event header
154+
void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override {
155+
/** update header **/
156+
using Key = o2::dataformats::MCInfoKeys;
157+
158+
eventHeader->putInfo<std::string>(Key::generator, "pythia8");
159+
eventHeader->putInfo<int>(Key::generatorVersion, PYTHIA_VERSION_INTEGER);
160+
eventHeader->putInfo<std::string>(Key::processName, name);
161+
eventHeader->putInfo<int>(Key::processCode, code);
162+
eventHeader->putInfo<float>(Key::weight, weight);
163+
164+
// Set PDF information
165+
eventHeader->putInfo<int>(Key::pdfParton1Id, id1pdf);
166+
eventHeader->putInfo<int>(Key::pdfParton2Id, id2pdf);
167+
eventHeader->putInfo<float>(Key::pdfX1, x1pdf);
168+
eventHeader->putInfo<float>(Key::pdfX2, x2pdf);
169+
eventHeader->putInfo<float>(Key::pdfScale, QFac);
170+
eventHeader->putInfo<float>(Key::pdfXF1, pdf1);
171+
eventHeader->putInfo<float>(Key::pdfXF2, pdf2);
172+
173+
// Set cross section
174+
eventHeader->putInfo<float>(Key::xSection, sigmaGen * 1e9);
175+
eventHeader->putInfo<float>(Key::xSectionError, sigmaErr * 1e9);
176+
177+
// Set event scale and nMPI
178+
eventHeader->putInfo<float>(Key::eventScale, QRen);
179+
eventHeader->putInfo<int>(Key::mpi, nMPI);
180+
181+
// Set accepted and attempted events
182+
eventHeader->putInfo<int>(Key::acceptedEvents, accepted);
183+
eventHeader->putInfo<int>(Key::attemptedEvents, attempted);
184+
185+
LOG(info) << "--- updated header weight = " << weight;
186+
187+
// The following is also set in the base class updateHeader function
188+
// but as far as I can tell, there is no Xsec weight in the default
189+
// header so this is not copied over for now
190+
191+
//size_t iw = 0;
192+
//auto xsecErr = info.weightContainerPtr->getTotalXsecErr();
193+
//for (auto w : info.weightContainerPtr->getTotalXsec()) {
194+
// std::string post = (iw == 0 ? "" : "_" + std::to_string(iw));
195+
// eventHeader->putInfo<float>(Key::weight + post, info.weightValueByIndex(iw));
196+
// eventHeader->putInfo<float>(Key::xSection + post, w * 1e9);
197+
// eventHeader->putInfo<float>(Key::xSectionError + post, xsecErr[iw] * 1e9);
198+
// iw++;
199+
//}
200+
}
201+
202+
void setEventHeaderProperties (Pythia8::Pythia &pythiaObject) {
203+
204+
auto& info = pythiaObject.info;
205+
206+
name = info.name();
207+
code = info.code();
208+
weight = info.weight();
209+
// Set PDF information
210+
id1pdf = info.id1pdf();
211+
id2pdf = info.id2pdf();
212+
x1pdf = info.x1pdf();
213+
x2pdf = info.x2pdf();
214+
QFac = info.QFac();
215+
pdf1 = info.pdf1();
216+
pdf2 = info.pdf2();
217+
// Set cross section
218+
sigmaGen = info.sigmaGen();
219+
sigmaErr = info.sigmaErr();
220+
// Set event scale and nMPI
221+
QRen = info.QRen();
222+
nMPI = info.nMPI();
223+
// Set accepted and attempted events
224+
accepted = info.nAccepted();
225+
attempted = info.nTried();
226+
}
227+
228+
private:
229+
// Interface to override import particles
230+
Pythia8::Event mOutputEvent;
231+
232+
// Properties of selection
233+
unsigned int mUsedSeed;
234+
235+
// Control gap-triggering
236+
unsigned long long mGeneratedEvents;
237+
int mInverseTriggerRatio;
238+
239+
// Handling generators
240+
Pythia8::Pythia pythiaObjectMinimumBias;
241+
Pythia8::Pythia pythiaObjectSignal;
242+
243+
// header info - needed to save event properties
244+
std::string name;
245+
int code;
246+
float weight;
247+
// PDF information
248+
int id1pdf;
249+
int id2pdf;
250+
float x1pdf;
251+
float x2pdf;
252+
float QFac;
253+
float pdf1;
254+
float pdf2;
255+
// cross section
256+
float sigmaGen;
257+
float sigmaErr;
258+
// event scale and nMPI
259+
float QRen;
260+
int nMPI;
261+
// accepted and attempted events
262+
int accepted;
263+
int attempted;
264+
};
265+
266+
} // namespace eventgen
267+
} // namespace o2
268+
269+
/** generator instance and settings **/
270+
271+
FairGenerator* getGeneratorPythia8GapGenJEPhoton(int inputTriggerRatio = 5, std::string pathMB = "",std::string pathSignal = "") {
272+
auto myGen = new o2::eventgen::GeneratorPythia8GapGenJEPhoton(inputTriggerRatio, pathMB, pathSignal);
273+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
274+
myGen->setUsedSeed(seed);
275+
myGen->readString("Random:setSeed on");
276+
myGen->readString("Random:seed " + std::to_string(seed));
277+
myGen->readString("HardQCD:all = on");
278+
return myGen;
279+
}
Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/generator_pythia8_gaptrigger_promptphotons_hook.C
3+
funcName = getGeneratorPythia8GapGenJEPhoton(2,"${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_minbias.cfg","${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_promptphoton.cfg")
4+
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_promptphoton.cfg
7+
includePartonEvent=true
Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/generator_pythia8_gaptrigger_jets_hook.C
3+
funcName = getGeneratorPythia8GapGenJE(2,"${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_minbias.cfg","${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_promptphoton.cfg")
4+
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_promptphoton.cfg
7+
includePartonEvent=true
Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
# 2 -> 2 prompt photon production, oversampling by pThat^4
2+
3+
### beams
4+
Beams:idA = 2212 # proton
5+
Beams:idB = 2212 # proton
6+
Beams:eCM = 13600. # GeV
7+
8+
### processes
9+
SoftQCD:inelastic = off
10+
HardQCD:all = off
11+
PromptPhoton:all = on
12+
13+
### decays
14+
ParticleDecays:limitTau0 = on
15+
ParticleDecays:tau0Max = 10.
16+
17+
### phase space cuts
18+
PhaseSpace:pTHatMin = 5
19+
PhaseSpace:pTHatMax = 600
20+
PhaseSpace:bias2Selection = on
21+
PhaseSpace:bias2SelectionPow = 4

0 commit comments

Comments
 (0)