|
| 1 | +int External() |
| 2 | +{ |
| 3 | + int checkPdgSignal[] = {443,100443}; |
| 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 nSignalJpsi{}; |
| 23 | + int nSignalPsi2S{}; |
| 24 | + int nSignalJpsiWithinAcc{}; |
| 25 | + int nSignalPsi2SWithinAcc{}; |
| 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]) { |
| 43 | + if(idMoth < 0){ |
| 44 | + // count signal PDG |
| 45 | + pdg == checkPdgSignal[0] ? nSignalJpsi++ : nSignalPsi2S++; |
| 46 | + // count signal PDG within acceptance |
| 47 | + if(rapidity > rapiditymin && rapidity < rapiditymax) { pdg == checkPdgSignal[0] ? nSignalJpsiWithinAcc++ : nSignalPsi2SWithinAcc++;} |
| 48 | + } |
| 49 | + auto child0 = o2::mcutils::MCTrackNavigator::getDaughter0(track, *tracks); |
| 50 | + auto child1 = o2::mcutils::MCTrackNavigator::getDaughter1(track, *tracks); |
| 51 | + if (child0 != nullptr && child1 != nullptr) { |
| 52 | + // check for parent-child relations |
| 53 | + auto pdg0 = child0->GetPdgCode(); |
| 54 | + auto pdg1 = child1->GetPdgCode(); |
| 55 | + std::cout << "First and last children of parent " << checkPdgSignal << " are PDG0: " << pdg0 << " PDG1: " << pdg1 << "\n"; |
| 56 | + if (std::abs(pdg0) == checkPdgDecay && std::abs(pdg1) == checkPdgDecay && pdg0 == -pdg1) { |
| 57 | + nLeptonPairs++; |
| 58 | + if (child0->getToBeDone() && child1->getToBeDone()) { |
| 59 | + nLeptonPairsToBeDone++; |
| 60 | + } |
| 61 | + } |
| 62 | + } |
| 63 | + } |
| 64 | + } |
| 65 | + } |
| 66 | + std::cout << "#events: " << nEvents << "\n" |
| 67 | + << "#leptons: " << nLeptons << "\n" |
| 68 | + << "#antileptons: " << nAntileptons << "\n" |
| 69 | + << "#signal (prompt Jpsi): " << nSignalJpsi << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalJpsiWithinAcc << "\n" |
| 70 | + << "#signal (prompt Psi(2S)): " << nSignalPsi2S << "; within acceptance " << rapiditymin << " < y < " << rapiditymax << " : " << nSignalPsi2SWithinAcc << "\n" |
| 71 | + << "#lepton pairs: " << nLeptonPairs << "\n" |
| 72 | + << "#lepton pairs to be done: " << nLeptonPairs << "\n"; |
| 73 | + |
| 74 | + |
| 75 | + if (nLeptonPairs == 0 || nLeptons == 0 || nAntileptons == 0) { |
| 76 | + std::cerr << "Number of leptons, number of anti-leptons as well as number of lepton pairs should all be greater than 1.\n"; |
| 77 | + return 1; |
| 78 | + } |
| 79 | + if (nLeptonPairs != nLeptonPairsToBeDone) { |
| 80 | + std::cerr << "The number of lepton pairs should be the same as the number of lepton pairs which should be transported.\n"; |
| 81 | + return 1; |
| 82 | + } |
| 83 | + |
| 84 | + return 0; |
| 85 | +} |
0 commit comments