macros/tracking_analysis/Plot-FastTrack_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 fasttrack_file_path
0026         {"/sphenix/user/gregtom3/data/Summer2018/fasttrack-efficiency-10GeV/fasttrack.root"};
0027 double MOMENTUM_MARGIN {0.01};
0028
0029 void Plot_FastTrack_Efficiency()
0030 {
0031     SetsPhenixStyle();
0032     gROOT->SetBatch(true);
0033
0034     TTree *const tracks {load_tree(fasttrack_file_path, "tracks")};
0035
0036     Long64_t ntracks {tracks->GetEntries()};
0037
0038     TH1F *h_true_count {new TH1F("h_true_count", "Psuedorapidity count", 100, -5, 5)};
0039     h_true_count->SetXTitle("#eta");
0040     h_true_count->SetYTitle("Count");
0041     h_true_count->SetLineColor(kBlue);
0042     TH1F *h_reco_count {new TH1F("h_reco_count", "Psuedorapidity count", 100, -5, 5)};
0043     h_reco_count->SetXTitle("#eta");
0044     h_reco_count->SetYTitle("Count");
0045     h_reco_count->SetLineColor(kRed);
0046     Float_t gpx, gpy, gpz, px, py, pz;
0047     tracks->SetBranchAddress("gpx", &gpx);
0048     tracks->SetBranchAddress("gpy", &gpy);
0049     tracks->SetBranchAddress("gpz", &gpz);
0050     tracks->SetBranchAddress("px", &px);
0051     tracks->SetBranchAddress("py", &py);
0052     tracks->SetBranchAddress("pz", &pz);
0053     for (Long64_t i {0}; i < ntracks; ++i) {
0054         if (tracks->LoadTree(i) < 0)
0055             break;
0056         tracks->GetEntry(i);
0057         const double mom {std::sqrt(px * px + py * py + pz * pz)};
0058         const double eta {0.5 *std::log((mom + pz) / (mom - pz))};
0059         const double gmom {std::sqrt(gpx * gpx + gpy * gpy + gpz * gpz)};
0060         const double geta {0.5 *std::log((gmom + gpz) / (gmom - gpz))};
0061
0062         h_true_count->Fill(geta);
0063         if (fabs((px - gpx) / gpx) < MOMENTUM_MARGIN
0064             && fabs((py - gpy) / gpy) < MOMENTUM_MARGIN
0065             && fabs((pz - gpz) / gpz) < MOMENTUM_MARGIN)
0066             h_reco_count->Fill(geta);
0067
0068     }
0069
0070     const Long64_t nbins {h_true_count->GetSize() - 2}; /* -2 for underflow and overflow */
0071     Double_t *x {new Double_t[nbins]};
0072     Double_t *y {new Double_t[nbins]};
0073     Double_t *ey {new Double_t[nbins]};
0074     Long64_t top {0};
0075     for (Long64_t i {0};i < nbins; ++i) {
0076         if (h_true_count->GetBinContent(i + 1) != 0) {
0077             const Double_t n {h_reco_count->GetBinContent(i + 1)};
0078             const Double_t N {h_true_count->GetBinContent(i + 1)};
0079
0080             x[top] = {h_true_count->GetBinCenter(i + 1)};
0081             y[top] = {n / N};
0082             ey[top++] = {sqrt(n / (N * N) + (n * n) / (N * N * N))};
0083         }
0084     }
0085     TGraphErrors *gr {new TGraphErrors(top, x, y, nullptr, ey)};
0086     gr->SetMarkerColor(kBlue);
0087     gr->SetMarkerStyle(21);
0088     gr->SetMarkerSize(0.5);
0089     gr->GetXaxis()->SetTitle("#eta");
0090     gr->GetYaxis()->SetTitle("Efficiency");
0091
0092     TCanvas *count {new TCanvas("count", "FastTrack Event Count",
0093                 gStyle->GetCanvasDefW(), gStyle->GetCanvasDefH())};
0094     h_true_count->GetYaxis()->SetRangeUser(0, 2000);
0095     h_true_count->Draw();
0096     h_reco_count->Draw("SAME");
0097     TLegend *l1 {new TLegend(0.825, .9, .95, 0.8, "Track")};
0098     l1->SetTextSize(0.03);
0099     l1->AddEntry(h_true_count, "True", "l");
0100     l1->AddEntry(h_reco_count, "Reco", "l");
0101     l1->Draw();
0102     gPad->RedrawAxis();
0103
0104     TCanvas *efficiency {new TCanvas("efficiency", "FastTrack Efficiency",
0105                 gStyle->GetCanvasDefW(), gStyle->GetCanvasDefH())};
0106     gr->GetYaxis()->SetRangeUser(-0.05, 1);
0107     gr->Draw("ALP");
0108     TLegend *l2 {new TLegend(0.825, 0.9, 0.95, 0.8, "Track")};
0109     l2->SetTextSize(0.03);
0110     l2->AddEntry(gr, "Efficiency", "l");
0111     l2->Draw();
0112     gPad->RedrawAxis();
0113
0114     TImage *const img {TImage::Create()};
0115     img->FromPad(count);
0116     std::stringstream name;
0117     name << "FastTrack_Event_Count-" << "margin=" << MOMENTUM_MARGIN <<
0118         ".png";
0119     img->WriteImage(strdup(name.str().c_str()));
0120
0121     name.str("");
0122     img->FromPad(efficiency);
0123     name << "FastTrack_Efficiency-" << "margin=" << MOMENTUM_MARGIN <<
0124         ".png";
0125     img->WriteImage(strdup(name.str().c_str()));
0126     gApplication->Terminate(0);
0127 }
0128
0129 TTree *load_tree(const char *const file_name, const char *const tree_name)
0130 {
0131     return (TTree *) (new TFile(file_name, "READ"))->Get(tree_name);
0132 }
0133
0134 int main(int argc, char *argv[])
0135 {
0136     TApplication app {"FastTrack Efficiency Plots", &argc, argv};
0137     if (argc > 1) {
0138         std::stringstream tmp {argv[1]};
0139         tmp >> MOMENTUM_MARGIN;
0140     }
0141
0142     Plot_FastTrack_Efficiency();
0143     app.Run();
0144     return 0;
0145 }