Skip to content

Commit 88d0f46

Browse files
motomioyaalcaliva
motomioya
authored and
alcaliva
committed
Add files to simulate HFs to dimuons with trigger gap (#1820)
* Add files to simulate HFs to dimuons with trigger gap * Modified requirements in test file * Modify O2DPG_ROOT to O2O2DPG_MC_CONFIG_ROOT (cherry picked from commit d1b9a1f)
1 parent 45374e7 commit 88d0f46

5 files changed

+385
-0
lines changed

‎MC/config/PWGDQ/external/generator/GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq.C‎

Lines changed: 203 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,203 @@
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_MC_CONFIG_ROOT/MC/config/PWGDQ/EvtGen)
7+
#include "GeneratorEvtGen.C"
8+
9+
#include <string>
10+
11+
using namespace o2::eventgen;
12+
13+
namespace o2
14+
{
15+
namespace eventgen
16+
{
17+
18+
class GeneratorHFToMu_EvtGenFwdY_gaptriggered : public o2::eventgen::GeneratorPythia8 {
19+
public:
20+
21+
/// constructor
22+
GeneratorHFToMu_EvtGenFwdY_gaptriggered(int inputTriggerRatio = 4) {
23+
24+
mGeneratedEvents = 0;
25+
mInverseTriggerRatio = inputTriggerRatio;
26+
// define minimum bias event generator
27+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
28+
TString pathconfigMB = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg");
29+
pythiaMBgen.readFile(pathconfigMB.Data());
30+
pythiaMBgen.readString("Random:setSeed on");
31+
pythiaMBgen.readString("Random:seed " + std::to_string(seed));
32+
mConfigMBdecays = "";
33+
mPDG = 4;
34+
mRapidityMin = -1.;
35+
mRapidityMax = 1.;
36+
mHadronMultiplicity = -1;
37+
mHadronRapidityMin = -1.;
38+
mHadronRapidityMax = 1.;
39+
mVerbose = false;
40+
}
41+
42+
/// Destructor
43+
~GeneratorHFToMu_EvtGenFwdY_gaptriggered() = default;
44+
45+
void setPDG(int val) { mPDG = val; };
46+
void addHadronPDGs(int pdg) { mHadronsPDGs.push_back(pdg); };
47+
void setHadronMultiplicity(int val) { mHadronMultiplicity = val; };
48+
void setRapidity(double valMin, double valMax)
49+
{
50+
mRapidityMin = valMin;
51+
mRapidityMax = valMax;
52+
};
53+
54+
void setRapidityHadron(double valMin, double valMax)
55+
{
56+
mHadronRapidityMin = valMin;
57+
mHadronRapidityMax = valMax;
58+
};
59+
60+
void setConfigMBdecays(TString val){mConfigMBdecays = val;}
61+
62+
void setVerbose(bool val) { mVerbose = val; };
63+
64+
protected:
65+
66+
Bool_t generateEvent() override
67+
{
68+
// reset event
69+
bool genOk = false;
70+
if (mGeneratedEvents % mInverseTriggerRatio == 0){
71+
bool ancestor = false;
72+
while (! (genOk && ancestor) ){
73+
/// reset event
74+
mPythia.event.reset();
75+
genOk = GeneratorPythia8::generateEvent();
76+
// find the q-qbar ancestor
77+
ancestor = findHeavyQuarkPair(mPythia.event);
78+
}
79+
} else {
80+
/// reset event
81+
pythiaMBgen.event.reset();
82+
while (!genOk) {
83+
genOk = pythiaMBgen.next();
84+
}
85+
mPythia.event = pythiaMBgen.event;
86+
}
87+
mGeneratedEvents++;
88+
if (mVerbose) mOutputEvent.list();
89+
return true;
90+
}
91+
92+
Bool_t Init() override
93+
{
94+
if(mConfigMBdecays.Contains("cfg")) pythiaMBgen.readFile(mConfigMBdecays.Data());
95+
GeneratorPythia8::Init();
96+
pythiaMBgen.init();
97+
return true;
98+
}
99+
100+
// search for q-qbar mother with at least one q in a selected rapidity window
101+
bool findHeavyQuarkPair(Pythia8::Event& event)
102+
{
103+
int countH[mHadronsPDGs.size()]; for(int ipdg=0; ipdg < mHadronsPDGs.size(); ipdg++) countH[ipdg]=0;
104+
bool hasq = false, hasqbar = false, atSelectedY = false, isOkAtPartonicLevel = false;
105+
for (int ipa = 0; ipa < event.size(); ++ipa) {
106+
107+
if(!isOkAtPartonicLevel){
108+
auto daughterList = event[ipa].daughterList();
109+
hasq = false; hasqbar = false; atSelectedY = false;
110+
for (auto ida : daughterList) {
111+
if (event[ida].id() == mPDG)
112+
hasq = true;
113+
if (event[ida].id() == -mPDG)
114+
hasqbar = true;
115+
if ((event[ida].y() > mRapidityMin) && (event[ida].y() < mRapidityMax))
116+
atSelectedY = true;
117+
}
118+
if (hasq && hasqbar && atSelectedY) isOkAtPartonicLevel = true;
119+
}
120+
121+
if( (mHadronMultiplicity <= 0) && isOkAtPartonicLevel) return true; // no selection at hadron level
122+
123+
/// check at hadron level if needed
124+
int ipdg=0;
125+
for (auto& pdgVal : mHadronsPDGs){
126+
if ( (TMath::Abs(event[ipa].id()) == pdgVal) && (event[ipa].y() > mHadronRapidityMin) && (event[ipa].y() < mHadronRapidityMax) ) countH[ipdg]++;
127+
if(isOkAtPartonicLevel && countH[ipdg] >= mHadronMultiplicity) return true;
128+
ipdg++;
129+
}
130+
}
131+
return false;
132+
};
133+
134+
135+
private:
136+
// Interface to override import particles
137+
Pythia8::Event mOutputEvent;
138+
139+
// Control gap-triggering
140+
unsigned long long mGeneratedEvents;
141+
int mInverseTriggerRatio;
142+
Pythia8::Pythia pythiaMBgen; // minimum bias event
143+
TString mConfigMBdecays;
144+
int mPDG;
145+
std::vector<int> mHadronsPDGs;
146+
int mHadronMultiplicity;
147+
double mRapidityMin;
148+
double mRapidityMax;
149+
double mHadronRapidityMin;
150+
double mHadronRapidityMax;
151+
bool mVerbose;
152+
};
153+
154+
}
155+
156+
}
157+
158+
FairGenerator*
159+
GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq(double rapidityMin = -4.3, double rapidityMax = -2.3, bool verbose = false, bool isbb = false, bool forceSemimuonicDecay = false)
160+
{
161+
TString pdgs;
162+
auto gen = new o2::eventgen::GeneratorEvtGen<o2::eventgen::GeneratorHFToMu_EvtGenFwdY_gaptriggered>();
163+
if (isbb == false) {
164+
gen->setPDG(4);
165+
pdgs = "411;421;431;4122;4232;4332";
166+
} else {
167+
gen->setPDG(5);
168+
pdgs = "511;521;531;541;5112;5122;5232;5132;5332";
169+
}
170+
gen->setRapidity(rapidityMin, rapidityMax);
171+
172+
gen->setRapidityHadron(rapidityMin,rapidityMax);
173+
gen->setHadronMultiplicity(1);
174+
TString pathO2table = gSystem->ExpandPathName("$O2DPG_MC_CONFIG_ROOT/MC/config/PWGDQ/pythia8/decayer/switchOffChadrons.cfg");
175+
gen->readFile(pathO2table.Data());
176+
gen->setConfigMBdecays(pathO2table);
177+
gen->setVerbose(verbose);
178+
179+
std::string spdg;
180+
TObjArray* obj = pdgs.Tokenize(";");
181+
gen->SetSizePdg(obj->GetEntriesFast());
182+
for (int i = 0; i < obj->GetEntriesFast(); i++) {
183+
spdg = obj->At(i)->GetName();
184+
gen->AddPdg(std::stoi(spdg), i);
185+
printf("PDG %d \n", std::stoi(spdg));
186+
gen->addHadronPDGs(std::stoi(spdg));
187+
}
188+
if (forceSemimuonicDecay == 1) gen->SetForceDecay(kEvtSemiMuonic);
189+
else gen->SetForceDecay(kEvtAll);
190+
191+
// set random seed
192+
gen->readString("Random:setSeed on");
193+
uint random_seed;
194+
unsigned long long int random_value = 0;
195+
ifstream urandom("/dev/urandom", ios::in|ios::binary);
196+
urandom.read(reinterpret_cast<char*>(&random_value), sizeof(random_seed));
197+
gen->readString(Form("Random:seed = %d", random_value % 900000001));
198+
199+
// print debug
200+
//gengen->PrintDebug();
201+
202+
return gen;
203+
}

