macros/tracking_analysis/Plot-SVTX_Efficiency.cpp

0001 #include <iostream>
0002 #include <cstdlib>
0003 #include <memory>
0004 #include <sstream>
0005 #include <string>
0006 #include <cstring>
0007 #include <vector>
0008 #include <unordered_map>
0009 #include <cmath>
0010 #include <utility>
0011
0012 #include "TFile.h"
0013 #include "TTree.h"
0014 #include "TCanvas.h"
0015 #include "TLegend.h"
0016 #include "TROOT.h"
0017 #include "TH1F.h"
0018 #include "TColor.h"
0019 #include "TImage.h"
0020 #include "TApplication.h"
0021 #include "TGraphErrors.h"
0022 #include "/sphenix/user/gregtom3/SBU/research/macros/macros/sPHENIXStyle/sPhenixStyle.C"
0023
0024 TTree *load_tree(const char *const file_name, const char *const tree_name);
0025 const char *const svtx_file_path {"/sphenix/user/giorgian/muons/svtx.root"};
0026
0027 double MOMENTUM_MARGIN {0.01};
0028
0029 void Plot_SVTX_Efficiency()
0030 {
0031     SetsPhenixStyle();
0032     gROOT->SetBatch(true);
0033
0034     TTree *const ntp_gtrack {load_tree(svtx_file_path, "ntp_gtrack")};
0035     TTree *const ntp_track {load_tree(svtx_file_path, "ntp_track")};
0036
0037     Long64_t ngtrack {ntp_gtrack->GetEntries()};
0038     Long64_t ntrack {ntp_track->GetEntries()};
0039
0040     TH1F *h_true_count {new TH1F("h_true_count", "Psuedorapidity count", 100, -5, 5)};
0041     h_true_count->SetXTitle("#eta");
0042     h_true_count->SetYTitle("Count");
0043     h_true_count->SetLineColor(kBlue);
0044     Float_t geta;
0045     ntp_gtrack->SetBranchAddress("geta", &geta);
0046     for (Long64_t i {0}; i < ngtrack; ++i) {
0047         if (ntp_gtrack->LoadTree(i) < 0)
0048             break;
0049         ntp_gtrack->GetEntry(i);
0050         h_true_count->Fill(geta);
0051     }
0052     TH1F *h_reco_count {new TH1F("h_reco_count", "Psuedorapidity count", 100, -5, 5)};
0053     h_reco_count->SetXTitle("#eta");
0054     h_reco_count->SetYTitle("Count");
0055     h_reco_count->SetLineColor(kRed);
0056     Float_t gpx, gpy, gpz, px, py, pz;
0057     ntp_track->SetBranchAddress("geta", &geta);
0058     ntp_track->SetBranchAddress("gpx", &gpx);
0059     ntp_track->SetBranchAddress("gpy", &gpy);
0060     ntp_track->SetBranchAddress("gpz", &gpz);
0061     ntp_track->SetBranchAddress("px", &px);
0062     ntp_track->SetBranchAddress("py", &py);
0063     ntp_track->SetBranchAddress("pz", &pz);
0064     for (Long64_t i {0}; i < ntrack; ++i) {
0065         if (ntp_track->LoadTree(i) < 0)
0066             break;
0067
0068         ntp_track->GetEntry(i);
0069         if (fabs((px - gpx) / gpx) < MOMENTUM_MARGIN
0070             && fabs((py - gpy) / gpy) < MOMENTUM_MARGIN
0071             && fabs((pz - gpz) / gpz) < MOMENTUM_MARGIN)
0072             h_reco_count->Fill(geta);
0073     }
0074
0075     const Long64_t nbins {h_true_count->GetSize() - 2}; /* -2 for underflow and overflow */
0076     Double_t *x {new Double_t[nbins]};
0077     Double_t *y {new Double_t[nbins]};
0078     Double_t *ey {new Double_t[nbins]};
0079     Long64_t top {0};
0080     for (Long64_t i {0}; i < nbins; ++i) {
0081         if (h_true_count->GetBinContent(i + 1) != 0) {
0082             const Double_t n {h_reco_count->GetBinContent(i + 1)};
0083             const Double_t N {h_true_count->GetBinContent(i + 1)};
0084
0085             x[top] = {h_true_count->GetBinCenter(i + 1)};
0086             y[top] = {n / N};
0087             ey[top++] = {sqrt(n / (N * N) + (n * n) / (N * N * N))};
0088         }
0089     }
0090     TGraphErrors *gr {new TGraphErrors(top, x, y, nullptr, ey)};
0091     gr->SetMarkerColor(kBlue);
0092     gr->SetMarkerStyle(21);
0093     gr->SetMarkerSize(0.5);
0094     gr->GetXaxis()->SetTitle("#eta");
0095     gr->GetYaxis()->SetTitle("Efficiency");
0096
0097     TCanvas *count {new TCanvas("count", "SVTX Event Count",
0098                 gStyle->GetCanvasDefW(), gStyle->GetCanvasDefH())};
0099     h_true_count->GetYaxis()->SetRangeUser(0, 2000);
0100     h_true_count->Draw();
0101     h_reco_count->Draw("SAME");
0102     TLegend *l1 {new TLegend(0.825, .90, .95, 0.80, "Track")};
0103     l1->SetTextSize(0.03);
0104     l1->AddEntry(h_true_count, "True", "l");
0105     l1->AddEntry(h_reco_count, "Reco", "l");
0106     l1->Draw();
0107     gPad->RedrawAxis();
0108
0109     TCanvas *efficiency {new TCanvas("efficiency", "SVTX Efficiency",
0110                                 gStyle->GetCanvasDefW(), gStyle->GetCanvasDefH())};
0111     gr->GetYaxis()->SetRangeUser(-0.05, 1);
0112     gr->Draw("ALP");
0113     TLegend *l2 {new TLegend(0.825, 0.90, 0.95, 0.8, "Track")};
0114     l2->SetTextSize(0.03);
0115     l2->AddEntry(gr, "Efficiency", "l");
0116     l2->Draw();
0117     gPad->RedrawAxis();
0118
0119     TImage *const img {TImage::Create()};
0120     img->FromPad(count);
0121     std::stringstream name;
0122     name << "SVTX_Event_Count-" << "margin=" << MOMENTUM_MARGIN << ".png";
0123     img->WriteImage(strdup(name.str().c_str()));
0124
0125     name.str("");
0126     img->FromPad(efficiency);
0127     name << "SVTX_Efficiency-" << "margin=" << MOMENTUM_MARGIN << ".png";
0128     img->WriteImage(strdup(name.str().c_str()));
0129     gApplication->Terminate(0);
0130 }
0131
0132 TTree *load_tree(const char *const file_name, const char *const tree_name)
0133 {
0134     return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
0135 }
0136
0137 int main(int argc, char *argv[])
0138 {
0139     TApplication app {"SVTX Efficiency Plots", &argc, argv};
0140     if (argc > 1) {
0141         std::stringstream tmp {argv[1]};
0142         tmp >> MOMENTUM_MARGIN;
0143     }
0144
0145     Plot_SVTX_Efficiency();
0146     app.Run();
0147     return 0;
0148 }