Skip to content

Commit 996d3b4

Browse files
jaimenormanjaimenorman
authored
add configuration and custom generator for jet-jet gap trigger setup (#1842)
* add configuration and custom generator for jet-jet gap trigger setup * update O2DPG_ROOT * add test for ini file
1 parent 324594a commit 996d3b4

File tree

5 files changed

+377
-0
lines changed

5 files changed

+377
-0
lines changed

‎MC/config/PWGGAJE/external/generator/generator_pythia8_gaptrigger_jets_hook.C‎

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

‎MC/config/PWGGAJE/ini/GeneratorJE_gapgen5_hook.ini‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
### jet-jet production with MB Gap 5
2+
### The external generator derives from GeneratorPythia8.
3+
[GeneratorExternal]
4+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/generator_pythia8_gaptrigger_jets_hook.C
5+
funcName = getGeneratorPythia8GapGenJE(5,"${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_minbias.cfg","${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg")

‎MC/config/PWGGAJE/ini/tests/GeneratorJE_gapgen5_hook.C‎

Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
4+
float ratioTrigger = 1./5; // one event triggered out of 5
5+
6+
7+
TFile file(path.c_str(), "READ");
8+
if (file.IsZombie()) {
9+
std::cerr << "Cannot open ROOT file " << path << "\n";
10+
return 1;
11+
}
12+
13+
auto tree = (TTree *)file.Get("o2sim");
14+
std::vector<o2::MCTrack> *tracks{};
15+
tree->SetBranchAddress("MCTrack", &tracks);
16+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
17+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
18+
19+
int nEventsMB{}, nEventsJetJet{};
20+
float sumWeightsMB{}, sumWeightsJetJet{};
21+
int sumTracks{};
22+
auto nEvents = tree->GetEntries();
23+
24+
for (int i = 0; i < nEvents; i++) {
25+
tree->GetEntry(i);
26+
27+
// check subgenerator information and event weights
28+
if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
29+
bool isValid = false;
30+
int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
31+
if (eventHeader->hasInfo(o2::dataformats::MCInfoKeys::weight)) {
32+
float weight = eventHeader->getInfo<float>(o2::dataformats::MCInfoKeys::weight,isValid);
33+
if (subGeneratorId == 0) {
34+
nEventsMB++;
35+
sumWeightsMB += weight;
36+
}
37+
else if (subGeneratorId == 1) {
38+
nEventsJetJet++;
39+
sumWeightsJetJet += weight;
40+
}
41+
}
42+
}
43+
sumTracks += tracks->size();
44+
}
45+
46+
std::cout << "--------------------------------\n";
47+
std::cout << "# Events: " << nEvents << "\n";
48+
std::cout << "# MB events: " << nEventsMB << "\n";
49+
std::cout << " sum of weights for MB events: " << sumWeightsMB << "\n";
50+
std::cout << "# Jet-jet events " << nEventsJetJet << "\n";
51+
std::cout << " sum of weights jet-jet events: " << sumWeightsJetJet << "\n";
52+
std::cout << "# tracks summed over all events (jet-jet + MB): " << sumTracks << "\n";
53+
54+
if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || nEventsMB > nEvents * (1 - ratioTrigger) * 1.05) { // we put some tolerance since the number of generated events is small
55+
std::cerr << "Number of generated MB events different than expected\n";
56+
return 1;
57+
}
58+
if (nEventsJetJet < nEvents * ratioTrigger * 0.95 || nEventsJetJet > nEvents * ratioTrigger * 1.05) {
59+
std::cerr << "Number of jet-jet generated events different than expected\n";
60+
return 1;
61+
}
62+
if(nEventsMB < sumWeightsMB * 0.95 || nEventsMB > sumWeightsMB * 1.05) {
63+
std::cerr << "Weights of MB events do not = 1 as expected\n";
64+
return 1;
65+
}
66+
if(sumTracks < 1) {
67+
std::cerr << "No tracks in simulated events\n";
68+
return 1;
69+
}
70+
return 0;
71+
}

‎MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# 2 -> 2 jet 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 = on
11+
12+
### decays
13+
ParticleDecays:limitTau0 = on
14+
ParticleDecays:tau0Max = 10.
15+
16+
### phase space cuts
17+
PhaseSpace:pTHatMin = 5
18+
PhaseSpace:pTHatMax = 600
19+
PhaseSpace:bias2Selection = on
20+
PhaseSpace:bias2SelectionPow = 4

‎MC/config/PWGGAJE/pythia8/generator/pythia8_minbias.cfg‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
# min. bias events
2+
3+
### beams
4+
Beams:idA = 2212 # proton
5+
Beams:idB = 2212 # proton
6+
Beams:eCM = 13600. # GeV
7+
8+
### processes
9+
SoftQCD:inelastic = on
10+
11+
### decays
12+
ParticleDecays:limitTau0 = on
13+
ParticleDecays:tau0Max = 10.

0 commit comments

Comments
 (0)