Skip to content

Commit 5329bae

Browse files
motomioyamotomioya
authored
Add scripts for LF to mumu (#1707)
* Add bbtomuons without forcing semileptonic decay * Add test file for bb->mumu with natural decay * Add LF to mumu * Change gap
1 parent dc99792 commit 5329bae

File tree

4 files changed

+410
-0
lines changed

4 files changed

+410
-0
lines changed

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

Lines changed: 286 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,286 @@
1+
R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/EvtGen)
2+
R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGDQ/external/generator)
3+
R__ADD_INCLUDE_PATH($O2DPG_ROOT/MC/config/PWGEM/external/generator)
4+
R__LOAD_LIBRARY(libpythia6)
5+
R__LOAD_LIBRARY(libEGPythia6)
6+
#include "GeneratorEvtGen.C"
7+
#include "GeneratorCocktail.C"
8+
#include "Generators/GeneratorPythia8.h"
9+
#include "Pythia8/Pythia.h"
10+
11+
using namespace std;
12+
using namespace Pythia8;
13+
14+
namespace o2 {
15+
namespace eventgen {
16+
17+
class CocktailParam : public GeneratorTGenerator {
18+
public:
19+
CocktailParam(GeneratorParam *thisGenerator)
20+
: GeneratorTGenerator(thisGenerator->GetName()) {
21+
setTGenerator(thisGenerator);
22+
};
23+
24+
~CocktailParam() { delete thisGenerator; };
25+
26+
private:
27+
GeneratorParam *thisGenerator = nullptr;
28+
};
29+
30+
31+
//my generator class
32+
class GeneratorPythia8GapTriggeredLFmumu : public GeneratorPythia8 {
33+
34+
public:
35+
GeneratorPythia8GapTriggeredLFmumu() : GeneratorPythia8() {
36+
mGeneratedEvents = 0;
37+
mInverseTriggerRatio = 1;
38+
fGeneratorCocktail = 0x0;
39+
mMode = -1;
40+
mTargetPDG = 0;
41+
};
42+
43+
GeneratorPythia8GapTriggeredLFmumu(int lInputTriggerRatio, float yMin, float yMax, int nPart, int mode) : GeneratorPythia8() {
44+
mGeneratedEvents = 0;
45+
mInverseTriggerRatio = lInputTriggerRatio;
46+
mMode = mode;
47+
// LMee cocktail settings:
48+
float minPt = 0;
49+
float maxPt = 25;
50+
float phiMin = 0.;
51+
float phiMax = 360.;
52+
Weighting_t weightMode = kNonAnalog;
53+
54+
//create cocktail generator : eta, eta', rho, omega, phi
55+
fGeneratorCocktail = new o2::eventgen::GeneratorEvtGen<GeneratorCocktail>();
56+
57+
// EXODUS decayer
58+
TString O2DPG_ROOT = TString(getenv("O2DPG_ROOT"));
59+
auto decayer = new PythiaDecayerConfig();
60+
decayer->SetDecayerExodus();
61+
TString useLMeeDecaytable = "$O2DPG_ROOT/MC/config/PWGEM/decaytables/decaytable_LMee.dat";
62+
useLMeeDecaytable=useLMeeDecaytable.ReplaceAll("$O2DPG_ROOT",O2DPG_ROOT);
63+
useLMeeDecaytable=useLMeeDecaytable.ReplaceAll("${O2DPG_ROOT}",O2DPG_ROOT);
64+
decayer->SetDecayTableFile(useLMeeDecaytable.Data());
65+
decayer->ReadDecayTable();
66+
67+
//Param
68+
GeneratorParamEMlib *emlib = new GeneratorParamEMlib();
69+
70+
// eta
71+
auto geneta = new GeneratorParam(nPart,emlib,GeneratorParamEMlib::kEta,"eta"); // 221
72+
geneta->SetName("eta");
73+
geneta->SetMomentumRange(0., 25.);
74+
geneta->SetPtRange(minPt, maxPt);
75+
geneta->SetYRange(yMin, yMax);
76+
geneta->SetPhiRange(phiMin, phiMax);
77+
geneta->SetWeighting(weightMode); // flat pt, y and v2 zero
78+
geneta->SetDecayer(decayer); // EXOUS;
79+
geneta->SetForceDecay(kDiMuon); // Dielectrons
80+
geneta->SetForceGammaConversion(kFALSE);
81+
geneta->SetSelectAll(kTRUE); // Store also the gamma in pi0->e+e-gamma
82+
geneta->Init();
83+
CocktailParam *newgeneta = new CocktailParam(geneta);
84+
85+
// etaprime
86+
auto genetaprime = new GeneratorParam(nPart,emlib,GeneratorParamEMlib::kEtaprime,"etaprime"); // 331
87+
genetaprime->SetName("etaprime");
88+
genetaprime->SetMomentumRange(0., 25.);
89+
genetaprime->SetPtRange(minPt, maxPt);
90+
genetaprime->SetYRange(yMin, yMax);
91+
genetaprime->SetPhiRange(phiMin, phiMax);
92+
genetaprime->SetWeighting(weightMode); // flat pt, y and v2 zero
93+
genetaprime->SetDecayer(decayer); // EXOUS;
94+
genetaprime->SetForceDecay(kDiMuon); // Dielectrons
95+
genetaprime->SetForceGammaConversion(kFALSE);
96+
genetaprime->SetSelectAll(kTRUE); // Store also the gamma in pi0->e+e-gamma
97+
genetaprime->Init();
98+
CocktailParam *newgenetaprime = new CocktailParam(genetaprime);
99+
100+
// rho
101+
auto genrho = new GeneratorParam(nPart,emlib,GeneratorParamEMlib::kRho0,"rho"); // 113
102+
genrho->SetName("rho");
103+
genrho->SetMomentumRange(0., 25.);
104+
genrho->SetPtRange(minPt, maxPt);
105+
genrho->SetYRange(yMin, yMax);
106+
genrho->SetPhiRange(phiMin, phiMax);
107+
genrho->SetWeighting(weightMode); // flat pt, y and v2 zero
108+
genrho->SetDecayer(decayer); // EXOUS;
109+
genrho->SetForceDecay(kDiMuon); // Dielectrons
110+
genrho->SetForceGammaConversion(kFALSE);
111+
genrho->SetSelectAll(kTRUE); // Store also the gamma in pi0->e+e-gamma
112+
genrho->Init();
113+
CocktailParam *newgenrho = new CocktailParam(genrho);
114+
115+
// Omega
116+
auto genomega = new GeneratorParam(nPart,emlib,GeneratorParamEMlib::kOmega,"omega"); //223
117+
genomega->SetName("omega");
118+
genomega->SetMomentumRange(0., 25.);
119+
genomega->SetPtRange(minPt, maxPt);
120+
genomega->SetYRange(yMin, yMax);
121+
genomega->SetPhiRange(phiMin, phiMax);
122+
genomega->SetWeighting(weightMode); // flat pt, y and v2 zero
123+
genomega->SetDecayer(decayer); // EXOUS;
124+
genomega->SetForceDecay(kDiMuon); // Dielectrons
125+
genomega->SetForceGammaConversion(kFALSE);
126+
genomega->SetSelectAll(kTRUE); // Store also the gamma in pi0->e+e-gamma
127+
genomega->Init();
128+
CocktailParam *newgenomega = new CocktailParam(genomega);
129+
130+
// phi
131+
auto genphi = new GeneratorParam(nPart,emlib,GeneratorParamEMlib::kPhi,"phi"); //333
132+
genphi->SetName("phi");
133+
genphi->SetMomentumRange(0., 25.);
134+
genphi->SetPtRange(minPt, maxPt);
135+
genphi->SetYRange(yMin, yMax);
136+
genphi->SetPhiRange(phiMin, phiMax);
137+
genphi->SetWeighting(weightMode); // flat pt, y and v2 zero
138+
genphi->SetDecayer(decayer); // EXOUS;
139+
genphi->SetForceDecay(kDiMuon); // Dielectrons
140+
genphi->SetForceGammaConversion(kFALSE);
141+
genphi->SetSelectAll(kTRUE); // Store also the gamma in pi0->e+e-gamma
142+
genphi->Init();
143+
CocktailParam *newgenphi = new CocktailParam(genphi);
144+
145+
int target_pdg = 1;
146+
147+
if (mMode < 0) {
148+
target_pdg = 1;
149+
cout << "all-particle mode is selected. all 5 mesons are injected in each event" << endl;
150+
cout << "add eta for signal" << endl;
151+
fGeneratorCocktail->AddGenerator(newgeneta, 1);
152+
cout << "add etaprime for signal" << endl;
153+
fGeneratorCocktail->AddGenerator(newgenetaprime, 1);
154+
cout << "add rho for signal" << endl;
155+
fGeneratorCocktail->AddGenerator(newgenrho, 1);
156+
cout << "add omega for signal" << endl;
157+
fGeneratorCocktail->AddGenerator(newgenomega, 1);
158+
cout << "add phi for signal" << endl;
159+
fGeneratorCocktail->AddGenerator(newgenphi, 1);
160+
} else if (mMode < 100) {
161+
cout << "1-particle Mode is selected. 1 meson selected randomly per job is injected in each event" << endl;
162+
TRandom3 *r3 = new TRandom3(0);
163+
double rndm = r3->Rndm();
164+
printf("rndm = %f\n", rndm);
165+
166+
if (rndm < 1/5.) {
167+
cout << "add eta for signal" << endl;
168+
fGeneratorCocktail->AddGenerator(newgeneta, 1);
169+
target_pdg = 221;
170+
} else if (rndm < 2/5.) {
171+
cout << "add etaprime for signal" << endl;
172+
fGeneratorCocktail->AddGenerator(newgenetaprime, 1);
173+
target_pdg = 331;
174+
} else if (rndm < 3/5.) {
175+
cout << "add rho for signal" << endl;
176+
fGeneratorCocktail->AddGenerator(newgenrho, 1);
177+
target_pdg = 113;
178+
} else if (rndm < 4/5.) {
179+
cout << "add omega for signal" << endl;
180+
fGeneratorCocktail->AddGenerator(newgenomega, 1);
181+
target_pdg = 223;
182+
} else if (rndm < 5/5.) {
183+
cout << "add phi for signal" << endl;
184+
fGeneratorCocktail->AddGenerator(newgenphi, 1);
185+
target_pdg = 333;
186+
}
187+
delete r3;
188+
} else { //directly select meson pdg
189+
target_pdg = mMode;
190+
switch (mMode) {
191+
case 221 :
192+
cout << "add eta for signal" << endl;
193+
fGeneratorCocktail->AddGenerator(newgeneta, 1);
194+
break;
195+
case 331 :
196+
cout << "add etaprime for signal" << endl;
197+
fGeneratorCocktail->AddGenerator(newgenetaprime, 1);
198+
break;
199+
case 113 :
200+
cout << "add rho for signal" << endl;
201+
fGeneratorCocktail->AddGenerator(newgenrho, 1);
202+
break;
203+
case 223 :
204+
cout << "add omega for signal" << endl;
205+
fGeneratorCocktail->AddGenerator(newgenomega, 1);
206+
break;
207+
case 333 :
208+
cout << "add phi for signal" << endl;
209+
fGeneratorCocktail->AddGenerator(newgenphi, 1);
210+
break;
211+
default:
212+
cout << "!WARNING! default : nothing is added to cocktail generator" << endl;
213+
target_pdg = 1;
214+
break;
215+
}
216+
}
217+
218+
// print debug
219+
fGeneratorCocktail->PrintDebug();
220+
fGeneratorCocktail->Init();
221+
222+
cout << "target_pdg for subGeneratorId is " << target_pdg << endl;
223+
addSubGenerator(0, "gap mb pythia");
224+
addSubGenerator(target_pdg, "injected cocktail");
225+
mTargetPDG = target_pdg;
226+
};
227+
228+
~GeneratorPythia8GapTriggeredLFmumu() = default;
229+
230+
protected:
231+
bool generateEvent() override
232+
{
233+
GeneratorPythia8::generateEvent();
234+
235+
if (mGeneratedEvents % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
236+
fGeneratorCocktail->generateEvent();
237+
notifySubGenerator(mTargetPDG);
238+
} else { // gap event
239+
notifySubGenerator(0);
240+
}
241+
mGeneratedEvents++;
242+
return true;
243+
}
244+
245+
bool importParticles() override
246+
{
247+
GeneratorPythia8::importParticles();
248+
249+
bool genOk = false;
250+
if ((mGeneratedEvents-1) % mInverseTriggerRatio == 0){ // add injected prompt signals to the stack
251+
fGeneratorCocktail->importParticles();
252+
int originalSize = mParticles.size();
253+
for(int ipart=0; ipart < fGeneratorCocktail->getParticles().size(); ipart++){
254+
TParticle part = TParticle(fGeneratorCocktail->getParticles().at(ipart));
255+
if(part.GetFirstMother() >= 0) part.SetFirstMother(part.GetFirstMother() + originalSize);
256+
if(part.GetFirstDaughter() >= 0) part.SetFirstDaughter(part.GetFirstDaughter() + originalSize);
257+
if(part.GetLastDaughter() >= 0) part.SetLastDaughter(part.GetLastDaughter() + originalSize);
258+
mParticles.push_back(part);
259+
// encodeParticleStatusAndTracking method already called in GeneratorEvtGen.C
260+
}
261+
fGeneratorCocktail->clearParticles();
262+
}
263+
264+
return true;
265+
}
266+
267+
private:
268+
GeneratorEvtGen<GeneratorCocktail> *fGeneratorCocktail;
269+
// Control gap-triggering
270+
unsigned long long mGeneratedEvents;
271+
int mInverseTriggerRatio;
272+
int mMode;
273+
int mTargetPDG;
274+
};
275+
276+
} // close eventgen
277+
} // close o2
278+
279+
// Predefined generators: // this function should be called in ini file.
280+
FairGenerator *GeneratorPythia8GapTriggeredLFmumu_ForEM(int inputTriggerRatio = 5, float yMin=-4.3, float yMax=-2.2, int nPart = 1, int mode = -1) {
281+
auto myGen = new GeneratorPythia8GapTriggeredLFmumu(inputTriggerRatio, yMin, yMax, nPart, mode);
282+
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
283+
myGen->readString("Random:setSeed on");
284+
myGen->readString("Random:seed " + std::to_string(seed));
285+
return myGen;
286+
}

