tracking_performances_dis/analysis/trk_dis_plots.cxx

0001 #include <fstream>
0002 #include <iostream>
0003 #include <string>
0004
0005 #include <TCanvas.h>
0006 #include <TChain.h>
0007 #include <TFile.h>
0008 #include <TH1D.h>
0009 #include <TH2D.h>
0010 #include <TMath.h>
0011 #include <TStyle.h>
0012 #include <TLegend.h>
0013
0014 #include "fmt/color.h"
0015 #include "fmt/core.h"
0016
0017 #include "nlohmann/json.hpp"
0018
0019 void trk_dis_plots(const std::string& config_name)
0020 {
0021     // Read our configuration

0022     std::ifstream  config_file{config_name};
0023     nlohmann::json config;
0024     config_file >> config;
0025
0026     const std::string hists_file    = config["hists_file"];
0027     const std::string detector      = config["detector"];
0028     const std::string output_prefix = config["output_prefix"];
0029     const int         ebeam         = config["ebeam"];
0030     const int         pbeam         = config["pbeam"];
0031     const int         Q2_min        = config["Min_Q2"];
0032     const int         nfiles        = config["nfiles"];
0033
0034     fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0035                 "Plotting DIS tracking analysis...\n");
0036     fmt::print(" - Detector package: {}\n", detector);
0037     fmt::print(" - input file for histograms: {}\n", hists_file);
0038     fmt::print(" - output prefix for plots: {}\n", output_prefix);
0039     fmt::print(" - ebeam: {}\n", ebeam);
0040     fmt::print(" - pbeam: {}\n", pbeam);
0041     fmt::print(" - Minimum Q2: {}\n", Q2_min);
0042     fmt::print(" - nfiles: {}\n", nfiles);
0043
0044     //--------------------------------------------------------------------------------------------------------------------------------------------

0045
0046     // Read file with histograms

0047     TFile* file = new TFile(hists_file.c_str());
0048
0049     std::cout<<"Reading histograms..."<<std::endl;
0050
0051     TH1D* h1a = (TH1D*) file->Get("h1a");
0052     TH1D* h1a1 = (TH1D*) file->Get("h1a1");
0053     TH1D* h1a2 = (TH1D*) file->Get("h1a2");
0054
0055     TH1D* h1b = (TH1D*) file->Get("h1b");
0056     TH1D* h1b1 = (TH1D*) file->Get("h1b1");
0057     TH1D* h1b2 = (TH1D*) file->Get("h1b2");
0058
0059     TH1D* h1c = (TH1D*) file->Get("h1c");
0060     TH1D* h1c1 = (TH1D*) file->Get("h1c1");
0061     TH1D* h1c2 = (TH1D*) file->Get("h1c2");
0062
0063     TH1D* h2a = (TH1D*) file->Get("h2a");
0064     TH1D* h2b = (TH1D*) file->Get("h2b");
0065
0066     //--------------------------------------------------------------------------------------------------------------------------------------------

0067
0068     // Make ratio histograms

0069     TH1 *h1rb1 = new TH1D("h1rb1","",100,-4,4); //Real-seeded tracks (Pt > 200 MeV/c cut)

0070     TH1 *h1rc1 = new TH1D("h1rc1","",100,-4,4); //Truth-seeded tracks (Pt > 200 MeV/c cut)

0071     TH1 *h1rb2 = new TH1D("h1rb2","",100,-4,4); //Real-seeded tracks (Pt > 500 MeV/c cut)

0072     TH1 *h1rc2 = new TH1D("h1rc2","",100,-4,4); //Truth-seeded tracks (Pt > 500 MeV/c cut)

