Skip to content

Commit a92d70b

Browse files
dsekihatdsekihat
authored
[MC/PWGEM] add HF and DY->dimuon MC in pp at 13.6 TeV (#2414)
1 parent fe572f3 commit a92d70b

11 files changed

Lines changed: 805 additions & 0 deletions

‎MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFMuon.C‎

Lines changed: 224 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,224 @@
1+
#include "Pythia8/Pythia.h"
2+
#include "Pythia8/HeavyIons.h"
3+
#include "FairGenerator.h"
4+
#include "FairPrimaryGenerator.h"
5+
#include "Generators/GeneratorPythia8.h"
6+
#include "TRandom3.h"
7+
#include "TParticlePDG.h"
8+
#include "TDatabasePDG.h"
9+
10+
#include <map>
11+
#include <unordered_set>
12+
//#include <utility> // for std::pair
13+
14+
using namespace Pythia8;
15+
16+
class GeneratorPythia8GapTriggeredHFLepton : public o2::eventgen::GeneratorPythia8
17+
{
18+
public:
19+
/// default constructor
20+
GeneratorPythia8GapTriggeredHFLepton() = default;
21+
22+
/// constructor
23+
GeneratorPythia8GapTriggeredHFLepton(TString configsignal, int quarkPdg = 4, int lInputTriggerRatio = 5, int lInputExternalID = 0)
24+
{
25+
26+
lGeneratedEvents = 0;
27+
lInverseTriggerRatio = lInputTriggerRatio;
28+
lExternalID = lInputExternalID;
29+
mQuarkPdg = quarkPdg;
30+
31+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
32+
33+
int offset = (int)(gRandom->Uniform(lInverseTriggerRatio)); // create offset to mitigate edge effects due to small number of events per job
34+
lGeneratedEvents += offset;
35+
36+
cout << "Initalizing extra PYTHIA object used to generate min-bias events..." << endl;
37+
TString pathconfigMB = gSystem->ExpandPathName("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_MB_gapevent.cfg");
38+
pythiaObjectMinimumBias.readFile(pathconfigMB.Data());
39+
pythiaObjectMinimumBias.readString("Random:setSeed on");
40+
pythiaObjectMinimumBias.readString("Random:seed " + std::to_string(seed));
41+
pythiaObjectMinimumBias.init();
42+
cout << "Initalization complete" << endl;
43+
cout << "Initalizing extra PYTHIA object used to generate signal events..." << endl;
44+
TString pathconfigSignal = gSystem->ExpandPathName(configsignal.Data());
45+
pythiaObjectSignal.readFile(pathconfigSignal.Data());
46+
pythiaObjectSignal.readString("Random:setSeed on");
47+
pythiaObjectSignal.readString("Random:seed " + std::to_string(seed));
48+
pythiaObjectSignal.init();
49+
cout << "Initalization complete" << endl;
50+
// flag the generators using type
51+
// addCocktailConstituent(type, "interesting");
52+
// addCocktailConstitent(0, "minbias");
53+
// Add Sub generators
54+
addSubGenerator(0, "default generator");
55+
addSubGenerator(1, "charm lepton");
56+
addSubGenerator(2, "beauty forced decay");
57+
addSubGenerator(3, "beauty no foced decay");
58+
}
59+
60+
/// Destructor
61+
~GeneratorPythia8GapTriggeredHFLepton() = default;
62+
63+
void addTriggerOnDaughter(int nb, int pdg)
64+
{
65+
mNbDaughter = nb;
66+
mPdgDaughter = pdg;
67+
};
68+
void setQuarkRapidity(float yMin, float yMax)
69+
{
70+
mQuarkRapidityMin = yMin;
71+
mQuarkRapidityMax = yMax;
72+
};
73+
void setDaughterRapidity(float yMin, float yMax)
74+
{
75+
mDaughterRapidityMin = yMin;
76+
mDaughterRapidityMax = yMax;
77+
};
78+
79+
protected:
80+
//__________________________________________________________________
81+
Bool_t generateEvent() override
82+
{
83+
/// reset event
84+
mPythia.event.reset();
85+
86+
// Simple straightforward check to alternate generators
87+
if (lGeneratedEvents % lInverseTriggerRatio == 0) {
88+
// Generate event of interest
89+
Bool_t lGenerationOK = kFALSE;
90+
while (!lGenerationOK) {
91+
if (pythiaObjectSignal.next()) {
92+
lGenerationOK = selectEvent(pythiaObjectSignal.event);
93+
}
94+
}
95+
mPythia.event = pythiaObjectSignal.event;
96+
notifySubGenerator(lExternalID);
97+
} else {
98+
// Generate minimum-bias event
99+
Bool_t lGenerationOK = kFALSE;
100+
while (!lGenerationOK) {
101+
lGenerationOK = pythiaObjectMinimumBias.next();
102+
}
103+
mPythia.event = pythiaObjectMinimumBias.event;
104+
notifySubGenerator(0);
105+
}
106+
107+
lGeneratedEvents++;
108+
// mPythia.next();
109+
110+
return true;
111+
}
112+
113+
bool selectEvent(const Pythia8::Event& event)
114+
{
115+
bool isGoodAtPartonLevel = false, isGoodAtDaughterLevel = (mPdgDaughter != 0) ? false : true;
116+
int nbDaughter = 0;
117+
for (auto iPart{0}; iPart < event.size(); ++iPart) {
118+
// search for Q-Qbar mother with at least one Q in rapidity window
119+
if (!isGoodAtPartonLevel) {
120+
auto daughterList = event[iPart].daughterList();
121+
bool hasQ = false, hasQbar = false, atSelectedY = false;
122+
for (auto iDau : daughterList) {
123+
if (event[iDau].id() == mQuarkPdg) {
124+
hasQ = true;
125+
}
126+
if (event[iDau].id() == -mQuarkPdg) {
127+
hasQbar = true;
128+
}
129+
if ((std::abs(event[iDau].id()) == mQuarkPdg) && (event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax))
130+
atSelectedY = true;
131+
}
132+
if (hasQ && hasQbar && atSelectedY) {
133+
isGoodAtPartonLevel = true;
134+
}
135+
}
136+
// search for mNbDaughter daughters of type mPdgDaughter in rapidity window
137+
if (!isGoodAtDaughterLevel) {
138+
int id = std::abs(event[iPart].id());
139+
float rap = event[iPart].y();
140+
if (id == mPdgDaughter) {
141+
int motherindexa = event[iPart].mother1();
142+
if (motherindexa > 0) {
143+
int idmother = std::abs(event[motherindexa].id());
144+
if (int(std::abs(idmother) / 100.) == 4 || int(std::abs(idmother) / 1000.) == 4 || int(std::abs(idmother) / 100.) == 5 || int(std::abs(idmother) / 1000.) == 5) {
145+
if (rap > mDaughterRapidityMin && rap < mDaughterRapidityMax) {
146+
nbDaughter++;
147+
if (nbDaughter >= mNbDaughter) isGoodAtDaughterLevel = true;
148+
}
149+
}
150+
}
151+
}
152+
}
153+
// we send the trigger
154+
if (isGoodAtPartonLevel && isGoodAtDaughterLevel) {
155+
return true;
156+
}
157+
}
158+
return false;
159+
};
160+
161+
private:
162+
// Interface to override import particles
163+
Pythia8::Event mOutputEvent;
164+
165+
// Properties of selection
166+
int mQuarkPdg;
167+
float mQuarkRapidityMin;
168+
float mQuarkRapidityMax;
169+
int mPdgDaughter;
170+
int mNbDaughter;
171+
float mDaughterRapidityMin;
172+
float mDaughterRapidityMax;
173+
174+
// Control gap-triggering
175+
Long64_t lGeneratedEvents;
176+
int lInverseTriggerRatio;
177+
// ID for different generators
178+
int lExternalID;
179+
180+
// Base event generators
181+
Pythia8::Pythia pythiaObjectMinimumBias; ///Minimum bias collision generator
182+
Pythia8::Pythia pythiaObjectSignal; ///Signal collision generator
183+
};
184+
185+
// Predefined generators:
186+
187+
// Charm-enriched forced decay
188+
FairGenerator* GeneratorPythia8GapTriggeredCharmLepton(int inputTriggerRatio, int inputExternalID, float yMin = -1.5, float yMax = 1.5)
189+
{
190+
auto myGen = new GeneratorPythia8GapTriggeredHFLepton("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_pp_cr2_forceddecayscharm_semimuonic.cfg", 4, inputTriggerRatio, inputExternalID);
191+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
192+
myGen->readString("Random:setSeed on");
193+
myGen->readString("Random:seed " + std::to_string(seed));
194+
myGen->setQuarkRapidity(yMin, yMax);
195+
myGen->addTriggerOnDaughter(2, 13);
196+
myGen->setDaughterRapidity(-6., -1.);
197+
return myGen;
198+
}
199+
200+
// Beauty-enriched forced decay
201+
FairGenerator* GeneratorPythia8GapTriggeredBeautyForcedDecays(int inputTriggerRatio, int inputExternalID, float yMin = -1.5, float yMax = 1.5)
202+
{
203+
auto myGen = new GeneratorPythia8GapTriggeredHFLepton("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar_forceddecayscharmbeauty_semimuonic.cfg", 5, inputTriggerRatio, inputExternalID);
204+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
205+
myGen->readString("Random:setSeed on");
206+
myGen->readString("Random:seed " + std::to_string(seed));
207+
myGen->setQuarkRapidity(yMin, yMax);
208+
myGen->addTriggerOnDaughter(2, 13);
209+
myGen->setDaughterRapidity(-6., -1.);
210+
return myGen;
211+
}
212+
213+
// Beauty-enriched no forced decay
214+
FairGenerator* GeneratorPythia8GapTriggeredBeautyNoForcedDecays(int inputTriggerRatio, int inputExternalID, float yMin = -1.5, float yMax = 1.5)
215+
{
216+
auto myGen = new GeneratorPythia8GapTriggeredHFLepton("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_bbbar.cfg", 5, inputTriggerRatio, inputExternalID);
217+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
218+
myGen->readString("Random:setSeed on");
219+
myGen->readString("Random:seed " + std::to_string(seed));
220+
myGen->setQuarkRapidity(yMin, yMax);
221+
myGen->addTriggerOnDaughter(2, 13);
222+
myGen->setDaughterRapidity(-6., -1.);
223+
return myGen;
224+
}

‎MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap1_pp13600GeV.ini‎

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/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C
3+
funcName = GeneratorPythia8GapTriggeredDYll(1, 1, 13, -6, -1, 2212, 2212, 13600.0)
4+
5+
[GeneratorPythia8]
6+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
7+
includePartonEvent = true

‎MC/config/PWGEM/ini/GeneratorHFGapTriggered_BeautyNoForcedDecay_Muon_pp13600GeV_Gap1.ini‎

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/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFMuon.C
3+
funcName = GeneratorPythia8GapTriggeredBeautyNoForcedDecays(1, 3, -6, -1)
4+
5+
[GeneratorPythia8]
6+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
7+
includePartonEvent = true

‎MC/config/PWGEM/ini/GeneratorHFGapTriggered_Charm_Muon_pp13600GeV_Gap1.ini‎

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/PWGEM/external/generator/Generator_pythia8_GapTriggered_HFMuon.C
3+
funcName = GeneratorPythia8GapTriggeredCharmLepton(1, 1, -6, -1)
4+
5+
[GeneratorPythia8]
6+
config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg
7+
includePartonEvent = true

‎MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap1_pp13600GeV.C‎

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
// Check that file exists, can be opened and has the correct tree
4+
TFile file(path.c_str(), "READ");
5+
if (file.IsZombie())
6+
{
7+
std::cerr << "Cannot open ROOT file " << path << "\n";
8+
return 1;
9+
}
10+
auto tree = (TTree *)file.Get("o2sim");
11+
if (!tree)
12+
{
13+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
14+
return 1;
15+
}
16+
std::vector<o2::MCTrack> *tracks{};
17+
tree->SetBranchAddress("MCTrack", &tracks);
18+
19+
// Check if all events are filled
20+
auto nEvents = tree->GetEntries();
21+
for (Long64_t i = 0; i < nEvents; ++i)
22+
{
23+
tree->GetEntry(i);
24+
if (tracks->empty())
25+
{
26+
std::cerr << "Empty entry found at event " << i << "\n";
27+
return 1;
28+
}
29+
}
30+
// check if each event has at least two oxygen ions
31+
for (int i = 0; i < nEvents; i++)
32+
{
33+
auto check = tree->GetEntry(i);
34+
int count = 0;
35+
for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack)
36+
{
37+
auto track = tracks->at(idxMCTrack);
38+
if (track.GetPdgCode() == 2212)
39+
{
40+
count++;
41+
}
42+
}
43+
if (count < 2)
44+
{
45+
std::cerr << "Event " << i << " has less than 2 oxygen ions\n";
46+
return 1;
47+
}
48+
}
49+
50+
return 0;
51+
}

0 commit comments

Comments
 (0)