Skip to content

Commit 49bda2f

Browse files
stefanopolitanostefanopolitano
authored
[PWGHF] Update correlated bkg config and ini files in pp and PbPb (#1987)
* [PWHGF] Update correlated bkg config and ini files in pp and PbPb * Minor fix * Implementing Marcello comments * Implementing Marcello comments * Implementing Marcello comments * Implementing Mattia comments * Implementing Mattia comments * Adding test for corr. bkg
1 parent 05fde30 commit 49bda2f

7 files changed

+839
-159
lines changed

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.ini‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
### The external generator derives from GeneratorPythia8.
2+
[GeneratorExternal]
3+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_hf.C
4+
funcName=GeneratorPythia8EmbedHFCharmAndBeauty()
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_hardQCD_5TeV_CorrBkg.cfg
8+
includePartonEvent=true

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.ini‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,5 +4,5 @@ fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_py
44
funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5)
55

66
[GeneratorPythia8]
7-
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_CorrBkg.cfg
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_Mode2_CorrBkg.cfg
88
includePartonEvent=true

‎MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_PbPb_corrBkg.C‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
//std::string path{"tf1/sgn_1_Kine.root"};
4+
5+
int checkPdgQuarkOne{4};
6+
int checkPdgQuarkTwo{5};
7+
float ratioTrigger = 1.; // one event triggered out of 1
8+
9+
std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4232};
10+
std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
11+
{411, {
12+
{-321, 211, 211}, // K- π+ π+ (non-resonant)
13+
{-313, 321}, // K*0(892) K+
14+
{-10311, 321}, // K*0(1430) K+
15+
{211, 333}, // φ π+
16+
{-321, 321, 211}, // K- K+ π+ (non-resonant)
17+
{113, 211}, // ρ0 π+
18+
{225, 211}, // f2(1270) π+
19+
{-211, 211, 211} // π- π+ π+ (non-resonant)
20+
}},
21+
{421, {
22+
{-321, 211}, // K- π+ (non-resonant)
23+
{-321, 111, 211}, // K- π+ π0
24+
{213, -321}, // ρ+ K-
25+
{-313, 111}, // antiK*0(892) π0
26+
{-323, 211}, // K*-(892) π+
27+
{-211, 211}, // π- π+
28+
{213, -211}, // ρ+ π-
29+
{-211, 211, 111}, // π- π+ π0
30+
{-321, 321} // K- K+
31+
}},
32+
{431, {
33+
{211, 333}, // φ π+
34+
{-313, 321}, // antiK*(892) K+
35+
{333, 213}, // φ ρ
36+
{113, 211}, // ρ π+
37+
{225, 211}, // f2(1270) π+
38+
{-211, 211, 211}, // π- π+ π+ (s-wave)
39+
{313, 211}, // K*(892)0 π+
40+
{10221, 321}, // f0(1370) K+
41+
{113, 321}, // ρ0 K+
42+
{-211, 321, 211}, // π- K+ π+ (non-resonant)
43+
{221, 211} // η π+
44+
}},
45+
{4122, {
46+
{2212, -321, 211}, // p K- π+ (non-resonant)
47+
{2212, -313}, // p K*0(892)
48+
{2224, -321}, // Δ++ K-
49+
{102134, 211}, // Λ(1520) K-
50+
{2212, -321, 211, 111}, // p K- π+ π0
51+
{2212, -211, 211}, // p π- π+
52+
{2212, 333} // p φ
53+
}},
54+
{4232, {
55+
{-313, 2212}, // antiK*0(892) p
56+
{2212, -321, 211}, // p K- π+
57+
{2212, 333}, // p φ
58+
{3222, -211, 211} // Σ+ π- π+
59+
}},
60+
};
61+
62+
63+
TFile file(path.c_str(), "READ");
64+
if (file.IsZombie()) {
65+
std::cerr << "Cannot open ROOT file " << path << "\n";
66+
return 1;
67+
}
68+
69+
auto tree = (TTree *)file.Get("o2sim");
70+
std::vector<o2::MCTrack> *tracks{};
71+
tree->SetBranchAddress("MCTrack", &tracks);
72+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
73+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
74+
75+
int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
76+
int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
77+
auto nEvents = tree->GetEntries();
78+
79+
for (int i = 0; i < nEvents; i++) {
80+
tree->GetEntry(i);
81+
82+
// check subgenerator information
83+
//if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
84+
// bool isValid = false;
85+
// int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
86+
// if (subGeneratorId == 0) {
87+
// nEventsMB++;
88+
// } else if (subGeneratorId == checkPdgQuarkOne) {
89+
// nEventsInjOne++;
90+
// } else if (subGeneratorId == checkPdgQuarkTwo) {
91+
// nEventsInjTwo++;
92+
// }
93+
//}
94+
95+
for (auto &track : *tracks) {
96+
auto pdg = track.GetPdgCode();
97+
if (std::abs(pdg) == checkPdgQuarkOne) {
98+
nQuarksOne++;
99+
continue;
100+
}
101+
if (std::abs(pdg) == checkPdgQuarkTwo) {
102+
nQuarksTwo++;
103+
continue;
104+
}
105+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
106+
nSignals++; // count signal PDG
107+
108+
std::vector<int> pdgsDecay{};
109+
std::vector<int> pdgsDecayAntiPart{};
110+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
111+
auto pdgDau = tracks->at(j).GetPdgCode();
112+
pdgsDecay.push_back(pdgDau);
113+
if (pdgDau != 333) { // phi is antiparticle of itself
114+
pdgsDecayAntiPart.push_back(-pdgDau);
115+
} else {
116+
pdgsDecayAntiPart.push_back(pdgDau);
117+
}
118+
}
119+
120+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
121+
std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
122+
123+
for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
124+
if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
125+
nSignalGoodDecay++;
126+
break;
127+
}
128+
}
129+
}
130+
}
131+
}
132+
133+
std::cout << "--------------------------------\n";
134+
std::cout << "# Events: " << nEvents << "\n";
135+
//std::cout << "# MB events: " << nEventsMB << "\n";
136+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
137+
//std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
138+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
139+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
140+
std::cout <<"# signal hadrons: " << nSignals << "\n";
141+
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
142+
143+
//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
144+
// std::cerr << "Number of generated MB events different than expected\n";
145+
// return 1;
146+
//}
147+
//if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
148+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
149+
// return 1;
150+
//}
151+
//if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
152+
// std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
153+
// return 1;
154+
//}
155+
156+
if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
157+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
158+
return 1;
159+
}
160+
if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
161+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
162+
return 1;
163+
}
164+
165+
float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
166+
if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
167+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
168+
return 1;
169+
}
170+
171+
return 0;
172+
}