0073
0074     h1rb1 = (TH1*) h1b1->Clone("h1rb1");
0075     h1rb1->Divide(h1a1);
0076     h1rb1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
0077     h1rb1->GetXaxis()->SetTitle("#eta");h1rb1->GetXaxis()->CenterTitle();
0078     h1rb1->GetYaxis()->SetTitle("Ratio");h1rb1->GetYaxis()->CenterTitle();
0079
0080     h1rc1 = (TH1*) h1c1->Clone("h1rc1");
0081     h1rc1->Divide(h1a1);
0082     h1rc1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
0083     h1rc1->GetXaxis()->SetTitle("#eta");h1rc1->GetXaxis()->CenterTitle();
0084     h1rc1->GetYaxis()->SetTitle("Ratio");h1rc1->GetYaxis()->CenterTitle();
0085
0086     h1rb2 = (TH1*) h1b2->Clone("h1rb2");
0087     h1rb2->Divide(h1a2);
0088     h1rb2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
0089     h1rb2->GetXaxis()->SetTitle("#eta");h1rb2->GetXaxis()->CenterTitle();
0090     h1rb2->GetYaxis()->SetTitle("Ratio");h1rb2->GetYaxis()->CenterTitle();
0091
0092     h1rc2 = (TH1*) h1c2->Clone("h1rc2");
0093     h1rc2->Divide(h1a2);
0094     h1rc2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
0095     h1rc2->GetXaxis()->SetTitle("#eta");h1rc2->GetXaxis()->CenterTitle();
0096     h1rc2->GetYaxis()->SetTitle("Ratio");h1rc2->GetYaxis()->CenterTitle();
0097
0098     //--------------------------------------------------------------------------------------------------------------------------------------------

0099
0100     // Make plots and save to PDF file

0101
0102     // Update Style

0103     gStyle->SetPadGridX(0);
0104     gStyle->SetPadGridY(0);
0105     gStyle->SetOptStat(0);
0106
0107     std::cout<<"Making plots..."<<std::endl;
0108
0109     //Generated charged particles

0110     TCanvas *c1a = new TCanvas("c1a");
0111     h1a->Draw();
0112     h1a1->Draw("same");
0113     h1a2->Draw("same");
0114
0115     TLegend *leg1a = new TLegend(0.25,0.6,0.6,0.875);
0116     leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
0117     leg1a->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0118     leg1a->AddEntry("h1a","All generated charged particles","l");
0119     leg1a->AddEntry("h1a1","+ P_{T} > 200 MeV/c","l");
0120     leg1a->AddEntry("h1a2","+ P_{T} > 500 MeV/c","l");
0121     leg1a->Draw();
0122
0123     //Real-seeded tracks

0124     TCanvas *c1b = new TCanvas("c1b");
0125     h1b->Draw();
0126     h1b1->Draw("same");
0127     h1b2->Draw("same");
0128
0129     TLegend *leg1b = new TLegend(0.25,0.6,0.6,0.875);
0130     leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
0131     leg1b->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0132     leg1b->AddEntry("h1b","All real-seeded tracks","l");
0133     leg1b->AddEntry("h1b1","+ P_{T} > 200 MeV/c","l");
0134     leg1b->AddEntry("h1b2","+ P_{T} > 500 MeV/c","l");
0135     leg1b->Draw();
0136
0137     //Truth-seeded tracks

0138     TCanvas *c1c = new TCanvas("c1c");
0139     h1c->Draw();
0140     h1c1->Draw("same");
0141     h1c2->Draw("same");
0142
0143     TLegend *leg1c = new TLegend(0.25,0.6,0.6,0.875);
0144     leg1c->SetBorderSize(0);leg1c->SetFillStyle(0);
0145     leg1c->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0146     leg1c->AddEntry("h1c","All truth-seeded tracks","l");
0147     leg1c->AddEntry("h1c1","+ P_{T} > 200 MeV/c","l");
0148     leg1c->AddEntry("h1c2","+ P_{T} > 500 MeV/c","l");
0149     leg1c->Draw();
0150
0151     //Comparison 1

0152     TCanvas *c1d = new TCanvas("c1d");
0153     auto frame_d1 = c1d->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
0154     frame_d1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_d1->GetXaxis()->CenterTitle();
0155     h1a1->Draw("same");
0156     h1b1->Draw("same");
0157     h1c1->Draw("same");
0158
0159     TLegend *leg1d = new TLegend(0.15,0.675,0.6,0.875);
0160     leg1d->SetBorderSize(0);leg1d->SetFillStyle(0);
0161     leg1d->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0162     leg1d->AddEntry("h1a1","Generated charged particles w/ P_{T} > 200 MeV/c","l");
0163     leg1d->AddEntry("h1b1","Real-seeded tracks w/ P_{T} > 200 MeV/c","l");
0164     leg1d->AddEntry("h1c1","Truth-seeded tracks w/ P_{T} > 200 MeV/c","l");
0165     leg1d->Draw();
0166
0167     //Comparison 2a

