File indexing completed on 2025-07-27 07:53:55
0001 #include <TChain.h>
0002 #include <TTreeReader.h>
0003 #include <TTreeReaderArray.h>
0004 #include <TTreeReaderValue.h>
0005 #include <TH2D.h>
0006 #include <TCanvas.h>
0007 #include <TFile.h>
0008 #include <TMath.h>
0009 #include <set>
0010 #include <iostream>
0011 #include <TString.h>
0012
0013 using namespace std;
0014
0015 int acceptance_analysis(TString filename, string outname_pdf, string outname_png)
0016 {
0017 TChain *chain = new TChain("events");
0018 chain->Add(filename);
0019
0020 TTreeReader reader(chain);
0021
0022 TTreeReaderArray<int> mc_pdg(reader, "MCParticles.PDG");
0023 TTreeReaderArray<int> mc_genStatus(reader, "MCParticles.generatorStatus");
0024 TTreeReaderArray<double> mc_px(reader, "MCParticles.momentum.x");
0025 TTreeReaderArray<double> mc_py(reader, "MCParticles.momentum.y");
0026 TTreeReaderArray<double> mc_pz(reader, "MCParticles.momentum.z");
0027
0028 TTreeReaderArray<int> contrib_particle_idx(reader, "_HcalEndcapNHitsContributions_particle.index");
0029 TTreeReaderArray<unsigned int> contrib_particle_cid(reader, "_HcalEndcapNHitsContributions_particle.collectionID");
0030
0031 int nEtaBins = 100;
0032 int nPhiBins = 100;
0033 double etaMin = -5, etaMax = 0;
0034
0035 TH2D* hEtaPhiAll = new TH2D("hEtaPhiAll", "All #pi- (status==1); #eta[1]; #phi[rad]",
0036 nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
0037
0038 TH2D* hEtaPhiDetected = new TH2D("hEtaPhiDetected", "#pi- detected in nHCal; #eta[1]; #phi[rad]",
0039 nEtaBins, etaMin, etaMax, nPhiBins, -TMath::Pi(), TMath::Pi());
0040
0041 while (reader.Next())
0042 {
0043 map<int, pair<double, double>> pi_minus_eta_phi;
0044 set<int> detected;
0045 for (size_t i = 0; i < mc_pdg.GetSize(); ++i)
0046 {
0047 if (mc_pdg[i] == -211 && mc_genStatus[i] == 1)
0048 {
0049 float px = mc_px[i];
0050 float py = mc_py[i];
0051 float pz = mc_pz[i];
0052
0053 float p = sqrt(px * px + py * py + pz * pz);
0054 float eta = 0.5 * log((p + pz) / (p - pz + 1e-8));
0055 float phi = atan2(py, px);
0056
0057 hEtaPhiAll->Fill(eta, phi);
0058 pi_minus_eta_phi[i] = make_pair(eta, phi);
0059 }
0060 }
0061
0062 for (size_t i = 0; i < contrib_particle_idx.GetSize(); i++) {
0063 int idx = contrib_particle_idx[i];
0064 if (pi_minus_eta_phi.count(idx)) {
0065 detected.insert(idx);
0066 }
0067 }
0068
0069 for (auto idx : detected) {
0070 auto [eta, phi] = pi_minus_eta_phi[idx];
0071 hEtaPhiDetected->Fill(eta, phi);
0072 }
0073 }
0074
0075 TH2D* hAcceptance = (TH2D*)hEtaPhiAll->Clone("hAcceptance");
0076 hAcceptance->Divide(hEtaPhiDetected);
0077 hAcceptance->SetTitle("#pi- detected/All");
0078 hAcceptance->SetMinimum(0);
0079 hAcceptance->SetMaximum(1);
0080
0081 TCanvas *canvas = new TCanvas("canvas", "pi- All", 1600, 600);
0082 canvas->Divide(3,1);
0083 canvas->cd(1);
0084 hEtaPhiAll->Draw("COLZ");
0085 canvas->cd(2);
0086 hEtaPhiDetected->Draw("COLZ");
0087 canvas->cd(3);
0088 hAcceptance->Draw("COLZ");
0089
0090 canvas->SaveAs(outname_pdf.c_str());
0091 canvas->SaveAs(outname_png.c_str());
0092
0093 return 0;
0094 }
0095
0096
