|
| 1 | +#include <TChain.h> |
| 2 | +#include <TTreeReader.h> |
| 3 | +#include <TTreeReaderArray.h> |
| 4 | +#include <TTreeReaderValue.h> |
| 5 | +#include <TH2D.h> |
| 6 | +#include <TCanvas.h> |
| 7 | +#include <TFile.h> |
| 8 | +#include <TMath.h> |
| 9 | +#include <set> |
| 10 | +#include <iostream> |
| 11 | +#include <TString.h> |
| 12 | + |
| 13 | +using namespace std; |
| 14 | + |
| 15 | +int acceptance_analysis(TString filename, string outname_pdf, string outname_png) |
| 16 | +{ |
| 17 | + TChain *chain = new TChain("events"); |
| 18 | + chain->Add(filename); |
| 19 | + |
| 20 | + TTreeReader reader(chain); |
| 21 | + |
| 22 | + TTreeReaderArray<int> mc_pdg(reader, "MCParticles.PDG"); |
| 23 | + TTreeReaderArray<int> mc_genStatus(reader, "MCParticles.generatorStatus"); |
| 24 | + TTreeReaderArray<double> mc_px(reader, "MCParticles.momentum.x"); |
| 25 | + TTreeReaderArray<double> mc_py(reader, "MCParticles.momentum.y"); |
| 26 | + TTreeReaderArray<double> mc_pz(reader, "MCParticles.momentum.z"); |
| 27 | + |
| 28 | + TTreeReaderArray<int> contrib_particle_idx(reader, "_HcalEndcapNHitsContributions_particle.index"); |
| 29 | + TTreeReaderArray<unsigned int> contrib_particle_cid(reader, "_HcalEndcapNHitsContributions_particle.collectionID"); |
| 30 | + |
| 31 | + int nEtaBins = 100; |
| 32 | + int nPhiBins = 100; |
| 33 | + double etaMin = -5, etaMax = 0; |
| 34 | + |
| 35 | + TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta[1]; #phi[rad]", |
| 36 | + nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi()); |
| 37 | + |
| 38 | + TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta[1]; #phi[rad]", |
| 39 | + nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi()); |
| 40 | + |
| 41 | + while (reader.Next()) |
| 42 | + { |
| 43 | + map<int, pair<double, double>> pi_minus_eta_phi; |
| 44 | + set<int> detected; |
| 45 | + for (size_t i = 0; i < mc_pdg.GetSize(); ++i) |
| 46 | + { |
| 47 | + if (mc_pdg[i] == -211 && mc_genStatus[i] == 1) |
| 48 | + { |
| 49 | + float px = mc_px[i]; |
| 50 | + float py = mc_py[i]; |
| 51 | + float pz = mc_pz[i]; |
| 52 | + |
| 53 | + float p = sqrt(px * px + py * py + pz * pz); |
| 54 | + float eta = 0.5 * log((p + pz) / (p - pz + 1e-8)); |
| 55 | + float phi = atan2(py, px); |
| 56 | + |
| 57 | + hEtaPhiAll->Fill(eta, phi); |
| 58 | + pi_minus_eta_phi[i] = make_pair(eta, phi); |
| 59 | + } |
| 60 | + } |
| 61 | + |
| 62 | + for (size_t i = 0; i < contrib_particle_idx.GetSize(); i++) { |
| 63 | + int idx = contrib_particle_idx[i]; |
| 64 | + if (pi_minus_eta_phi.count(idx)) { |
| 65 | + detected.insert(idx); |
| 66 | + } |
| 67 | + } |
| 68 | + |
| 69 | + for (auto idx : detected) { |
| 70 | + auto [eta, phi] = pi_minus_eta_phi[idx]; |
| 71 | + hEtaPhiDetected->Fill(eta, phi); |
| 72 | + } |
| 73 | + } |
| 74 | + |
| 75 | + TH2D* hAcceptance = (TH2D*)hEtaPhiAll->Clone("hAcceptance"); |
| 76 | + hAcceptance->Divide(hEtaPhiDetected); |
| 77 | + hAcceptance->SetTitle("#pi- detected/All"); |
| 78 | + hAcceptance->SetMinimum(0); |
| 79 | + hAcceptance->SetMaximum(1); |
| 80 | + |
| 81 | + TCanvas *canvas = new TCanvas("canvas", "pi- All", 1600, 600); |
| 82 | + canvas->Divide(3,1); |
| 83 | + canvas->cd(1); |
| 84 | + hEtaPhiAll->Draw("COLZ"); |
| 85 | + canvas->cd(2); |
| 86 | + hEtaPhiDetected->Draw("COLZ"); |
| 87 | + canvas->cd(3); |
| 88 | + hAcceptance->Draw("COLZ"); |
| 89 | + |
| 90 | + canvas->SaveAs(outname_pdf.c_str()); |
| 91 | + canvas->SaveAs(outname_png.c_str()); |
| 92 | + |
| 93 | + return 0; |
| 94 | +} |
| 95 | + |
| 96 | + |
0 commit comments