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