‎MC/config/PWGHF/ini/tests/GeneratorHF_D2H_ccbar_and_bbbar_gap5_Mode2_corrBkg.C‎

Lines changed: 170 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,170 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
4+
int checkPdgQuarkOne{4};
5+
int checkPdgQuarkTwo{5};
6+
float ratioTrigger = 1./5; // one event triggered out of 5
7+
8+
std::vector<int> checkPdgHadron{411, 421, 431, 4122, 4232};
9+
std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted pdg of daughters
10+
{411, {
11+
{-321, 211, 211}, // K- π+ π+ (non-resonant)
12+
{-313, 321}, // K*0(892) K+
13+
{-10311, 321}, // K*0(1430) K+
14+
{211, 333}, // φ π+
15+
{-321, 321, 211}, // K- K+ π+ (non-resonant)
16+
{113, 211}, // ρ0 π+
17+
{225, 211}, // f2(1270) π+
18+
{-211, 211, 211} // π- π+ π+ (non-resonant)
19+
}},
20+
{421, {
21+
{-321, 211}, // K- π+ (non-resonant)
22+
{-321, 111, 211}, // K- π+ π0
23+
{213, -321}, // ρ+ K-
24+
{-313, 111}, // antiK*0(892) π0
25+
{-323, 211}, // K*-(892) π+
26+
{-211, 211}, // π- π+
27+
{213, -211}, // ρ+ π-
28+
{-211, 211, 111}, // π- π+ π0
29+
{-321, 321} // K- K+
30+
}},
31+
{431, {
32+
{211, 333}, // φ π+
33+
{-313, 321}, // antiK*(892) K+
34+
{333, 213}, // φ ρ
35+
{113, 211}, // ρ π+
36+
{225, 211}, // f2(1270) π+
37+
{-211, 211, 211}, // π- π+ π+ (s-wave)
38+
{313, 211}, // K*(892)0 π+
39+
{10221, 321}, // f0(1370) K+
40+
{113, 321}, // ρ0 K+
41+
{-211, 321, 211}, // π- K+ π+ (non-resonant)
42+
{221, 211} // η π+
43+
}},
44+
{4122, {
45+
{2212, -321, 211}, // p K- π+ (non-resonant)
46+
{2212, -313}, // p K*0(892)
47+
{2224, -321}, // Δ++ K-
48+
{102134, 211}, // Λ(1520) K-
49+
{2212, -321, 211, 111}, // p K- π+ π0
50+
{2212, -211, 211}, // p π- π+
51+
{2212, 333} // p φ
52+
}},
53+
{4232, {
54+
{-313, 2212}, // antiK*0(892) p
55+
{2212, -321, 211}, // p K- π+
56+
{2212, 333}, // p φ
57+
{3222, -211, 211} // Σ+ π- π+
58+
}},
59+
};
60+
61+
TFile file(path.c_str(), "READ");
62+
if (file.IsZombie()) {
63+
std::cerr << "Cannot open ROOT file " << path << "\n";
64+
return 1;
65+
}
66+
67+
auto tree = (TTree *)file.Get("o2sim");
68+
std::vector<o2::MCTrack> *tracks{};
69+
tree->SetBranchAddress("MCTrack", &tracks);
70+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
71+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
72+
73+
int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
74+
int nQuarksOne{}, nQuarksTwo{}, nSignals{}, nSignalGoodDecay{};
75+
auto nEvents = tree->GetEntries();
76+
77+
for (int i = 0; i < nEvents; i++) {
78+
tree->GetEntry(i);
79+
80+
// check subgenerator information
81+
if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
82+
bool isValid = false;
83+
int subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
84+
if (subGeneratorId == 0) {
85+
nEventsMB++;
86+
} else if (subGeneratorId == checkPdgQuarkOne) {
87+
nEventsInjOne++;
88+
} else if (subGeneratorId == checkPdgQuarkTwo) {
89+
nEventsInjTwo++;
90+
}
91+
}
92+
93+
for (auto &track : *tracks) {
94+
auto pdg = track.GetPdgCode();
95+
if (std::abs(pdg) == checkPdgQuarkOne) {
96+
nQuarksOne++;
97+
continue;
98+
}
99+
if (std::abs(pdg) == checkPdgQuarkTwo) {
100+
nQuarksTwo++;
101+
continue;
102+
}
103+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) { // found signal
104+
nSignals++; // count signal PDG
105+
106+
std::vector<int> pdgsDecay{};
107+
std::vector<int> pdgsDecayAntiPart{};
108+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
109+
auto pdgDau = tracks->at(j).GetPdgCode();
110+
pdgsDecay.push_back(pdgDau);
111+
if (pdgDau != 333) { // phi is antiparticle of itself
112+
pdgsDecayAntiPart.push_back(-pdgDau);
113+
} else {
114+
pdgsDecayAntiPart.push_back(pdgDau);
115+
}
116+
}
117+
118+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
119+
std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
120+
121+
for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
122+
if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
123+
nSignalGoodDecay++;
124+
break;
125+
}
126+
}
127+
}
128+
}
129+
}
130+
131+
std::cout << "--------------------------------\n";
132+
std::cout << "# Events: " << nEvents << "\n";
133+
std::cout << "# MB events: " << nEventsMB << "\n";
134+
std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
135+
std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
136+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkOne) << nQuarksOne << "\n";
137+
std::cout << Form("# %d (anti)quarks: ", checkPdgQuarkTwo) << nQuarksTwo << "\n";
138+
std::cout <<"# signal hadrons: " << nSignals << "\n";
139+
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
140+
141+
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
142+
std::cerr << "Number of generated MB events different than expected\n";
143+
return 1;
144+
}
145+
if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
146+
std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
147+
return 1;
148+
}
149+
if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
150+
std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
151+
return 1;
152+
}
153+
154+
if (nQuarksOne < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
155+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkOne << " lower than expected\n";
156+
return 1;
157+
}
158+
if (nQuarksTwo < nEvents * ratioTrigger) { // we expect anyway more because the same quark is repeated several time, after each gluon radiation
159+
std::cerr << "Number of generated (anti)quarks " << checkPdgQuarkTwo << " lower than expected\n";
160+
return 1;
161+
}
162+
163+
float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
164+
if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
165+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
166+
return 1;
167+
}
168+
169+
return 0;
170+
}

0 commit comments

Comments
 (0)