Skip to content

Commit 2eb0a5e

Browse files
author
Mattia Faggin
committed
Add test macro.
1 parent 9ee68f5 commit 2eb0a5e

File tree

1 file changed

+200
-0
lines changed

1 file changed

+200
-0
lines changed

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

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
#include "TFile.h"
2+
#include "TTree.h"
3+
4+
#include <string>
5+
#include <iostream>
6+
7+
int External() {
8+
std::string path{"/home/mattia/Documenti/cernbox/Documents/PostDoc/D2H/MC/corrBkgSigmaC/tf1/genevents_Kine.root"};
9+
10+
int checkPdgQuarkOne{4};
11+
int checkPdgQuarkTwo{5};
12+
float ratioTrigger = 1./5.; // one event triggered out of 5
13+
std::array<std::array<int, 2>, 2> pdgReplParticles = {std::array{413, 14122}, std::array{413, 4124}};
14+
std::array<std::array<int, 2>, 2> pdgReplPartCounters = {std::array{0, 0}, std::array{0, 0}};
15+
std::array<float, 2> freqRepl = {0.5, 0.5};
16+
std::map<int, int> sumOrigReplacedParticles = {{413, 0}};
17+
18+
std::array<int, 2> checkPdgHadron{14122, 4124};
19+
std::map<int, std::vector<std::vector<int>>> checkHadronDecays{ // sorted (!) pdg of daughters
20+
//{14122, {{4222, -211}, {4112, 211}, {4122, 211, -211}}}, // Lc(2595)+
21+
//{4124, {{4222, -211}, {4112, 211}, {4122, 211, -211}}} // Lc(2625)+
22+
{14122, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}}, // Lc(2595)+
23+
{4124, {{-211, 4222}, {211, 4112}, {-211, 211, 4122}}} // Lc(2625)+
24+
};
25+
26+
TFile file(path.c_str(), "READ");
27+
if (file.IsZombie()) {
28+
std::cerr << "Cannot open ROOT file " << path << "\n";
29+
return 1;
30+
}
31+
32+
auto tree = (TTree *)file.Get("o2sim");
33+
std::vector<o2::MCTrack> *tracks{};
34+
tree->SetBranchAddress("MCTrack", &tracks);
35+
o2::dataformats::MCEventHeader *eventHeader = nullptr;
36+
tree->SetBranchAddress("MCEventHeader.", &eventHeader);
37+
38+
int nEventsMB{}, nEventsInjOne{}, nEventsInjTwo{};
39+
int nSignals{}, nSignalGoodDecay{};
40+
auto nEvents = tree->GetEntries();
41+
42+
for (int i = 0; i < nEvents; i++) {
43+
44+
std::cout << std::endl;
45+
46+
tree->GetEntry(i);
47+
48+
// check subgenerator information
49+
int subGeneratorId{-1};
50+
if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) {
51+
bool isValid = false;
52+
subGeneratorId = eventHeader->getInfo<int>(o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid);
53+
if (subGeneratorId == 0) {
54+
nEventsMB++;
55+
} else if (subGeneratorId == checkPdgQuarkOne) {
56+
nEventsInjOne++;
57+
} else if (subGeneratorId == checkPdgQuarkTwo) {
58+
nEventsInjTwo++;
59+
}
60+
}
61+
62+
for (auto &track : *tracks) {
63+
auto pdg = track.GetPdgCode();
64+
auto absPdg = std::abs(pdg);
65+
if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), absPdg) != checkPdgHadron.end()) { // found signal
66+
nSignals++; // count signal PDG
67+
std::cout << "==> signal " << absPdg << " found!" << std::endl;
68+
69+
if (subGeneratorId == checkPdgQuarkOne) { // replacement only for prompt ---> BUT ALSO NON-PROMPT D* SEEM TO BE REPLACED
70+
for (int iRepl{0}; iRepl<2; ++iRepl) {
71+
if (absPdg == pdgReplParticles[iRepl][0]) {
72+
pdgReplPartCounters[iRepl][0]++;
73+
sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
74+
} else if (absPdg == pdgReplParticles[iRepl][1]) {
75+
pdgReplPartCounters[iRepl][1]++;
76+
sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
77+
}
78+
}
79+
} else if (subGeneratorId == checkPdgQuarkTwo) {
80+
std::cout << " NB: we have a " << absPdg << " also in event with quark " << checkPdgQuarkTwo << std::endl;
81+
std::cout << " ### mother indices: ";
82+
int idFirstMother = track.getMotherTrackId();
83+
int idSecondMother = track.getSecondMotherTrackId();
84+
std::vector<int> motherIds = {};
85+
for(int i=idFirstMother; i<=idSecondMother; i++) {
86+
std::cout << i << " ";
87+
motherIds.push_back(i);
88+
}
89+
bool partonicEventOn = false;
90+
if(motherIds != std::vector<int>{-1, -1}) {
91+
std::cout << "The " << absPdg << " particle has mothers. This should mean that it comes directly from parton hadronization, and that the partonic event was kept in the MC production " << std::endl;
92+
partonicEventOn = true;
93+
}
94+
std::cout << " ### mother PDG codes: ";
95+
std::vector<int> motherPdgCodes = {};
96+
if(partonicEventOn) {
97+
for(int i=idFirstMother; i<=idSecondMother; i++) {
98+
motherPdgCodes.push_back(tracks->at(i).GetPdgCode());
99+
std::cout << motherPdgCodes.back() << " ";
100+
}
101+
102+
/// check that among the mothers there is a c/cbar quark
103+
/// This means that the charm hadron comes from the c-quark hadronization, where the c/cbar quark
104+
/// comes from a c-cbar pair present in the current event, tagged with a b-bbar (e.g. double-parton scattering)
105+
if(std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 4) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -4) == motherPdgCodes.end()) {
106+
/// if the partinc event is not really saved and we arrive here, it means that motherIds != {-1, -1} because
107+
/// the hadron comes from the decay of a beauty hadron. This can happen if and only if this is not a replaced one (i.e. native from Lambdab0 decay)
108+
if (std::find(motherPdgCodes.begin(), motherPdgCodes.end(), 5122) == motherPdgCodes.end() && std::find(motherPdgCodes.begin(), motherPdgCodes.end(), -5122) == motherPdgCodes.end()) {
109+
std::cerr << "The particle " << absPdg << " does not originate neither from a c/c-bar quark (replaced) nor from a Lambda_b0 decay. There is something wrong, aborting..." << std::endl;
110+
return 1;
111+
}
112+
}
113+
}
114+
std::cout << std::endl;
115+
116+
/// only if we arrive here it means that everything is ok, and we can safely update the counters for the final statistics
117+
for (int iRepl{0}; iRepl<2; ++iRepl) {
118+
if (absPdg == pdgReplParticles[iRepl][0]) {
119+
pdgReplPartCounters[iRepl][0]++;
120+
sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
121+
} else if (absPdg == pdgReplParticles[iRepl][1]) {
122+
pdgReplPartCounters[iRepl][1]++;
123+
sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]++;
124+
}
125+
}
126+
}
127+
128+
129+
std::vector<int> pdgsDecay{};
130+
std::vector<int> pdgsDecayAntiPart{};
131+
if (track.getFirstDaughterTrackId() >= 0 && track.getLastDaughterTrackId() >= 0) {
132+
for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) {
133+
auto pdgDau = tracks->at(j).GetPdgCode();
134+
pdgsDecay.push_back(pdgDau);
135+
std::cout << " -- daughter " << j << ": " << pdgDau << std::endl;
136+
if (pdgDau != 333) { // phi is antiparticle of itself
137+
pdgsDecayAntiPart.push_back(-pdgDau);
138+
} else {
139+
pdgsDecayAntiPart.push_back(pdgDau);
140+
}
141+
}
142+
}
143+
144+
std::sort(pdgsDecay.begin(), pdgsDecay.end());
145+
std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end());
146+
147+
for (auto &decay : checkHadronDecays[std::abs(pdg)]) {
148+
if (pdgsDecay == decay || pdgsDecayAntiPart == decay) {
149+
nSignalGoodDecay++;
150+
std::cout << " !!! GOOD DECAY FOUND !!!" << std::endl;
151+
break;
152+
}
153+
}
154+
}
155+
} // end loop over tracks
156+
}
157+
158+
std::cout << "--------------------------------\n";
159+
std::cout << "# Events: " << nEvents << "\n";
160+
std::cout << "# MB events: " << nEventsMB << "\n";
161+
std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkOne) << nEventsInjOne << "\n";
162+
std::cout << Form("# events injected with %d quark pair: ", checkPdgQuarkTwo) << nEventsInjTwo << "\n";
163+
std::cout <<"# signal hadrons: " << nSignals << "\n";
164+
std::cout <<"# signal hadrons decaying in the correct channel: " << nSignalGoodDecay << "\n";
165+
166+
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
167+
std::cerr << "Number of generated MB events different than expected\n";
168+
return 1;
169+
}
170+
if (nEventsInjOne < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjOne > nEvents * ratioTrigger * 0.5 * 1.05) {
171+
std::cerr << "Number of generated events injected with " << checkPdgQuarkOne << " different than expected\n";
172+
return 1;
173+
}
174+
if (nEventsInjTwo < nEvents * ratioTrigger * 0.5 * 0.95 || nEventsInjTwo > nEvents * ratioTrigger * 0.5 * 1.05) {
175+
std::cerr << "Number of generated events injected with " << checkPdgQuarkTwo << " different than expected\n";
176+
return 1;
177+
}
178+
179+
float fracForcedDecays = float(nSignalGoodDecay) / nSignals;
180+
if (fracForcedDecays < 0.9) { // we put some tolerance (e.g. due to oscillations which might change the final state)
181+
std::cerr << "Fraction of signals decaying into the correct channel " << fracForcedDecays << " lower than expected\n";
182+
return 1;
183+
}
184+
185+
for (int iRepl{0}; iRepl<2; ++iRepl) {
186+
187+
std::cout << " --- pdgReplPartCounters[" << iRepl << "][1] = " << pdgReplPartCounters[iRepl][1] << ", freqRepl[" << iRepl <<"] = " << freqRepl[iRepl] << ", sumOrigReplacedParticles[pdgReplParticles[" << iRepl << "][0]] =" << sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] << std::endl;
188+
189+
if (std::abs(pdgReplPartCounters[iRepl][1] - freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]]) > 2 * std::sqrt(freqRepl[iRepl] * sumOrigReplacedParticles[pdgReplParticles[iRepl][0]])) { // 2 sigma compatibility
190+
float fracMeas = 0.;
191+
if (sumOrigReplacedParticles[pdgReplParticles[iRepl][0]] > 0.) {
192+
fracMeas = float(pdgReplPartCounters[iRepl][1]) / sumOrigReplacedParticles[pdgReplParticles[iRepl][0]];
193+
}
194+
std::cerr << "Fraction of replaced " << pdgReplParticles[iRepl][0] << " into " << pdgReplParticles[iRepl][1] << " is " << fracMeas <<" (expected "<< freqRepl[iRepl] << ")\n";
195+
return 1;
196+
}
197+
}
198+
199+
return 0;
200+
}

0 commit comments

Comments
 (0)