Skip to content

Commit 90aa3cd

Browse files
BongHwialcaliva
BongHwi
authored and
alcaliva
committed
Add the exotic injection script for resonances
(cherry picked from commit 5c07c41)
1 parent 81dcf48 commit 90aa3cd

File tree

6 files changed

+580
-0
lines changed

6 files changed

+580
-0
lines changed

‎MC/config/PWGLF/ini/GeneratorLF_Resonances_PbPb_exotic.ini‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
[GeneratorExternal]
2+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_LF.C
3+
funcName=generateLF("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/resonancelistgun_exotic_pbpb.json", true, 4)
4+
5+
# [GeneratorPythia8]
6+
# config=${O2_ROOT}/share/Generators/egconfig/pythia8_hi.cfg
7+
8+
[DecayerPythia8] # after for transport code!
9+
config[0]=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
10+
config[1]=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/resonances.cfg

‎MC/config/PWGLF/ini/GeneratorLF_Resonances_pp_exotic.ini‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
[GeneratorExternal]
2+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_LF.C
3+
funcName=generateLF("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/resonancelistgun_exotic.json", true, 4)
4+
5+
# [GeneratorPythia8] # if triggered then this will be used as the background event
6+
# config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg
7+
8+
[DecayerPythia8] # after for transport code!
9+
config[0]=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
10+
config[1]=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/resonances.cfg

‎MC/config/PWGLF/ini/tests/GeneratorLF_Resonances_PbPb_exotic.C‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
int External()
2+
{
3+
std::string path{"o2sim_Kine.root"};
4+
int numberOfInjectedSignalsPerEvent{1};
5+
int numberOfGapEvents{4};
6+
int numberOfEventsProcessed{0};
7+
int numberOfEventsProcessedWithoutInjection{0};
8+
std::vector<int> injectedPDGs = {
9+
9010221, // f_0(980)
10+
10221, // f_0(1370)
11+
9030221, // f_0(1500)
12+
10331, // f_0(1710)
13+
20223, // f_1(1285)
14+
20333, // f_1(1420)
15+
335, // f_2(1525)
16+
10323, // K1(1270)+
17+
-10323, // K1(1270)-bar
18+
123314, // Xi(1820)-
19+
-123314, // Xi(1820)+
20+
123324, // Xi(1820)0
21+
-123324 // Xi(1820)0bar
22+
};
23+
std::vector<std::vector<int>> decayDaughters = {
24+
{211, -211}, // f_0(980)
25+
{310, 310}, // f_0(1370)
26+
{310, 310}, // f_0(1500)
27+
{310, 310}, // f_0(1710)
28+
{310, -321, 211}, // f_1(1285)
29+
{310, -321, 211}, // f_1(1420)
30+
{310, 310}, // f_2(1525)
31+
{321, 211}, // K1(1270)+
32+
{-321, -211}, // K1(1270)-bar
33+
{2212, 211}, // Delta(1232)+
34+
{3122, -311}, // Xi(1820)-
35+
{3122, 311}, // Xi(1820)+
36+
{3122, 310}, // Xi(1820)0
37+
{-3122, 310} // Xi(1820)0bar
38+
};
39+
40+
auto nInjection = injectedPDGs.size();
41+
42+
TFile file(path.c_str(), "READ");
43+
if (file.IsZombie())
44+
{
45+
std::cerr << "Cannot open ROOT file " << path << "\n";
46+
return 1;
47+
}
48+
49+
auto tree = (TTree *)file.Get("o2sim");
50+
if (!tree)
51+
{
52+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
53+
return 1;
54+
}
55+
std::vector<o2::MCTrack> *tracks{};
56+
tree->SetBranchAddress("MCTrack", &tracks);
57+
58+
std::vector<int> nSignal;
59+
for (int i = 0; i < nInjection; i++)
60+
{
61+
nSignal.push_back(0);
62+
}
63+
std::vector<std::vector<int>> nDecays;
64+
std::vector<int> nNotDecayed;
65+
for (int i = 0; i < nInjection; i++)
66+
{
67+
std::vector<int> nDecay;
68+
for (int j = 0; j < decayDaughters[i].size(); j++)
69+
{
70+
nDecay.push_back(0);
71+
}
72+
nDecays.push_back(nDecay);
73+
nNotDecayed.push_back(0);
74+
}
75+
auto nEvents = tree->GetEntries();
76+
bool hasInjection = false;
77+
for (int i = 0; i < nEvents; i++)
78+
{
79+
hasInjection = false;
80+
numberOfEventsProcessed++;
81+
auto check = tree->GetEntry(i);
82+
for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack)
83+
{
84+
auto track = tracks->at(idxMCTrack);
85+
auto pdg = track.GetPdgCode();
86+
auto it = std::find(injectedPDGs.begin(), injectedPDGs.end(), pdg);
87+
int index = std::distance(injectedPDGs.begin(), it); // index of injected PDG
88+
if (it != injectedPDGs.end()) // found
89+
{
90+
// count signal PDG
91+
nSignal[index]++;
92+
if (track.getFirstDaughterTrackId() < 0)
93+
{
94+
nNotDecayed[index]++;
95+
continue;
96+
}
97+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j)
98+
{
99+
auto pdgDau = tracks->at(j).GetPdgCode();
100+
bool foundDau = false;
101+
// count decay PDGs
102+
for (int idxDaughter = 0; idxDaughter < decayDaughters[index].size(); ++idxDaughter)
103+
{
104+
if (pdgDau == decayDaughters[index][idxDaughter])
105+
{
106+
nDecays[index][idxDaughter]++;
107+
foundDau = true;
108+
hasInjection = true;
109+
break;
110+
}
111+
}
112+
if (!foundDau)
113+
{
114+
std::cerr << "Decay daughter not found: " << pdg << " -> " << pdgDau << "\n";
115+
}
116+
}
117+
}
118+
}
119+
if (!hasInjection)
120+
{
121+
numberOfEventsProcessedWithoutInjection++;
122+
}
123+
}
124+
std::cout << "--------------------------------\n";
125+
std::cout << "# Events: " << nEvents << "\n";
126+
for (int i = 0; i < nInjection; i++)
127+
{
128+
std::cout << "# Mother \n";
129+
std::cout << injectedPDGs[i] << " generated: " << nSignal[i] << ", " << nNotDecayed[i] << " did not decay\n";
130+
if (nSignal[i] == 0)
131+
{
132+
std::cerr << "No generated: " << injectedPDGs[i] << "\n";
133+
// return 1; // At least one of the injected particles should be generated
134+
}
135+
for (int j = 0; j < decayDaughters[i].size(); j++)
136+
{
137+
std::cout << "# Daughter " << decayDaughters[i][j] << ": " << nDecays[i][j] << "\n";
138+
}
139+
// if (nSignal[i] != nEvents * numberOfInjectedSignalsPerEvent)
140+
// {
141+
// std::cerr << "Number of generated: " << injectedPDGs[i] << ", lower than expected\n";
142+
// // return 1; // Don't need to return 1, since the number of generated particles is not the same for each event
143+
// }
144+
}
145+
std::cout << "--------------------------------\n";
146+
std::cout << "Number of events processed: " << numberOfEventsProcessed << "\n";
147+
std::cout << "Number of input for the gap events: " << numberOfGapEvents << "\n";
148+
std::cout << "Number of events processed without injection: " << numberOfEventsProcessedWithoutInjection << "\n";
149+
// injected event + numberOfGapEvents*gap events + injected event + numberOfGapEvents*gap events + ...
150+
// total fraction of the gap event: numberOfEventsProcessedWithoutInjection/numberOfEventsProcessed
151+
float ratioOfNormalEvents = numberOfEventsProcessedWithoutInjection / numberOfEventsProcessed;
152+
if (ratioOfNormalEvents > 0.75)
153+
{
154+
std::cout << "The number of injected event is loo low!!" << std::endl;
155+
return 1;
156+
}
157+
158+
return 0;
159+
}
160+
161+
void GeneratorLF_Resonances_PbPb5360_injection() { External(); }

