Skip to content

Commit 9dc9ea3

Browse files
klsmith15kklsmith15k
authored
Add files via upload (#2199)
missing test macro for PR#2186
1 parent 0e607be commit 9dc9ea3

File tree

1 file changed

+190
-0
lines changed

1 file changed

+190
-0
lines changed

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

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

0 commit comments

Comments
 (0)