‎MC/config/PWGDQ/ini/GeneratorHF_bbbarToDDbarToMuons_fwdy_TriggerGap.ini‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
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_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq.C
6+
funcName = GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq(-4.3, -2.3, false, true, false)
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_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_hf.cfg
13+
hooksFileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
14+
hooksFuncName = pythia8_userhooks_bbbar(-4.3, -2.3)

‎MC/config/PWGDQ/ini/GeneratorHF_ccbarToMuonsSemileptonic_fwdy_inel_TriggerGap.ini‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
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_MC_CONFIG_ROOT}/MC/config/PWGDQ/external/generator/GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq.C
6+
funcName = GeneratorHFToMu_EvtGenFwdY_gaptriggered_dq(-4.3, -2.3, false, false, true)
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_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_inel.cfg
13+
hooksFileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/hooks/pythia8_userhooks_qqbar.C
14+
hooksFuncName = pythia8_userhooks_ccbar(-4.3, -2.3)

‎MC/config/PWGDQ/ini/tests/GeneratorHF_bbbarToDDbarToMuons_fwdy_TriggerGap.C‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
int External()
2+
{
3+
int checkPdgDecay = 13;
4+
std::string path{"o2sim_Kine.root"};
5+
TFile file(path.c_str(), "READ");
6+
if (file.IsZombie()) {
7+
std::cerr << "Cannot open ROOT file " << path << "\n";
8+
return 1;
9+
}
10+
auto tree = (TTree*)file.Get("o2sim");
11+
std::vector<o2::MCTrack>* tracks{};
12+
tree->SetBranchAddress("MCTrack", &tracks);
13+
14+
int nLeptons{};
15+
int nLeptonsInAcceptance{};
16+
int nLeptonsToBeDone{};
17+
int nSignalPairs{};
18+
int nLeptonPairs{};
19+
int nLeptonPairsInAcceptance{};
20+
int nLeptonPairsToBeDone{};
21+
auto nEvents = tree->GetEntries();
22+
23+
for (int i = 0; i < nEvents; i++) {
24+
tree->GetEntry(i);
25+
int nleptonseinacc = 0;
26+
int nleptonse = 0;
27+
int nleptonseToBeDone = 0;
28+
int nopenHeavy = 0;
29+
for (auto& track : *tracks) {
30+
auto pdg = track.GetPdgCode();
31+
auto y = track.GetRapidity();
32+
if (std::abs(pdg) == checkPdgDecay) {
33+
int igmother = track.getMotherTrackId();
34+
if (igmother > 0) {
35+
auto gmTrack = (*tracks)[igmother];
36+
int gmpdg = gmTrack.GetPdgCode();
37+
if ( int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 ) {
38+
nLeptons++;
39+
nleptonse++;
40+
if (-4.3 < y && y < -2.2) {
41+
nleptonseinacc++;
42+
nLeptonsInAcceptance++;
43+
}
44+
if (track.getToBeDone()) {
45+
nLeptonsToBeDone++;
46+
nleptonseToBeDone++;
47+
}
48+
}
49+
}
50+
} else if (std::abs(pdg) == 411 || std::abs(pdg) == 421 || std::abs(pdg) == 431 || std::abs(pdg) == 4122 || std::abs(pdg) == 4132 || std::abs(pdg) == 4232 || std::abs(pdg) == 4332 || std::abs(pdg) == 511 || std::abs(pdg) == 521 || std::abs(pdg) == 531 || std::abs(pdg) == 541 || std::abs(pdg) == 5112 || std::abs(pdg) == 5122 || std::abs(pdg) == 5232 || std::abs(pdg) == 5132 || std::abs(pdg) == 5332) {
51+
nopenHeavy++;
52+
}
53+
}
54+
if (nopenHeavy > 1) nSignalPairs++;
55+
if (nleptonse > 1) nLeptonPairs++;
56+
if (nleptonseToBeDone > 1) nLeptonPairsToBeDone++;
57+
if (nleptonseinacc > 1) nLeptonPairsInAcceptance++;
58+
}
59+
std::cout << "#events: " << nEvents << "\n"
60+
<< "#muons in acceptance: " << nLeptonsInAcceptance << "\n"
61+
<< "#muon pairs in acceptance: " << nLeptonPairsInAcceptance << "\n"
62+
<< "#muons: " << nLeptons << "\n"
63+
<< "#muons to be done: " << nLeptonsToBeDone << "\n"
64+
<< "#signal pairs: " << nSignalPairs << "\n"
65+
<< "#muon pairs: " << nLeptonPairs << "\n"
66+
<< "#muon pairs to be done: " << nLeptonPairsToBeDone << "\n";
67+
if (nLeptonPairs != nLeptonPairsToBeDone) {
68+
std::cerr << "The number of muon pairs should be the same as the number of muon pairs which should be transported.\n";
69+
return 1;
70+
}
71+
if (nLeptons != nLeptonsToBeDone) {
72+
std::cerr << "The number of muons should be the same as the number of muons which should be transported.\n";
73+
return 1;
74+
}
75+
76+
return 0;
77+
}

