|
| 1 | +// Test of the single-channel ML signal samples GeneratorLF_ResonancesML_*_pp1360.ini that redirect here. |
| 2 | +// Every event carries exactly one injected parent (no gap); the parent and the intermediate resonances decayed by the |
| 3 | +// generator must follow the forced chain (charge conjugated for the negative-PDG parent). Parents with a mother come |
| 4 | +// from the underlying Pythia event, decay with the default table, and are only counted. |
| 5 | +int External() |
| 6 | +{ |
| 7 | + const std::string path{"o2sim_Kine.root"}; |
| 8 | + const int parentPDG{10323}; |
| 9 | + const std::vector<int> daughtersOfPositive{313, 211}; |
| 10 | + const std::map<int, std::vector<int>> intermediateDecaysOfPositive{{313, {321, -211}}}; |
| 11 | + const std::set<int> selfConjugate{111, 113, 310}; |
| 12 | + auto conjugate = [&](int pdg, int sign) { return (sign > 0 || selfConjugate.count(pdg)) ? pdg : -pdg; }; |
| 13 | + auto sorted = [](std::vector<int> v) { std::sort(v.begin(), v.end()); return v; }; |
| 14 | + |
| 15 | + TFile file(path.c_str(), "READ"); |
| 16 | + if (file.IsZombie()) { |
| 17 | + std::cerr << "Cannot open ROOT file " << path << "\n"; |
| 18 | + return 1; |
| 19 | + } |
| 20 | + auto tree = (TTree*)file.Get("o2sim"); |
| 21 | + if (!tree) { |
| 22 | + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; |
| 23 | + return 1; |
| 24 | + } |
| 25 | + std::vector<o2::MCTrack>* tracks{}; |
| 26 | + tree->SetBranchAddress("MCTrack", &tracks); |
| 27 | + |
| 28 | + int nEventsWrongParentCount{0}, nFromUnderlyingEvent{0}, nParents{0}, nNotDecayed{0}, nWrongDecay{0}, nIntermediates{0}, nWrongIntermediate{0}; |
| 29 | + std::map<int, int> nParentsBySign; |
| 30 | + for (Long64_t i = 0; i < tree->GetEntries(); ++i) { |
| 31 | + tree->GetEntry(i); |
| 32 | + int nParentsInEvent{0}; |
| 33 | + for (const auto& track : *tracks) { |
| 34 | + const int pdg = track.GetPdgCode(); |
| 35 | + if (std::abs(pdg) != parentPDG) { |
| 36 | + continue; |
| 37 | + } |
| 38 | + if (track.getMotherTrackId() >= 0) { |
| 39 | + ++nFromUnderlyingEvent; |
| 40 | + continue; |
| 41 | + } |
| 42 | + ++nParentsInEvent; |
| 43 | + ++nParents; |
| 44 | + const int sign = pdg > 0 ? 1 : -1; |
| 45 | + ++nParentsBySign[sign]; |
| 46 | + if (track.getFirstDaughterTrackId() < 0) { |
| 47 | + ++nNotDecayed; |
| 48 | + continue; |
| 49 | + } |
| 50 | + std::vector<int> expected, found; |
| 51 | + for (int d : daughtersOfPositive) { |
| 52 | + expected.push_back(conjugate(d, sign)); |
| 53 | + } |
| 54 | + for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j) { |
| 55 | + const auto& daughter = tracks->at(j); |
| 56 | + found.push_back(daughter.GetPdgCode()); |
| 57 | + const auto it = intermediateDecaysOfPositive.find(sign > 0 ? daughter.GetPdgCode() : conjugate(daughter.GetPdgCode(), -1)); |
| 58 | + if (it == intermediateDecaysOfPositive.end()) { |
| 59 | + continue; |
| 60 | + } |
| 61 | + ++nIntermediates; |
| 62 | + std::vector<int> expectedGrand, foundGrand; |
| 63 | + for (int g : it->second) { |
| 64 | + expectedGrand.push_back(conjugate(g, sign)); |
| 65 | + } |
| 66 | + for (int k = daughter.getFirstDaughterTrackId(); k >= 0 && k <= daughter.getLastDaughterTrackId(); ++k) { |
| 67 | + foundGrand.push_back(tracks->at(k).GetPdgCode()); |
| 68 | + } |
| 69 | + if (sorted(expectedGrand) != sorted(foundGrand)) { |
| 70 | + ++nWrongIntermediate; |
| 71 | + } |
| 72 | + } |
| 73 | + if (sorted(expected) != sorted(found)) { |
| 74 | + ++nWrongDecay; |
| 75 | + } |
| 76 | + } |
| 77 | + if (nParentsInEvent != 1) { |
| 78 | + ++nEventsWrongParentCount; |
| 79 | + } |
| 80 | + } |
| 81 | + |
| 82 | + std::cout << "--------------------------------\n"; |
| 83 | + std::cout << "# Events: " << tree->GetEntries() << "\n"; |
| 84 | + std::cout << "# Parents " << parentPDG << ": " << nParentsBySign[1] << ", anti: " << nParentsBySign[-1] << "\n"; |
| 85 | + std::cout << "# Parents from the underlying event (not checked): " << nFromUnderlyingEvent << "\n"; |
| 86 | + std::cout << "# Events without exactly one parent: " << nEventsWrongParentCount << "\n"; |
| 87 | + std::cout << "# Parents not decayed: " << nNotDecayed << ", with wrong daughters: " << nWrongDecay << "\n"; |
| 88 | + std::cout << "# Intermediate resonances: " << nIntermediates << ", with wrong daughters: " << nWrongIntermediate << "\n"; |
| 89 | + std::cout << "--------------------------------\n"; |
| 90 | + if (tree->GetEntries() == 0 || nEventsWrongParentCount || nNotDecayed || nWrongDecay || nWrongIntermediate || |
| 91 | + (!intermediateDecaysOfPositive.empty() && nIntermediates != nParents)) { |
| 92 | + std::cerr << "Forced decay chain not reproduced\n"; |
| 93 | + return 1; |
| 94 | + } |
| 95 | + return 0; |
| 96 | +} |
0 commit comments