Skip to content

Commit 58feaef

Browse files
Adding missing test file
1 parent 3afd31e commit 58feaef

File tree

1 file changed

+101
-0
lines changed

1 file changed

+101
-0
lines changed
Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
int External()
2+
{
3+
int checkPdgSignal[] = {553,100553,200553};
4+
int checkPdgDecay = 13;
5+
double rapiditymin = -4.3; double rapiditymax = -2.3;
6+
std::string path{"o2sim_Kine.root"};
7+
std::cout << "Check for\nsignal PDG " << checkPdgSignal << "\ndecay PDG " << checkPdgDecay << "\n";
8+
TFile file(path.c_str(), "READ");
9+
if (file.IsZombie()) {
10+
std::cerr << "Cannot open ROOT file " << path << "\n";
11+
return 1;
12+
}
13+
14+
auto tree = (TTree*)file.Get("o2sim");
15+
std::vector<o2::MCTrack>* tracks{};
16+
tree->SetBranchAddress("MCTrack", &tracks);
17+
18+
int nLeptons{};
19+
int nAntileptons{};
20+
int nLeptonPairs{};
21+
int nLeptonPairsToBeDone{};
22+
int nSignalUpsilon1S{};
23+
int nSignalUpsilon2S{};
24+
int nSignalUpsilon1SWithinAcc{};
25+
int nSignalUpsilon2SWithinAcc{};
26+
auto nEvents = tree->GetEntries();
27+
o2::steer::MCKinematicsReader mcreader("o2sim", o2::steer::MCKinematicsReader::Mode::kMCKine);
28+
Bool_t isInjected = kFALSE;
29+
30+
for (int i = 0; i < nEvents; i++) {
31+
tree->GetEntry(i);
32+
for (auto& track : *tracks) {
33+
auto pdg = track.GetPdgCode();
34+
auto rapidity = track.GetRapidity();
35+
auto idMoth = track.getMotherTrackId();
36+
if (pdg == checkPdgDecay) {
37+
// count leptons
38+
nLeptons++;
39+
} else if(pdg == -checkPdgDecay) {
40+
// count anti-leptons
41+
nAntileptons++;
42+
} else if (pdg == checkPdgSignal[0] || pdg == checkPdgSignal[1] || pdg == checkPdgSignal[2]) {
43+
if(idMoth < 0){
44+
// count signal PDG
45+
if (pdg == checkPdgSignal[0]) {
46+
nSignalUpsilon1S++;
47+
} else if (pdg == checkPdgSignal[1]) {
48+
nSignalUpsilon2S++;
49+
} else {
50+
nSignalUpsilon3S++;
51+
}
52+
53+
// count signal PDG within acceptance
54+
if (rapidity > rapiditymin && rapidity < rapiditymax) {
55+
if (pdg == checkPdgSignal[0]) {
56+
nSignalUpsilon1SWithinAcc++;
57+
} else if (pdg == checkPdgSignal[1]) {
58+
nSignalUpsilon2SWithinAcc++;
59+
} else {
60+
nSignalUpsilon3SWithinAcc++;
61+
}
62+
}
63+
}
64+
auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks);
65+
auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks);
66+
if (child0 != nullptr && child1 != nullptr) {
67+
// check for parent-child relations
68+
auto pdg0 = child0->GetPdgCode();
69+
auto pdg1 = child1->GetPdgCode();
70+
std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n";
71+
if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) {
72+
nLeptonPairs++;
73+
if (child0->getToBeDone() && child1->getToBeDone()) {
74+
nLeptonPairsToBeDone++;
75+
}
76+
}
77+
}
78+
}
79+
}
80+
}
81+
std::cout << "#events: " << nEvents << "\n"
82+
<< "#leptons: " << nLeptons << "\n"
83+
<< "#antileptons: " << nAntileptons << "\n"
84+
<< "#signal (Upsilon(1S)): " << nSignalUpsilon1S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon1SWithinAcc << "\n"
85+
<< "#signal (Upsilon(2S)): " << nSignalUpsilon2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon2SWithinAcc << "\n"
86+
<< "#signal (Upsilon(2S)): " << nSignalUpsilon3S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalUpsilon3SWithinAcc << "\n"
87+
<< "#lepton pairs: " << nLeptonPairs << "\n"
88+
<< "#lepton pairs to be done: " << nLeptonPairs << "\n";
89+
90+
91+
if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) {
92+
std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n";
93+
return 1;
94+
}
95+
if (nLeptonPairs != nLeptonPairsToBeDone) {
96+
std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n";
97+
return 1;
98+
}
99+
100+
return 0;
101+
}

0 commit comments

Comments
 (0)