‎MC/config/PWGDQ/ini/tests/GeneratorHF_ccbarToMuonsSemileptonic_fwdy_inel_TriggerGap.C‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
int External()
2+
{
3+
int checkPdgDecay = 13;
4+
std::string path{"o2sim_Kine.root"};
5+
TFile file(path.c_str(), "READ");
6+
if (file.IsZombie()) {
7+
std::cerr << "Cannot open ROOT file " << path << "\n";
8+
return 1;
9+
}
10+
auto tree = (TTree*)file.Get("o2sim");
11+
std::vector<o2::MCTrack>* tracks{};
12+
tree->SetBranchAddress("MCTrack", &tracks);
13+
14+
int nLeptons{};
15+
int nLeptonsInAcceptance{};
16+
int nLeptonsToBeDone{};
17+
int nSignalPairs{};
18+
int nLeptonPairs{};
19+
int nLeptonPairsInAcceptance{};
20+
int nLeptonPairsToBeDone{};
21+
auto nEvents = tree->GetEntries();
22+
23+
for (int i = 0; i < nEvents; i++) {
24+
tree->GetEntry(i);
25+
int nleptonseinacc = 0;
26+
int nleptonse = 0;
27+
int nleptonseToBeDone = 0;
28+
int nopenHeavy = 0;
29+
for (auto& track : *tracks) {
30+
auto pdg = track.GetPdgCode();
31+
auto y = track.GetRapidity();
32+
if (std::abs(pdg) == checkPdgDecay) {
33+
int igmother = track.getMotherTrackId();
34+
if (igmother > 0) {
35+
auto gmTrack = (*tracks)[igmother];
36+
int gmpdg = gmTrack.GetPdgCode();
37+
if ( int(std::abs(gmpdg)/100.) == 4 || int(std::abs(gmpdg)/1000.) == 4 || int(std::abs(gmpdg)/100.) == 5 || int(std::abs(gmpdg)/1000.) == 5 ) {
38+
nLeptons++;
39+
nleptonse++;
40+
if (-4.3 < y && y < -2.2) {
41+
nleptonseinacc++;
42+
nLeptonsInAcceptance++;
43+
}
44+
if (track.getToBeDone()) {
45+
nLeptonsToBeDone++;
46+
nleptonseToBeDone++;
47+
}
48+
}
49+
}
50+
} else if (std::abs(pdg) == 411 || std::abs(pdg) == 421 || std::abs(pdg) == 431 || std::abs(pdg) == 4122 || std::abs(pdg) == 4132 || std::abs(pdg) == 4232 || std::abs(pdg) == 4332) {
51+
nopenHeavy++;
52+
}
53+
}
54+
if (nopenHeavy > 1) nSignalPairs++;
55+
if (nleptonse > 1) nLeptonPairs++;
56+
if (nleptonseToBeDone > 1) nLeptonPairsToBeDone++;
57+
if (nleptonseinacc > 1) nLeptonPairsInAcceptance++;
58+
}
59+
std::cout << "#events: " << nEvents << "\n"
60+
<< "#muons in acceptance: " << nLeptonsInAcceptance << "\n"
61+
<< "#muon pairs in acceptance: " << nLeptonPairsInAcceptance << "\n"
62+
<< "#muons: " << nLeptons << "\n"
63+
<< "#muons to be done: " << nLeptonsToBeDone << "\n"
64+
<< "#signal pairs: " << nSignalPairs << "\n"
65+
<< "#muon pairs: " << nLeptonPairs << "\n"
66+
<< "#muon pairs to be done: " << nLeptonPairsToBeDone << "\n";
67+
if (nLeptonPairs != nLeptonPairsToBeDone) {
68+
std::cerr << "The number of muon pairs should be the same as the number of muon pairs which should be transported.\n";
69+
return 1;
70+
}
71+
if (nLeptons != nLeptonsToBeDone) {
72+
std::cerr << "The number of muons should be the same as the number of muons which should be transported.\n";
73+
return 1;
74+
}
75+
76+
return 0;
77+
}

0 commit comments

Comments
 (0)