tracking_performances_dis/analysis/vtx_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 <TF1.h>
0011 #include <TMath.h>
0012 #include <TStyle.h>
0013 #include <TLegend.h>
0014
0015 #include "fmt/color.h"
0016 #include "fmt/core.h"
0017
0018 #include "nlohmann/json.hpp"
0019
0020 double StudentT(double *x, double *par);
0021
0022 void vtx_dis_plots(const std::string& config_name)
0023 {
0024     // Read our configuration
0025     std::ifstream  config_file{config_name};
0026     nlohmann::json config;
0027     config_file >> config;
0028
0029     const std::string hists_file    = config["hists_file"];
0030     const std::string detector      = config["detector"];
0031     const std::string output_prefix = config["output_prefix"];
0032     const int         ebeam         = config["ebeam"];
0033     const int         pbeam         = config["pbeam"];
0034     const int         Q2_min        = config["Min_Q2"];
0035     const int         nfiles        = config["nfiles"];
0036
0037     fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0038                 "Plotting DIS tracking analysis...\n");
0039     fmt::print(" - Detector package: {}\n", detector);
0040     fmt::print(" - input file for histograms: {}\n", hists_file);
0041     fmt::print(" - output prefix for plots: {}\n", output_prefix);
0042     fmt::print(" - ebeam: {}\n", ebeam);
0043     fmt::print(" - pbeam: {}\n", pbeam);
0044     fmt::print(" - Minimum Q2: {}\n", Q2_min);
0045     fmt::print(" - nfiles: {}\n", nfiles);
0046
0047     //--------------------------------------------------------------------------------------------------------------------------------------------
0048
0049     // Read file with histograms
0050     TFile* file = new TFile(hists_file.c_str());
0051
0052     std::cout<<"Reading histograms..."<<std::endl;
0053
0054     TH2D* hgVr = (TH2D*) file->Get("recoVsMCTracks");
0055     TH2D* heff = (TH2D*) file->Get("recoVtxEff");
0056
0057     TH2D* hg1 = (TH2D*) file->Get("genVtxYvsXHist");
0058     TH2D* hg2 = (TH2D*) file->Get("genVtxRvsZHist");
0059
0060     TH2D* hr1 = (TH2D*) file->Get("recoVtxYvsX");
0061     TH2D* hr2 = (TH2D*) file->Get("recoVtxRvsZ");
0062
0063     TH2D* hres1r = (TH2D*) file->Get("vtxResXvsGenTrk");
0064     TH2D* hres2r = (TH2D*) file->Get("vtxResYvsGenTrk");
0065     TH2D* hres3r = (TH2D*) file->Get("vtxResZvsGenTrk");
0066
0067     TH2D* hres1g = (TH2D*) file->Get("vtxResXvsRecoTrk");
0068     TH2D* hres2g = (TH2D*) file->Get("vtxResYvsRecoTrk");
0069     TH2D* hres3g = (TH2D*) file->Get("vtxResZvsRecoTrk");
0070
0071     TH1D* hng1 = (TH1D*) file->Get("numGenTracks");
0072     TH1D* hng2 = (TH1D*) file->Get("numGenTrkswithVtx");
0073
0074     TH1D* hnr1 = (TH1D*) file->Get("numRecoTracks");
0075     TH1D* hnr2 = (TH1D*) file->Get("numRecoTrkswithVtx");
0076
0077     //--------------------------------------------------------------------------------------------------------------------------------------------
0078
0079     //Fit Function
0080     TF1 *myfunction = new TF1("fit", StudentT, -2, 2, 4);
0081
0082     //Fitting res v/s MC histograms
0083     myfunction->SetParameters(hres1g->GetMaximum(), 0, 0.05, 1);
0084     hres1g->FitSlicesY(myfunction, 0, -1, 10);
0085
0086     myfunction->SetParameters(hres2g->GetMaximum(), 0, 0.05, 1);
0087     hres2g->FitSlicesY(myfunction, 0, -1, 10);
0088
0089     myfunction->SetParameters(hres3g->GetMaximum(), 0, 0.05, 1);
0090     hres3g->FitSlicesY(myfunction, 0, -1, 10);
0091
0092     //Fitting res v/s RC histograms
0093     myfunction->SetParameters(hres1r->GetMaximum(), 0, 0.05, 1);
0094     hres1r->FitSlicesY(myfunction, 0, -1, 10);
0095
0096     myfunction->SetParameters(hres2r->GetMaximum(), 0, 0.05, 1);
0097     hres2r->FitSlicesY(myfunction, 0, -1, 10);
0098
0099     myfunction->SetParameters(hres3r->GetMaximum(), 0, 0.05, 1);
0100     hres3r->FitSlicesY(myfunction, 0, -1, 10);
0101
0102     //--------------------------------------------------------------------------------------------------------------------------------------------
0103
0104     // Make new efficiency histograms
0105     TH1D *vtxEffVsGenTrkHist = new TH1D("vtxEffVsGenTrk","",31,-0.5,30.5);
0106     TH1D *vtxEffVsRecoTrkHist = new TH1D("vtxEffVsRecoTrk","",31,-0.5,30.5);
0107
0108     // Fill efficiency v/s MC histogram
0109     for(int i=0; i<=hng1->GetNbinsX(); i++)
0110       {
0111     float neventsMC = hng1->GetBinContent(i);
0112     float nvtxevtsMC = hng2->GetBinContent(i);
0113
0114     if(neventsMC != 0)
0115     {
0116         vtxEffVsGenTrkHist->SetBinContent(i, nvtxevtsMC/neventsMC);
0117         vtxEffVsGenTrkHist->SetBinError(i, sqrt((nvtxevtsMC+1)/(neventsMC+2)*((nvtxevtsMC+2)/(neventsMC+3)-(nvtxevtsMC+1)/(neventsMC+2))));
0118     }
0119       }
0120
0121     vtxEffVsGenTrkHist->SetMarkerColor(kRed);
0122     vtxEffVsGenTrkHist->SetMarkerStyle(8); vtxEffVsGenTrkHist->SetMarkerSize(1.2);
0123     vtxEffVsGenTrkHist->SetTitle("Vertexing Efficiency vs MC Tracks");
0124     vtxEffVsGenTrkHist->GetXaxis()->SetTitle("N_{MC}");
0125     vtxEffVsGenTrkHist->GetYaxis()->SetTitle("Vertexing Efficiency");
0126     vtxEffVsGenTrkHist->GetYaxis()->SetRangeUser(0, 1.2);
0127
0128     //Fill efficiency v/s RC histogram
0129     for(int i=0; i<=hnr1->GetNbinsX(); i++)
0130       {
0131     float neventsRC = hnr1->GetBinContent(i);
0132     float nvtxevtsRC = hnr2->GetBinContent(i);
0133
0134     if(neventsRC != 0)
0135     {
0136         vtxEffVsRecoTrkHist->SetBinContent(i, nvtxevtsRC/neventsRC);
0137         vtxEffVsRecoTrkHist->SetBinError(i, sqrt((nvtxevtsRC+1)/(neventsRC+2)*((nvtxevtsRC+2)/(neventsRC+3)-(nvtxevtsRC+1)/(neventsRC+2))));
0138     }
0139       }
0140
0141     vtxEffVsRecoTrkHist->SetMarkerColor(kRed);
0142     vtxEffVsRecoTrkHist->SetMarkerStyle(8); vtxEffVsRecoTrkHist->SetMarkerSize(1.2);
0143     vtxEffVsRecoTrkHist->SetTitle("Vertexing Efficiency vs RC Tracks");
0144     vtxEffVsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}");
0145     vtxEffVsRecoTrkHist->GetYaxis()->SetTitle("Vertexing Efficiency");
0146     vtxEffVsRecoTrkHist->GetYaxis()->SetRangeUser(0, 1.2);
0147
0148     //--------------------------------------------------------------------------------------------------------------------------------------------
0149
0150     // Make plots and save to PDF file
0151
0152     // Update Style
0153     gStyle->SetPadGridX(0);
0154     gStyle->SetPadGridY(0);
0155     gStyle->SetOptStat(0);
0156
0157     std::cout<<"Making plots..."<<std::endl;
0158
0159     // MC vs RC associated particles
0160     TCanvas *c1 = new TCanvas("c1","MC vs Reco Tracks",800,600);
0161     c1->cd(1);
0162     auto func = new TF1("func","x",-1,40);
0163     func->SetLineStyle(2); func->SetLineColor(1);
0164
0165     hgVr->Draw();
0166     func->Draw("SAMEL");
0167
0168     //Vertexing Efficiency
0169     TCanvas *c2 = new TCanvas("c2","Vertexing Efficiency",800,600);
0170     int nevents(0);
0171     for(int i=0; i<=heff->GetNbinsX(); i++) nevents = nevents + heff->GetBinContent(i);
0172     heff->Scale(100./nevents);
0173
0174     c2->cd(1);
0175     heff->Draw("pe");
0176
0177     // MC vx versus vy Vertex
0178     TCanvas *c3 = new TCanvas("c3","MC Vertices",800,600);
0179     c3->cd(1);
0180     hg1->Draw();
0181
0182     // MC vr versus vz Vertex
0183     TCanvas *c4 = new TCanvas("c4","MC Vertices",800,600);
0184     c4->cd(1);
0185     hg2->Draw();
0186
0187     // Reconstructed Vertex vx versus vy
0188     TCanvas *c5 = new TCanvas("c5","Reconstructed Vertices",800,600);
0189     c5->cd(1);
0190     hr1->Draw();
0191
0192     // Reconstructed Vertex vr versus vz
0193     TCanvas *c6 = new TCanvas("c6","Reconstructed Vertices",800,600);
0194     c6->cd(1);
0195     hr2->Draw();
0196
0197     //Vertex Resolution vs MC Tracks
0198     TCanvas *c7 = new TCanvas("c7","VtxResX vs MCTrks",800,600);
0199     c7->cd(1);
0200     hres1g->Draw("COLZ");
0201
0202     TCanvas *c8 = new TCanvas("c8","VtxResY vs MCTrks",800,600);
0203     c8->cd(1);
0204     hres2g->Draw("COLZ");
0205
0206     TCanvas *c9 = new TCanvas("c9","VtxResZ vs MCTrks",800,600);
0207     c9->cd(1);
0208     hres3g->Draw("COLZ");
0209
0210     //Vertex Resolution vs RC Tracks
0211     TCanvas *c10 = new TCanvas("c10","VtxResX vs RCTrks",800,600);
0212     c10->cd(1);
0213     hres1r->Draw("COLZ");
0214
0215     TCanvas *c11 = new TCanvas("c11","VtxResY vs RCTrks",800,600);
0216     c11->cd(1);
0217     hres2r->Draw("COLZ");
0218
0219     TCanvas *c12 = new TCanvas("c12","VtxResZ vs RCTrks",800,600);
0220     c12->cd(1);
0221     hres3r->Draw("COLZ");
0222
0223     // Res Sigma v/s MC tracks
0224     TCanvas *c13 = new TCanvas("c13","Vertex Resolution vs MC Tracks",800,600);
0225     c13->cd(1);
0226
0227     TH1D *resXsigmaM = (TH1D*)gDirectory->Get("vtxResXvsGenTrk_2");
0228     TH1D *resYsigmaM = (TH1D*)gDirectory->Get("vtxResYvsGenTrk_2");
0229     TH1D *resZsigmaM = (TH1D*)gDirectory->Get("vtxResZvsGenTrk_2");
0230
0231     resXsigmaM->SetMarkerStyle(20); resYsigmaM->SetMarkerStyle(21); resZsigmaM->SetMarkerStyle(22);
0232     resXsigmaM->SetMarkerSize(1.2); resYsigmaM->SetMarkerSize(1.2); resZsigmaM->SetMarkerSize(1.2);
0233     resXsigmaM->SetMarkerColor(kBlue); resYsigmaM->SetMarkerColor(kRed); resZsigmaM->SetMarkerColor(kBlack);
0234     resXsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks"); resYsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks");
0235     resZsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks");
0236     resZsigmaM->GetXaxis()->SetTitle("N_{MC}");
0237     resXsigmaM->GetYaxis()->SetTitle("#sigma (mm)"); resYsigmaM->GetYaxis()->SetTitle("#sigma (mm)"); resZsigmaM->GetYaxis()->SetTitle("#sigma (mm)");
0238     resXsigmaM->GetYaxis()->SetRangeUser(0, 1); resYsigmaM->GetYaxis()->SetRangeUser(0, 1); resZsigmaM->GetYaxis()->SetRangeUser(0, 1);
0239
0240     resXsigmaM->Draw("P");
0241     resYsigmaM->Draw("PSAME");
0242     resZsigmaM->Draw("PSAME");
0243
0244     TLegend *legend1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0245     legend1->AddEntry(resXsigmaM, "v_{x}", "lep");
0246     legend1->AddEntry(resYsigmaM, "v_{y}", "lep");
0247     legend1->AddEntry(resZsigmaM, "v_{z}", "lep");
0248     legend1->Draw();
0249
0250     // Res Mean v/s MC tracks
0251     TCanvas *c14 = new TCanvas("c14","Vertex Resolution Mean vs MC Tracks",800,600);
0252     c14->cd(1);
0253
0254     TH1D *resXmeanM = (TH1D*)gDirectory->Get("vtxResXvsGenTrk_1");
0255     TH1D *resYmeanM = (TH1D*)gDirectory->Get("vtxResYvsGenTrk_1");
0256     TH1D *resZmeanM = (TH1D*)gDirectory->Get("vtxResZvsGenTrk_1");
0257
0258     resXmeanM->SetMarkerStyle(20); resYmeanM->SetMarkerStyle(21); resZmeanM->SetMarkerStyle(22);
0259     resXmeanM->SetMarkerSize(1.2); resYmeanM->SetMarkerSize(1.2); resZmeanM->SetMarkerSize(1.2);
0260     resXmeanM->SetMarkerColor(kBlue); resYmeanM->SetMarkerColor(kRed); resZmeanM->SetMarkerColor(kBlack);
0261     resXmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks"); resYmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks");
0262     resZmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks");
0263     resZmeanM->GetXaxis()->SetTitle("N_{MC}");
0264     resXmeanM->GetYaxis()->SetTitle("#mu (mm)"); resYmeanM->GetYaxis()->SetTitle("#mu (mm)"); resZmeanM->GetYaxis()->SetTitle("#mu (mm)");
0265     resXmeanM->GetYaxis()->SetRangeUser(-1, 1); resYmeanM->GetYaxis()->SetRangeUser(-1, 1); resZmeanM->GetYaxis()->SetRangeUser(-1, 1);
0266
0267     resXmeanM->Draw("P");
0268     resYmeanM->Draw("PSAME");
0269     resZmeanM->Draw("PSAME");
0270
0271     TLegend *legend2 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0272     legend2->AddEntry(resXmeanM, "v_{x}", "lep");
0273     legend2->AddEntry(resYmeanM, "v_{y}", "lep");
0274     legend2->AddEntry(resZmeanM, "v_{z}", "lep");
0275     legend2->Draw();
0276
0277     //Res Sigma v/s RC tracks
0278     TCanvas *c15 = new TCanvas("c15","Vertex Resolution vs RC Tracks",800,600);
0279     c15->cd(1);
0280
0281     TH1D *resXsigmaR = (TH1D*)gDirectory->Get("vtxResXvsRecoTrk_2");
0282     TH1D *resYsigmaR = (TH1D*)gDirectory->Get("vtxResYvsRecoTrk_2");
0283     TH1D *resZsigmaR = (TH1D*)gDirectory->Get("vtxResZvsRecoTrk_2");
0284
0285     resXsigmaR->SetMarkerStyle(20); resYsigmaR->SetMarkerStyle(21); resZsigmaR->SetMarkerStyle(22);
0286     resXsigmaR->SetMarkerSize(1.2); resYsigmaR->SetMarkerSize(1.2); resZsigmaR->SetMarkerSize(1.2);
0287     resXsigmaR->SetMarkerColor(kBlue); resYsigmaR->SetMarkerColor(kRed); resZsigmaR->SetMarkerColor(kBlack);
0288     resXsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks"); resYsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks");
0289     resZsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks");
0290     resZsigmaR->GetXaxis()->SetTitle("N_{RC}");
0291     resXsigmaR->GetYaxis()->SetTitle("#sigma (mm)"); resYsigmaR->GetYaxis()->SetTitle("#sigma (mm)"); resZsigmaR->GetYaxis()->SetTitle("#sigma (mm)");
0292     resXsigmaR->GetYaxis()->SetRangeUser(0, 1); resYsigmaR->GetYaxis()->SetRangeUser(0, 1); resZsigmaR->GetYaxis()->SetRangeUser(0, 1);
0293
0294     resXsigmaR->Draw("P");
0295     resYsigmaR->Draw("PSAME");
0296     resZsigmaR->Draw("PSAME");
0297
0298     TLegend *legend3 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0299     legend3->AddEntry(resXsigmaR, "v_{x}", "lep");
0300     legend3->AddEntry(resYsigmaR, "v_{y}", "lep");
0301     legend3->AddEntry(resZsigmaR, "v_{z}", "lep");
0302     legend3->Draw();
0303
0304     //Res Mean v/s RC tracks
0305     TCanvas *c16 = new TCanvas("c16","Vertex Resolution Mean vs RC Tracks",800,600);
0306     c16->cd(1);
0307
0308     TH1D *resXmeanR = (TH1D*)gDirectory->Get("vtxResXvsRecoTrk_1");
0309     TH1D *resYmeanR = (TH1D*)gDirectory->Get("vtxResYvsRecoTrk_1");
0310     TH1D *resZmeanR = (TH1D*)gDirectory->Get("vtxResZvsRecoTrk_1");
0311
0312     resXmeanR->SetMarkerStyle(20); resYmeanR->SetMarkerStyle(21); resZmeanR->SetMarkerStyle(22);
0313     resXmeanR->SetMarkerSize(1.2); resYmeanR->SetMarkerSize(1.2); resZmeanR->SetMarkerSize(1.2);
0314     resXmeanR->SetMarkerColor(kBlue); resYmeanR->SetMarkerColor(kRed); resZmeanR->SetMarkerColor(kBlack);
0315     resXmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks"); resYmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks");
0316     resZmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks");
0317     resZmeanR->GetXaxis()->SetTitle("N_{RC}");
0318     resXmeanR->GetYaxis()->SetTitle("#mu (mm)"); resYmeanR->GetYaxis()->SetTitle("#mu (mm)"); resZmeanR->GetYaxis()->SetTitle("#mu (mm)");
0319     resXmeanR->GetYaxis()->SetRangeUser(-1, 1); resYmeanR->GetYaxis()->SetRangeUser(-1, 1); resZmeanR->GetYaxis()->SetRangeUser(-1, 1);
0320
0321     resXmeanR->Draw("P");
0322     resYmeanR->Draw("PSAME");
0323     resZmeanR->Draw("PSAME");
0324
0325     TLegend *legend4 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0326     legend4->AddEntry(resXmeanR, "v_{x}", "lep");
0327     legend4->AddEntry(resYmeanR, "v_{y}", "lep");
0328     legend4->AddEntry(resZmeanR, "v_{z}", "lep");
0329     legend4->Draw();
0330
0331     //Vertexing Efficiency vs MC Tracks
0332     TCanvas *c17 = new TCanvas("c17","Vertexing Efficiency vs MC Tracks",800,600);
0333     c17->cd(1);
0334
0335     vtxEffVsGenTrkHist->Draw("P");
0336
0337     //Vertexing Efficiency vs RC Tracks
0338     TCanvas *c18 = new TCanvas("c18","Vertexing Efficiency vs RC Tracks",800,600);
0339     c18->cd(1);
0340
0341     vtxEffVsRecoTrkHist->Draw("P");
0342
0343     //--------------------------------------------------------------------------------------------------------------------------------------------
0344
0345     // Print plots to pdf file
0346     c1->Print(fmt::format("{}.pdf[", output_prefix).c_str());
0347     c1->Print(fmt::format("{}.pdf", output_prefix).c_str());
0348     c2->Print(fmt::format("{}.pdf", output_prefix).c_str());
0349     c3->Print(fmt::format("{}.pdf", output_prefix).c_str());
0350     c4->Print(fmt::format("{}.pdf", output_prefix).c_str());
0351     c5->Print(fmt::format("{}.pdf", output_prefix).c_str());
0352     c6->Print(fmt::format("{}.pdf", output_prefix).c_str());
0353     c7->Print(fmt::format("{}.pdf", output_prefix).c_str());
0354     c8->Print(fmt::format("{}.pdf", output_prefix).c_str());
0355     c9->Print(fmt::format("{}.pdf", output_prefix).c_str());
0356     c10->Print(fmt::format("{}.pdf", output_prefix).c_str());
0357     c11->Print(fmt::format("{}.pdf", output_prefix).c_str());
0358     c12->Print(fmt::format("{}.pdf", output_prefix).c_str());
0359     c13->Print(fmt::format("{}.pdf", output_prefix).c_str());
0360     c14->Print(fmt::format("{}.pdf", output_prefix).c_str());
0361     c15->Print(fmt::format("{}.pdf", output_prefix).c_str());
0362     c16->Print(fmt::format("{}.pdf", output_prefix).c_str());
0363     c17->Print(fmt::format("{}.pdf", output_prefix).c_str());
0364     c18->Print(fmt::format("{}.pdf", output_prefix).c_str());
0365     c18->Print(fmt::format("{}.pdf]", output_prefix).c_str());
0366
0367 }
0368
0369 //Defining Fitting Function
0370 double StudentT(double *x, double *par){
0371     double norm = par[0];
0372     double mean = par[1];
0373     double sigma = par[2];
0374     double nu = par[3];
0375
0376     double pi = 3.14;
0377     double st = norm * (TMath::Gamma((nu+1.0)/2.0)/(TMath::Gamma(nu/2.0)*TMath::Sqrt(pi*nu)*sigma)) * TMath::Power( (1.0+TMath::Power((x[0]-mean)/sigma,2.0)/nu), (-(nu+1.0)/2.0) );
0378     return st;
0379 }