‎MC/config/PWGLF/ini/tests/GeneratorLF_Resonances_pp_exotic.C‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
1+
int External()
2+
{
3+
std::string path{"o2sim_Kine.root"};
4+
int numberOfInjectedSignalsPerEvent{1};
5+
int numberOfGapEvents{4};
6+
int numberOfEventsProcessed{0};
7+
int numberOfEventsProcessedWithoutInjection{0};
8+
std::vector<int> injectedPDGs = {
9+
9010221, // f_0(980)
10+
10221, // f_0(1370)
11+
9030221, // f_0(1500)
12+
10331, // f_0(1710)
13+
20223, // f_1(1285)
14+
20333, // f_1(1420)
15+
335, // f_2(1525)
16+
10323, // K1(1270)+
17+
-10323, // K1(1270)-bar
18+
123314, // Xi(1820)-
19+
-123314, // Xi(1820)+
20+
123324, // Xi(1820)0
21+
-123324 // Xi(1820)0bar
22+
};
23+
std::vector<std::vector<int>> decayDaughters = {
24+
{211, -211}, // f_0(980)
25+
{310, 310}, // f_0(1370)
26+
{310, 310}, // f_0(1500)
27+
{310, 310}, // f_0(1710)
28+
{310, -321, 211}, // f_1(1285)
29+
{310, -321, 211}, // f_1(1420)
30+
{310, 310}, // f_2(1525)
31+
{321, 211}, // K1(1270)+
32+
{-321, -211}, // K1(1270)-bar
33+
{2212, 211}, // Delta(1232)+
34+
{3122, -311}, // Xi(1820)-
35+
{3122, 311}, // Xi(1820)+
36+
{3122, 310}, // Xi(1820)0
37+
{-3122, 310} // Xi(1820)0bar
38+
};
39+
40+
auto nInjection = injectedPDGs.size();
41+
42+
TFile file(path.c_str(), "READ");
43+
if (file.IsZombie())
44+
{
45+
std::cerr << "Cannot open ROOT file " << path << "\n";
46+
return 1;
47+
}
48+
49+
auto tree = (TTree *)file.Get("o2sim");
50+
if (!tree)
51+
{
52+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
53+
return 1;
54+
}
55+
std::vector<o2::MCTrack> *tracks{};
56+
tree->SetBranchAddress("MCTrack", &tracks);
57+
58+
std::vector<int> nSignal;
59+
for (int i = 0; i < nInjection; i++)
60+
{
61+
nSignal.push_back(0);
62+
}
63+
std::vector<std::vector<int>> nDecays;
64+
std::vector<int> nNotDecayed;
65+
for (int i = 0; i < nInjection; i++)
66+
{
67+
std::vector<int> nDecay;
68+
for (int j = 0; j < decayDaughters[i].size(); j++)
69+
{
70+
nDecay.push_back(0);
71+
}
72+
nDecays.push_back(nDecay);
73+
nNotDecayed.push_back(0);
74+
}
75+
auto nEvents = tree->GetEntries();
76+
bool hasInjection = false;
77+
for (int i = 0; i < nEvents; i++)
78+
{
79+
hasInjection = false;
80+
numberOfEventsProcessed++;
81+
auto check = tree->GetEntry(i);
82+
for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack)
83+
{
84+
auto track = tracks->at(idxMCTrack);
85+
auto pdg = track.GetPdgCode();
86+
auto it = std::find(injectedPDGs.begin(), injectedPDGs.end(), pdg);
87+
int index = std::distance(injectedPDGs.begin(), it); // index of injected PDG
88+
if (it != injectedPDGs.end()) // found
89+
{
90+
// count signal PDG
91+
nSignal[index]++;
92+
if (track.getFirstDaughterTrackId() < 0)
93+
{
94+
nNotDecayed[index]++;
95+
continue;
96+
}
97+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j)
98+
{
99+
auto pdgDau = tracks->at(j).GetPdgCode();
100+
bool foundDau = false;
101+
// count decay PDGs
102+
for (int idxDaughter = 0; idxDaughter < decayDaughters[index].size(); ++idxDaughter)
103+
{
104+
if (pdgDau == decayDaughters[index][idxDaughter])
105+
{
106+
nDecays[index][idxDaughter]++;
107+
foundDau = true;
108+
hasInjection = true;
109+
break;
110+
}
111+
}
112+
if (!foundDau)
113+
{
114+
std::cerr << "Decay daughter not found: " << pdg << " -> " << pdgDau << "\n";
115+
}
116+
}
117+
}
118+
}
119+
if (!hasInjection)
120+
{
121+
numberOfEventsProcessedWithoutInjection++;
122+
}
123+
}
124+
std::cout << "--------------------------------\n";
125+
std::cout << "# Events: " << nEvents << "\n";
126+
for (int i = 0; i < nInjection; i++)
127+
{
128+
std::cout << "# Mother \n";
129+
std::cout << injectedPDGs[i] << " generated: " << nSignal[i] << ", " << nNotDecayed[i] << " did not decay\n";
130+
if (nSignal[i] == 0)
131+
{
132+
std::cerr << "No generated: " << injectedPDGs[i] << "\n";
133+
// return 1; // At least one of the injected particles should be generated
134+
}
135+
for (int j = 0; j < decayDaughters[i].size(); j++)
136+
{
137+
std::cout << "# Daughter " << decayDaughters[i][j] << ": " << nDecays[i][j] << "\n";
138+
}
139+
// if (nSignal[i] != nEvents * numberOfInjectedSignalsPerEvent)
140+
// {
141+
// std::cerr << "Number of generated: " << injectedPDGs[i] << ", lower than expected\n";
142+
// // return 1; // Don't need to return 1, since the number of generated particles is not the same for each event
143+
// }
144+
}
145+
std::cout << "--------------------------------\n";
146+
std::cout << "Number of events processed: " << numberOfEventsProcessed << "\n";
147+
std::cout << "Number of input for the gap events: " << numberOfGapEvents << "\n";
148+
std::cout << "Number of events processed without injection: " << numberOfEventsProcessedWithoutInjection << "\n";
149+
// injected event + numberOfGapEvents*gap events + injected event + numberOfGapEvents*gap events + ...
150+
// total fraction of the gap event: numberOfEventsProcessedWithoutInjection/numberOfEventsProcessed
151+
float ratioOfNormalEvents = numberOfEventsProcessedWithoutInjection / numberOfEventsProcessed;
152+
if (ratioOfNormalEvents > 0.75)
153+
{
154+
std::cout << "The number of injected event is loo low!!" << std::endl;
155+
return 1;
156+
}
157+
158+
return 0;
159+
}
160+
161+
void GeneratorLF_Resonances_pp1360_injection() { External(); }

0 commit comments

Comments
 (0)