|
| 1 | +int External() |
| 2 | +{ |
| 3 | + std::string path{"o2sim_Kine.root"}; |
| 4 | + int numberOfInjectedSignalsPerEvent{1}; |
| 5 | + int numberOfGapEvents{4}; |
| 6 | + int numberOfEventsProcessed{0}; |
| 7 | + int numberOfEventsProcessedWithoutInjection{0}; |
| 8 | + std::vector<int> injectedPDGs = { |
| 9 | + 102134, // Lambda(1520)0 |
| 10 | + -102134, // Lambda(1520)0bar |
| 11 | + 3324, // Xi(1530)0 |
| 12 | + -3324, // Xi(1530)0bar |
| 13 | + 123314, // Xi(1820)- |
| 14 | + -123314, // Xi(1820)+ |
| 15 | + 123324, // Xi(1820)0 |
| 16 | + -123324 // Xi(1820)0bar |
| 17 | + }; |
| 18 | + std::vector<std::vector<int>> decayDaughters = { |
| 19 | + {2212, -321}, // Lambda(1520)0 |
| 20 | + {-2212, 321}, // Lambda(1520)0bar |
| 21 | + {3312, 211}, // Xi(1530)0 |
| 22 | + {-3312, -211}, // Xi(1530)0bar |
| 23 | + {3122, -321}, // Xi(1820)- |
| 24 | + {-3122, 321}, // Xi(1820)+ |
| 25 | + {3122, 310}, // Xi(1820)0 |
| 26 | + {-3122, -310} // Xi(1820)0bar |
| 27 | + }; |
| 28 | + |
| 29 | + auto nInjection = injectedPDGs.size(); |
| 30 | + |
| 31 | + TFile file(path.c_str(), "READ"); |
| 32 | + if (file.IsZombie()) |
| 33 | + { |
| 34 | + std::cerr << "Cannot open ROOT file " << path << "\n"; |
| 35 | + return 1; |
| 36 | + } |
| 37 | + |
| 38 | + auto tree = (TTree *)file.Get("o2sim"); |
| 39 | + if (!tree) |
| 40 | + { |
| 41 | + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; |
| 42 | + return 1; |
| 43 | + } |
| 44 | + std::vector<o2::MCTrack> *tracks{}; |
| 45 | + tree->SetBranchAddress("MCTrack", &tracks); |
| 46 | + |
| 47 | + std::vector<int> nSignal; |
| 48 | + for (int i = 0; i < nInjection; i++) |
| 49 | + { |
| 50 | + nSignal.push_back(0); |
| 51 | + } |
| 52 | + std::vector<std::vector<int>> nDecays; |
| 53 | + std::vector<int> nNotDecayed; |
| 54 | + for (int i = 0; i < nInjection; i++) |
| 55 | + { |
| 56 | + std::vector<int> nDecay; |
| 57 | + for (int j = 0; j < decayDaughters[i].size(); j++) |
| 58 | + { |
| 59 | + nDecay.push_back(0); |
| 60 | + } |
| 61 | + nDecays.push_back(nDecay); |
| 62 | + nNotDecayed.push_back(0); |
| 63 | + } |
| 64 | + auto nEvents = tree->GetEntries(); |
| 65 | + bool hasInjection = false; |
| 66 | + for (int i = 0; i < nEvents; i++) |
| 67 | + { |
| 68 | + hasInjection = false; |
| 69 | + numberOfEventsProcessed++; |
| 70 | + auto check = tree->GetEntry(i); |
| 71 | + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) |
| 72 | + { |
| 73 | + auto track = tracks->at(idxMCTrack); |
| 74 | + auto pdg = track.GetPdgCode(); |
| 75 | + auto it = std::find(injectedPDGs.begin(), injectedPDGs.end(), pdg); |
| 76 | + int index = std::distance(injectedPDGs.begin(), it); // index of injected PDG |
| 77 | + if (it != injectedPDGs.end()) // found |
| 78 | + { |
| 79 | + // count signal PDG |
| 80 | + nSignal[index]++; |
| 81 | + if (track.getFirstDaughterTrackId() < 0) |
| 82 | + { |
| 83 | + nNotDecayed[index]++; |
| 84 | + continue; |
| 85 | + } |
| 86 | + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) |
| 87 | + { |
| 88 | + auto pdgDau = tracks->at(j).GetPdgCode(); |
| 89 | + bool foundDau = false; |
| 90 | + // count decay PDGs |
| 91 | + for (int idxDaughter = 0; idxDaughter < decayDaughters[index].size(); ++idxDaughter) |
| 92 | + { |
| 93 | + if (pdgDau == decayDaughters[index][idxDaughter]) |
| 94 | + { |
| 95 | + nDecays[index][idxDaughter]++; |
| 96 | + foundDau = true; |
| 97 | + hasInjection = true; |
| 98 | + break; |
| 99 | + } |
| 100 | + } |
| 101 | + if (!foundDau) |
| 102 | + { |
| 103 | + std::cerr << "Decay daughter not found: " << pdg << " -> " << pdgDau << "\n"; |
| 104 | + } |
| 105 | + } |
| 106 | + } |
| 107 | + } |
| 108 | + if (!hasInjection) |
| 109 | + { |
| 110 | + numberOfEventsProcessedWithoutInjection++; |
| 111 | + } |
| 112 | + } |
| 113 | + std::cout << "--------------------------------\n"; |
| 114 | + std::cout << "# Events: " << nEvents << "\n"; |
| 115 | + for (int i = 0; i < nInjection; i++) |
| 116 | + { |
| 117 | + std::cout << "# Mother \n"; |
| 118 | + std::cout << injectedPDGs[i] << " generated: " << nSignal[i] << ", " << nNotDecayed[i] << " did not decay\n"; |
| 119 | + if (nSignal[i] == 0) |
| 120 | + { |
| 121 | + std::cerr << "No generated: " << injectedPDGs[i] << "\n"; |
| 122 | + // return 1; // At least one of the injected particles should be generated |
| 123 | + } |
| 124 | + for (int j = 0; j < decayDaughters[i].size(); j++) |
| 125 | + { |
| 126 | + std::cout << "# Daughter " << decayDaughters[i][j] << ": " << nDecays[i][j] << "\n"; |
| 127 | + } |
| 128 | + // if (nSignal[i] != nEvents * numberOfInjectedSignalsPerEvent) |
| 129 | + // { |
| 130 | + // std::cerr << "Number of generated: " << injectedPDGs[i] << ", lower than expected\n"; |
| 131 | + // // return 1; // Don't need to return 1, since the number of generated particles is not the same for each event |
| 132 | + // } |
| 133 | + } |
| 134 | + std::cout << "--------------------------------\n"; |
| 135 | + std::cout << "Number of events processed: " << numberOfEventsProcessed << "\n"; |
| 136 | + std::cout << "Number of input for the gap events: " << numberOfGapEvents << "\n"; |
| 137 | + std::cout << "Number of events processed without injection: " << numberOfEventsProcessedWithoutInjection << "\n"; |
| 138 | + // injected event + numberOfGapEvents*gap events + injected event + numberOfGapEvents*gap events + ... |
| 139 | + // total fraction of the gap event: numberOfEventsProcessedWithoutInjection/numberOfEventsProcessed |
| 140 | + float ratioOfNormalEvents = numberOfEventsProcessedWithoutInjection / numberOfEventsProcessed; |
| 141 | + if (ratioOfNormalEvents > 0.75) |
| 142 | + { |
| 143 | + std::cout << "The number of injected event is loo low!!" << std::endl; |
| 144 | + return 1; |
| 145 | + } |
| 146 | + |
| 147 | + return 0; |
| 148 | +} |
| 149 | + |
| 150 | +void GeneratorLF_Resonances_pp1360_injection() { External(); } |
0 commit comments