‎MC/config/PWGEM/ini/Generator_GapTriggered_LFmumu.ini‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
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_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_LFmumu.C
6+
funcName = GeneratorPythia8GapTriggeredLFmumu_ForEM(5+1, -4.3, -2.2, 1, 1)
7+
8+
[GeneratorPythia8]
9+
config = ${O2DPG_ROOT}/MC/config/PWGEM/pythia8/generator/pythia8_MB_gapevent.cfg

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

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
int External()
2+
{
3+
std::string path{"o2sim_Kine.root"};
4+
TFile file(path.c_str(), "READ");
5+
if (file.IsZombie()) {
6+
std::cerr << "Cannot open ROOT file " << path << "\n";
7+
return 1;
8+
}
9+
10+
auto tree = (TTree*)file.Get("o2sim");
11+
std::vector<o2::MCTrack>* tracks{};
12+
tree->SetBranchAddress("MCTrack", &tracks);
13+
14+
int nMesons{};
15+
int nMesonsDiMuonDecay{};
16+
auto nEvents = tree->GetEntries();
17+
18+
for (int i = 0; i < nEvents; i++) {
19+
tree->GetEntry(i);
20+
for (auto& track : *tracks) {
21+
auto pdg = track.GetPdgCode();
22+
auto y = track.GetRapidity();
23+
if ((pdg == 221) || (pdg == 331) || (pdg == 223) || (pdg == 113) || (pdg == 333)) {
24+
if ((y>-4.3) && (y<-2.2)) {
25+
nMesons++;
26+
Int_t counterel = 0;
27+
Int_t counterpos = 0;
28+
int k1 = track.getFirstDaughterTrackId();
29+
int k2 = track.getLastDaughterTrackId();
30+
// k1 < k2 and no -1 for k2
31+
for (int d=k1; d <= k2; d++) {
32+
if (d>0) {
33+
auto decay = (*tracks)[d];
34+
int pdgdecay = decay.GetPdgCode();
35+
if (pdgdecay == 13) {
36+
counterel++;
37+
}
38+
if (pdgdecay == -13) {
39+
counterpos++;
40+
}
41+
}
42+
}
43+
if ((counterel>0) && (counterpos>0)) nMesonsDiMuonDecay++;
44+
}
45+
}
46+
}
47+
}
48+
49+
std::cout << "#events: " << nEvents << "\n"
50+
<< "#mesons: " << nMesons << "\n";
51+
52+
if (nMesons < (nEvents/5)) {
53+
std::cerr << "One should have at least one mesons in forward region per 5 events.\n";
54+
return 1;
55+
}
56+
//if (nMesonsDiMuonDecay < nEvents) {
57+
//std::cerr << "One meson to dimuon decay per event should be produced.\n";
58+
//return 1;
59+
//}
60+
61+
return 0;
62+
}

0 commit comments

Comments
 (0)