Skip to content

Commit 4fb6550

Browse files
author
Florian Jonas
committed
[PWGHF] add configuration for D2H ccbar and bbbar gap2 for pO collisions
1 parent 35d151b commit 4fb6550

File tree

3 files changed

+497
-0
lines changed

3 files changed

+497
-0
lines changed

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap2_SCCR_pO.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_gaptriggered_hf.C
4+
funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(2, -1.5, 1.5)
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_SCCR_pO.cfg
8+
includePartonEvent=true

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

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

0 commit comments

Comments
 (0)