|
| 1 | + |
| 2 | +int External() |
| 3 | +{ |
| 4 | + const std::string path{"o2sim_Kine.root"}; |
| 5 | + const int numberOfGapEvents{4}; // generateLFRapidity(..., gap=5) |
| 6 | + const std::vector<int> injectedPDGs = { |
| 7 | + 9010221, // f0(980) |
| 8 | + 3324, -3324, // Xi(1530)0 and anti-Xi(1530)0 |
| 9 | + 123324, -123324, // Xi(1820)0 and anti-Xi(1820)0 |
| 10 | + 123314, -123314, // Xi(1820)- and Xi(1820)+ |
| 11 | + 123334, -123334 // Omega(2012)- and Omega(2012)+ |
| 12 | + }; |
| 13 | + const std::vector<std::vector<int>> decayDaughters = { |
| 14 | + {211, -211}, |
| 15 | + {3312, 211}, {-3312, -211}, |
| 16 | + {3122, 310}, {-3122, 310}, |
| 17 | + {3122, -321}, {-3122, 321}, |
| 18 | + {3312, 310}, {-3312, 310} |
| 19 | + }; |
| 20 | + |
| 21 | + TFile file(path.c_str(), "READ"); |
| 22 | + if (file.IsZombie()) { |
| 23 | + std::cerr << "Cannot open ROOT file " << path << "\n"; |
| 24 | + return 1; |
| 25 | + } |
| 26 | + auto tree = (TTree*)file.Get("o2sim"); |
| 27 | + if (!tree) { |
| 28 | + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; |
| 29 | + return 1; |
| 30 | + } |
| 31 | + std::vector<o2::MCTrack>* tracks{}; |
| 32 | + tree->SetBranchAddress("MCTrack", &tracks); |
| 33 | + |
| 34 | + std::vector<int> nSignal(injectedPDGs.size(), 0); |
| 35 | + std::vector<int> nNotDecayed(injectedPDGs.size(), 0); |
| 36 | + std::vector<std::vector<int>> nDecays; |
| 37 | + for (const auto& daughters : decayDaughters) { |
| 38 | + nDecays.emplace_back(daughters.size(), 0); |
| 39 | + } |
| 40 | + |
| 41 | + int numberOfEventsProcessed{0}; |
| 42 | + int numberOfEventsProcessedWithoutInjection{0}; |
| 43 | + for (Long64_t i = 0; i < tree->GetEntries(); ++i) { |
| 44 | + tree->GetEntry(i); |
| 45 | + ++numberOfEventsProcessed; |
| 46 | + bool hasInjection{false}; |
| 47 | + for (const auto& track : *tracks) { |
| 48 | + const auto pdg = track.GetPdgCode(); |
| 49 | + const auto it = std::find(injectedPDGs.begin(), injectedPDGs.end(), pdg); |
| 50 | + if (it == injectedPDGs.end()) { |
| 51 | + continue; |
| 52 | + } |
| 53 | + const auto index = static_cast<size_t>(std::distance(injectedPDGs.begin(), it)); |
| 54 | + ++nSignal[index]; |
| 55 | + if (track.getFirstDaughterTrackId() < 0) { |
| 56 | + ++nNotDecayed[index]; |
| 57 | + continue; |
| 58 | + } |
| 59 | + for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j) { |
| 60 | + const auto pdgDau = tracks->at(j).GetPdgCode(); |
| 61 | + bool foundDau{false}; |
| 62 | + for (size_t k = 0; k < decayDaughters[index].size(); ++k) { |
| 63 | + if (pdgDau == decayDaughters[index][k]) { |
| 64 | + ++nDecays[index][k]; |
| 65 | + foundDau = true; |
| 66 | + hasInjection = true; |
| 67 | + break; |
| 68 | + } |
| 69 | + } |
| 70 | + if (!foundDau) { |
| 71 | + std::cerr << "Decay daughter not found: " << pdg << " -> " << pdgDau |
| 72 | + << " (mother=" << track.getMotherTrackId() |
| 73 | + << ", secondMother=" << track.getSecondMotherTrackId() << ")\n"; |
| 74 | + } |
| 75 | + } |
| 76 | + } |
| 77 | + if (!hasInjection) { |
| 78 | + ++numberOfEventsProcessedWithoutInjection; |
| 79 | + } |
| 80 | + } |
| 81 | + |
| 82 | + std::cout << "--------------------------------\n"; |
| 83 | + std::cout << "# Events: " << tree->GetEntries() << "\n"; |
| 84 | + for (size_t i = 0; i < injectedPDGs.size(); ++i) { |
| 85 | + std::cout << "# Mother\n"; |
| 86 | + std::cout << injectedPDGs[i] << " generated: " << nSignal[i] |
| 87 | + << ", " << nNotDecayed[i] << " did not decay\n"; |
| 88 | + for (size_t j = 0; j < decayDaughters[i].size(); ++j) { |
| 89 | + std::cout << "# Daughter " << decayDaughters[i][j] << ": " << nDecays[i][j] << "\n"; |
| 90 | + } |
| 91 | + } |
| 92 | + std::cout << "--------------------------------\n"; |
| 93 | + std::cout << "Number of events processed: " << numberOfEventsProcessed << "\n"; |
| 94 | + std::cout << "Number of input for the gap events: " << numberOfGapEvents << "\n"; |
| 95 | + std::cout << "Number of events processed without injection: " |
| 96 | + << numberOfEventsProcessedWithoutInjection << "\n"; |
| 97 | + const double ratioOfNormalEvents = numberOfEventsProcessed |
| 98 | + ? static_cast<double>(numberOfEventsProcessedWithoutInjection) / |
| 99 | + numberOfEventsProcessed |
| 100 | + : 0.0; |
| 101 | + std::cout << "Fraction without injection: " << ratioOfNormalEvents << "\n"; |
| 102 | + const double expectedRatio = static_cast<double>(numberOfGapEvents) / (numberOfGapEvents + 1); |
| 103 | + std::cout << "Expected fraction for 1+" << numberOfGapEvents << " pattern: " << expectedRatio << "\n"; |
| 104 | + |
| 105 | + // Same basic gap sanity check as the referenced O2DPG test. |
| 106 | + if (ratioOfNormalEvents > 0.90 || ratioOfNormalEvents < 0.70) { |
| 107 | + std::cerr << "The number of injected events is too low or too high\n"; |
| 108 | + return 1; |
| 109 | + } |
| 110 | + return 0; |
| 111 | +} |
0 commit comments