0168     TCanvas *c1e = new TCanvas("c1e");
0169     auto frame_e1 = c1e->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
0170     frame_e1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_e1->GetXaxis()->CenterTitle();
0171     h1a2->Draw("same");
0172     h1b2->Draw("P same");
0173
0174     TLegend *leg1e = new TLegend(0.15,0.675,0.6,0.875);
0175     leg1e->SetBorderSize(0);leg1e->SetFillStyle(0);
0176     leg1e->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0177     leg1e->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
0178     leg1e->AddEntry("h1b2","Real-seeded tracks w/ P_{T} > 500 MeV/c","p");
0179     leg1e->Draw();
0180
0181     //Comparison 2b

0182     TCanvas *c1e1 = new TCanvas("c1e1");
0183     frame_e1->Draw();
0184     h1a2->Draw("same");
0185     h1c2->Draw("P same");
0186
0187     TLegend *leg1e1 = new TLegend(0.15,0.675,0.6,0.875);
0188     leg1e1->SetBorderSize(0);leg1e1->SetFillStyle(0);
0189     leg1e1->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0190     leg1e1->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
0191     leg1e1->AddEntry("h1c2","Truth-seeded tracks w/ P_{T} > 500 MeV/c","p");
0192     leg1e1->Draw();
0193
0194     //Comparison 1 ratio

0195     TCanvas *c1f = new TCanvas("c1f");
0196     h1rb1->Draw();
0197     h1rc1->Draw("same");
0198
0199     TLegend *leg1f = new TLegend(0.575,0.25,0.875,0.45);
0200     leg1f->SetBorderSize(0);leg1f->SetFillStyle(0);
0201     leg1f->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0202     leg1f->AddEntry("h1rb1","Real-seeded tracking","l");
0203     leg1f->AddEntry("h1rc1","Truth-seeded tracking","l");
0204     leg1f->Draw();
0205
0206     //Comparison 2 ratio

0207     TCanvas *c1g = new TCanvas("c1g");
0208     h1rb2->Draw();
0209     h1rc2->Draw("same");
0210
0211     TLegend *leg1g = new TLegend(0.575,0.25,0.875,0.45);
0212     leg1g->SetBorderSize(0);leg1g->SetFillStyle(0);
0213     leg1g->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0214     leg1g->AddEntry("h1rb2","Real-seeded tracking","l");
0215     leg1g->AddEntry("h1rc2","Truth-seeded tracking","l");
0216     leg1g->Draw();
0217
0218     //Hit-based associations -- real-seeded tracks

0219     TCanvas *c2a = new TCanvas("c2a");
0220     c2a->SetLogy();
0221     h2a->Draw();
0222
0223     TLegend *leg2a = new TLegend(0.25,0.6,0.6,0.875);
0224     leg2a->SetBorderSize(0);leg2a->SetFillStyle(0);
0225     leg2a->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0226     leg2a->Draw();
0227
0228     //Hit-based associations -- truth-seeded tracks

0229     TCanvas *c2b = new TCanvas("c2b");
0230     c2b->SetLogy();
0231     h2b->Draw();
0232
0233     TLegend *leg2b = new TLegend(0.25,0.6,0.6,0.875);
0234     leg2b->SetBorderSize(0);leg2a->SetFillStyle(0);
0235     leg2b->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
0236     leg2b->Draw();
0237
0238     //--------------------------------------------------------------------------------------------------------------------------------------------

0239
0240     // Print plots to pdf file

0241     c1a->Print(fmt::format("{}.pdf[", output_prefix).c_str());
0242     c1a->Print(fmt::format("{}.pdf", output_prefix).c_str());
0243     c1b->Print(fmt::format("{}.pdf", output_prefix).c_str());
0244     c1c->Print(fmt::format("{}.pdf", output_prefix).c_str());
0245     c1d->Print(fmt::format("{}.pdf", output_prefix).c_str());
0246     c1e->Print(fmt::format("{}.pdf", output_prefix).c_str());
0247     c1e1->Print(fmt::format("{}.pdf", output_prefix).c_str());
0248     c1f->Print(fmt::format("{}.pdf", output_prefix).c_str());
0249     c1g->Print(fmt::format("{}.pdf", output_prefix).c_str());
0250     c2a->Print(fmt::format("{}.pdf", output_prefix).c_str());
0251     c2b->Print(fmt::format("{}.pdf", output_prefix).c_str());
0252     c2b->Print(fmt::format("{}.pdf]", output_prefix).c_str());
0253
0254 }