diff --git a/benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx b/benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx new file mode 100644 index 00000000..9493db67 --- /dev/null +++ b/benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx @@ -0,0 +1,307 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fmt/color.h" +#include "fmt/core.h" + +#include "nlohmann/json.hpp" + +void vtx_dis_analysis(const std::string& config_name) +{ + // Read our configuration + std::ifstream config_file{config_name}; + nlohmann::json config; + config_file >> config; + + const std::string rec_file = config["rec_file"]; + const std::string detector = config["detector"]; + const std::string output_prefix = config["output_prefix"]; + const int ebeam = config["ebeam"]; + const int pbeam = config["pbeam"]; + const int Q2_min = config["Min_Q2"]; + + fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), + "Running DIS tracking analysis...\n"); + fmt::print(" - Detector package: {}\n", detector); + fmt::print(" - input file: {}\n", rec_file); + fmt::print(" - output prefix for histograms: {}\n", output_prefix); + fmt::print(" - ebeam: {}\n", ebeam); + fmt::print(" - pbeam: {}\n", pbeam); + fmt::print(" - Minimum Q2: {}\n", Q2_min); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Set output file for the histograms + std::string output_name_hists = fmt::format("{}.root", output_prefix); + cout << "Output file for histograms = " << output_name_hists << endl; + TFile* ofile = new TFile(output_name_hists.c_str(), "RECREATE"); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Set up input file chain + TChain *mychain = new TChain("events"); + mychain->Add(rec_file.c_str()); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // TTreeReader + TTreeReader tree_reader(mychain); + + // Reco Vertex + TTreeReaderArray recoVtxType = {tree_reader, "CentralTrackVertices.type"}; + TTreeReaderArray recoVtxX = {tree_reader, "CentralTrackVertices.position.x"}; + TTreeReaderArray recoVtxY = {tree_reader, "CentralTrackVertices.position.y"}; + TTreeReaderArray recoVtxZ = {tree_reader, "CentralTrackVertices.position.z"}; + + TTreeReaderArray assoPartBegin = {tree_reader, "CentralTrackVertices.associatedParticles_begin"}; + TTreeReaderArray assoPartEnd = {tree_reader, "CentralTrackVertices.associatedParticles_end"}; + TTreeReaderArray assoPartIndex = {tree_reader, "_CentralTrackVertices_associatedParticles.index"}; + + // MC + TTreeReaderArray mcGenStat = {tree_reader, "MCParticles.generatorStatus"}; + TTreeReaderArray mcPDG = {tree_reader, "MCParticles.PDG"}; + TTreeReaderArray mcCharge = {tree_reader, "MCParticles.charge"}; + TTreeReaderArray mcMomX = {tree_reader, "MCParticles.momentum.x"}; + TTreeReaderArray mcMomY = {tree_reader, "MCParticles.momentum.y"}; + TTreeReaderArray mcMomZ = {tree_reader, "MCParticles.momentum.z"}; + + TTreeReaderArray mcVtxX = {tree_reader, "MCParticles.vertex.x"}; + TTreeReaderArray mcVtxY = {tree_reader, "MCParticles.vertex.y"}; + TTreeReaderArray mcVtxZ = {tree_reader, "MCParticles.vertex.z"}; + + TTreeReaderArray mcParentBegin = {tree_reader, "MCParticles.parents_begin"}; + TTreeReaderArray mcParentEnd = {tree_reader, "MCParticles.parents_end"}; + TTreeReaderArray mcParentIndex = {tree_reader, "_MCParticles_parents.index"}; + + // Reco + TTreeReaderArray recoType = {tree_reader, "ReconstructedChargedParticles.type"}; + + + //-------------------------------------------------------------------------------------------------------------------------------------------- + // Define Histograms + + // Gen + TH1D *numGenTracksHist = new TH1D("numGenTracks","",31,-0.5,30.5); + + TH2D *genVtxYvsXHist = new TH2D("genVtxYvsXHist","",200,-1.,1.,200,-1,1); + genVtxYvsXHist->Draw("COLZ"); + genVtxYvsXHist->SetTitle("MC Vertex: v_{x} versus v_{y}"); + genVtxYvsXHist->GetXaxis()->SetTitle("x-coordinate (in mm)"); + genVtxYvsXHist->GetYaxis()->SetTitle("y-coordinate (in mm)"); + genVtxYvsXHist->GetXaxis()->CenterTitle(); genVtxYvsXHist->GetYaxis()->CenterTitle(); + + TH2D *genVtxRvsZHist = new TH2D("genVtxRvsZHist","",200,-100.,100.,100,0,0.5); + genVtxRvsZHist->Draw("COLZ"); + genVtxRvsZHist->SetTitle("MC Vertex: v_{r} versus v_{z}"); + genVtxRvsZHist->GetXaxis()->SetTitle("z-coordinate (in mm)"); + genVtxRvsZHist->GetYaxis()->SetTitle("#sqrt{x^{2} + y^{2}} (in mm)"); + genVtxRvsZHist->GetXaxis()->CenterTitle(); genVtxRvsZHist->GetYaxis()->CenterTitle(); + + // Reco + TH1D *numRecoTracksHist = new TH1D("numRecoTracks","",31,-0.5,30.5); + + TH1D *recoVtxEffHist = new TH1D("recoVtxEff","",4,-0.5,3.5); + recoVtxEffHist->Draw("P"); + recoVtxEffHist->SetLineColor(kRed); recoVtxEffHist->SetMarkerColor(kRed); + recoVtxEffHist->SetMarkerStyle(8); recoVtxEffHist->SetMarkerSize(1.2); + recoVtxEffHist->SetTitle("Vertexing Efficiency"); + recoVtxEffHist->GetXaxis()->SetTitle("Vertex Count"); recoVtxEffHist->GetYaxis()->SetTitle("nEvents/total_events %"); + recoVtxEffHist->GetXaxis()->CenterTitle(); recoVtxEffHist->GetYaxis()->CenterTitle(); + + TH2D *recoVtxYvsXHist = new TH2D("recoVtxYvsX","",200,-1.,1.,200,-1,1); + recoVtxYvsXHist->Draw("COLZ"); + recoVtxYvsXHist->SetTitle("Reconstructed Vertex: v_{x} versus v_{y}"); + recoVtxYvsXHist->GetXaxis()->SetTitle("x-coordinate (in mm)"); recoVtxYvsXHist->GetYaxis()->SetTitle("y-coordinate (in mm)"); + recoVtxYvsXHist->GetXaxis()->CenterTitle(); recoVtxYvsXHist->GetYaxis()->CenterTitle(); + + TH2D *recoVtxRvsZHist = new TH2D("recoVtxRvsZ","",200,-100.,100.,100,0.,0.8); + recoVtxRvsZHist->Draw("COLZ"); + recoVtxRvsZHist->SetTitle("Reconstructed Vertex: v_{r} versus v_{z}"); + recoVtxRvsZHist->GetXaxis()->SetTitle("z-coordinate (in mm)"); + recoVtxRvsZHist->GetYaxis()->SetTitle("#sqrt{x^{2} + y^{2}} (in mm)"); + recoVtxRvsZHist->GetXaxis()->CenterTitle(); recoVtxRvsZHist->GetYaxis()->CenterTitle(); + + //Reco Vs MC + TH2D *recoVsMCTracksHist = new TH2D("recoVsMCTracks","",31,-0.5,30.5,31,-0.5,30.5); + recoVsMCTracksHist->Draw("COLZ"); + recoVsMCTracksHist->SetTitle("Number of particles associated with vertex"); + recoVsMCTracksHist->GetXaxis()->SetTitle("N_{MC}"); recoVsMCTracksHist->GetYaxis()->SetTitle("N_{RC}"); + recoVsMCTracksHist->GetXaxis()->CenterTitle(); recoVsMCTracksHist->GetYaxis()->CenterTitle(); + + // Resolution + TH2D *vtxResXvsGenTrkHist = new TH2D("vtxResXvsGenTrk","",31,-0.5,30.5,200,-1,1); + vtxResXvsGenTrkHist->Draw("COLZ"); + vtxResXvsGenTrkHist->SetTitle("Vertex Resolution X vs MC Tracks"); + vtxResXvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}"); vtxResXvsGenTrkHist->GetYaxis()->SetTitle("recVtx_x - mcVtx_x (in mm)"); + vtxResXvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResXvsGenTrkHist->GetYaxis()->CenterTitle(); + + TH2D *vtxResYvsGenTrkHist = new TH2D("vtxResYvsGenTrk","",31,-0.5,30.5,200,-1,1); + vtxResYvsGenTrkHist->Draw("COLZ"); + vtxResYvsGenTrkHist->SetTitle("Vertex Resolution Y vs MC Tracks"); + vtxResYvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}"); + vtxResYvsGenTrkHist->GetYaxis()->SetTitle("recVtx_y - mcVtx_y (in mm)"); + vtxResYvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResYvsGenTrkHist->GetYaxis()->CenterTitle(); + + TH2D *vtxResZvsGenTrkHist = new TH2D("vtxResZvsGenTrk","",31,-0.5,30.5,200,-1,1); + vtxResZvsGenTrkHist->Draw("COLZ"); + vtxResZvsGenTrkHist->SetTitle("Vertex Resolution Z vs MC Tracks"); + vtxResZvsGenTrkHist->GetXaxis()->SetTitle("N_{MC}"); + vtxResZvsGenTrkHist->GetYaxis()->SetTitle("recVtx_z - mcVtx_z (in mm)"); + vtxResZvsGenTrkHist->GetXaxis()->CenterTitle(); vtxResZvsGenTrkHist->GetYaxis()->CenterTitle(); + + TH2D *vtxResXvsRecoTrkHist = new TH2D("vtxResXvsRecoTrk","",31,-0.5,30.5,200,-1,1); + vtxResXvsRecoTrkHist->Draw("COLZ"); + vtxResXvsRecoTrkHist->SetTitle("Vertex Resolution X vs Reconstructed Tracks"); + vtxResXvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}"); + vtxResXvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_x - mcVtx_x (in mm)"); + vtxResXvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResXvsRecoTrkHist->GetYaxis()->CenterTitle(); + + TH2D *vtxResYvsRecoTrkHist = new TH2D("vtxResYvsRecoTrk","",31,-0.5,30.5,200,-1,1); + vtxResYvsRecoTrkHist->Draw("COLZ"); + vtxResYvsRecoTrkHist->SetTitle("Vertex Resolution Y vs Reconstructed Tracks"); + vtxResYvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}"); + vtxResYvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_y - mcVtx_y (in mm)"); + vtxResYvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResYvsRecoTrkHist->GetYaxis()->CenterTitle(); + + TH2D *vtxResZvsRecoTrkHist = new TH2D("vtxResZvsRecoTrk","",31,-0.5,30.5,200,-1,1); + vtxResZvsRecoTrkHist->Draw("COLZ"); + vtxResZvsRecoTrkHist->SetTitle("Vertex Resolution Z vs Reconstructed Tracks"); + vtxResZvsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}"); + vtxResZvsRecoTrkHist->GetYaxis()->SetTitle("recVtx_z - mcVtx_z (in mm)"); + vtxResZvsRecoTrkHist->GetXaxis()->CenterTitle(); vtxResZvsRecoTrkHist->GetYaxis()->CenterTitle(); + + TH1D *numGenTrkswithVtxHist = new TH1D("numGenTrkswithVtx","",31,-0.5,30.5); + TH1D *numRecoTrkswithVtxHist = new TH1D("numRecoTrkswithVtx","",31,-0.5,30.5); + + + int counter(0); + + //Loop over events + std::cout<<"Analyzing "<GetEntries()<<" events!"<Fill(mcEvtVtx.x(), mcEvtVtx.y()); + TVector3 mcRadius(mcEvtVtx.x(), mcEvtVtx.y(), 0); + genVtxRvsZHist->Fill(mcEvtVtx.z(), mcRadius.Mag()); + + //Filtering MC Tracks + int numMCTracks=0; + for(unsigned int i=0; i 1e-4) continue; + + TVector3 mcPartMom(mcMomX[i], mcMomY[i], mcMomZ[i]); + if(fabs(mcPartMom.Eta()) > 3.5) continue; + + numMCTracks++; + } + numGenTracksHist->Fill(numMCTracks); + + ////////////////////////////////////////////////////////////////////////// + ////////////////////// Analyze Reconstructed Tracks ////////////////////// + ////////////////////////////////////////////////////////////////////////// + + numRecoTracksHist->Fill(recoType.GetSize()); + + //Finding Reconstructed Vertex and Vertexing Efficiency + int nVtx=0; + float diff=999.; + int nAssoPart=0; + TVector3 recoEvtVtx(-999., -999., -999.); + for(unsigned int i=0; iFill(nVtx); + recoVtxYvsXHist->Fill(recoEvtVtx.x(), recoEvtVtx.y()); + TVector3 recoRadius(recoEvtVtx.x(), recoEvtVtx.y(), 0); + recoVtxRvsZHist->Fill(recoEvtVtx.z(), recoRadius.Mag()); + + vtxResXvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.x() - mcEvtVtx.x()); + vtxResYvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.y() - mcEvtVtx.y()); + vtxResZvsGenTrkHist->Fill(numMCTracks, recoEvtVtx.z() - mcEvtVtx.z()); + + vtxResXvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.x() - mcEvtVtx.x()); + vtxResYvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.y() - mcEvtVtx.y()); + vtxResZvsRecoTrkHist->Fill(nAssoPart, recoEvtVtx.z() - mcEvtVtx.z()); + + if(nVtx !=0) { + numGenTrkswithVtxHist->Fill(numMCTracks); + numRecoTrkswithVtxHist->Fill(recoType.GetSize()); + + recoVsMCTracksHist->Fill(numMCTracks, nAssoPart);} + + } + + //-------------------------------------------------------------------------------------------------------------------------------------------- + ofile->Write(); // Write histograms to file + ofile->Close(); // Close output file + +} + diff --git a/benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx b/benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx new file mode 100644 index 00000000..ba206ec9 --- /dev/null +++ b/benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx @@ -0,0 +1,379 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fmt/color.h" +#include "fmt/core.h" + +#include "nlohmann/json.hpp" + +double StudentT(double *x, double *par); + +void vtx_dis_plots(const std::string& config_name) +{ + // Read our configuration + std::ifstream config_file{config_name}; + nlohmann::json config; + config_file >> config; + + const std::string hists_file = config["hists_file"]; + const std::string detector = config["detector"]; + const std::string output_prefix = config["output_prefix"]; + const int ebeam = config["ebeam"]; + const int pbeam = config["pbeam"]; + const int Q2_min = config["Min_Q2"]; + const int nfiles = config["nfiles"]; + + fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green), + "Plotting DIS tracking analysis...\n"); + fmt::print(" - Detector package: {}\n", detector); + fmt::print(" - input file for histograms: {}\n", hists_file); + fmt::print(" - output prefix for plots: {}\n", output_prefix); + fmt::print(" - ebeam: {}\n", ebeam); + fmt::print(" - pbeam: {}\n", pbeam); + fmt::print(" - Minimum Q2: {}\n", Q2_min); + fmt::print(" - nfiles: {}\n", nfiles); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Read file with histograms + TFile* file = new TFile(hists_file.c_str()); + + std::cout<<"Reading histograms..."<Get("recoVsMCTracks"); + TH2D* heff = (TH2D*) file->Get("recoVtxEff"); + + TH2D* hg1 = (TH2D*) file->Get("genVtxYvsXHist"); + TH2D* hg2 = (TH2D*) file->Get("genVtxRvsZHist"); + + TH2D* hr1 = (TH2D*) file->Get("recoVtxYvsX"); + TH2D* hr2 = (TH2D*) file->Get("recoVtxRvsZ"); + + TH2D* hres1r = (TH2D*) file->Get("vtxResXvsGenTrk"); + TH2D* hres2r = (TH2D*) file->Get("vtxResYvsGenTrk"); + TH2D* hres3r = (TH2D*) file->Get("vtxResZvsGenTrk"); + + TH2D* hres1g = (TH2D*) file->Get("vtxResXvsRecoTrk"); + TH2D* hres2g = (TH2D*) file->Get("vtxResYvsRecoTrk"); + TH2D* hres3g = (TH2D*) file->Get("vtxResZvsRecoTrk"); + + TH1D* hng1 = (TH1D*) file->Get("numGenTracks"); + TH1D* hng2 = (TH1D*) file->Get("numGenTrkswithVtx"); + + TH1D* hnr1 = (TH1D*) file->Get("numRecoTracks"); + TH1D* hnr2 = (TH1D*) file->Get("numRecoTrkswithVtx"); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + //Fit Function + TF1 *myfunction = new TF1("fit", StudentT, -2, 2, 4); + + //Fitting res v/s MC histograms + myfunction->SetParameters(hres1g->GetMaximum(), 0, 0.05, 1); + hres1g->FitSlicesY(myfunction, 0, -1, 10); + + myfunction->SetParameters(hres2g->GetMaximum(), 0, 0.05, 1); + hres2g->FitSlicesY(myfunction, 0, -1, 10); + + myfunction->SetParameters(hres3g->GetMaximum(), 0, 0.05, 1); + hres3g->FitSlicesY(myfunction, 0, -1, 10); + + //Fitting res v/s RC histograms + myfunction->SetParameters(hres1r->GetMaximum(), 0, 0.05, 1); + hres1r->FitSlicesY(myfunction, 0, -1, 10); + + myfunction->SetParameters(hres2r->GetMaximum(), 0, 0.05, 1); + hres2r->FitSlicesY(myfunction, 0, -1, 10); + + myfunction->SetParameters(hres3r->GetMaximum(), 0, 0.05, 1); + hres3r->FitSlicesY(myfunction, 0, -1, 10); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Make new efficiency histograms + TH1D *vtxEffVsGenTrkHist = new TH1D("vtxEffVsGenTrk","",31,-0.5,30.5); + TH1D *vtxEffVsRecoTrkHist = new TH1D("vtxEffVsRecoTrk","",31,-0.5,30.5); + + // Fill efficiency v/s MC histogram + for(int i=0; i<=hng1->GetNbinsX(); i++) + { + float neventsMC = hng1->GetBinContent(i); + float nvtxevtsMC = hng2->GetBinContent(i); + + if(neventsMC != 0) + { + vtxEffVsGenTrkHist->SetBinContent(i, nvtxevtsMC/neventsMC); + vtxEffVsGenTrkHist->SetBinError(i, sqrt((nvtxevtsMC+1)/(neventsMC+2)*((nvtxevtsMC+2)/(neventsMC+3)-(nvtxevtsMC+1)/(neventsMC+2)))); + } + } + + vtxEffVsGenTrkHist->SetMarkerColor(kRed); + vtxEffVsGenTrkHist->SetMarkerStyle(8); vtxEffVsGenTrkHist->SetMarkerSize(1.2); + vtxEffVsGenTrkHist->SetTitle("Vertexing Efficiency vs MC Tracks"); + vtxEffVsGenTrkHist->GetXaxis()->SetTitle("N_{MC}"); + vtxEffVsGenTrkHist->GetYaxis()->SetTitle("Vertexing Efficiency"); + vtxEffVsGenTrkHist->GetYaxis()->SetRangeUser(0, 1.2); + + //Fill efficiency v/s RC histogram + for(int i=0; i<=hnr1->GetNbinsX(); i++) + { + float neventsRC = hnr1->GetBinContent(i); + float nvtxevtsRC = hnr2->GetBinContent(i); + + if(neventsRC != 0) + { + vtxEffVsRecoTrkHist->SetBinContent(i, nvtxevtsRC/neventsRC); + vtxEffVsRecoTrkHist->SetBinError(i, sqrt((nvtxevtsRC+1)/(neventsRC+2)*((nvtxevtsRC+2)/(neventsRC+3)-(nvtxevtsRC+1)/(neventsRC+2)))); + } + } + + vtxEffVsRecoTrkHist->SetMarkerColor(kRed); + vtxEffVsRecoTrkHist->SetMarkerStyle(8); vtxEffVsRecoTrkHist->SetMarkerSize(1.2); + vtxEffVsRecoTrkHist->SetTitle("Vertexing Efficiency vs RC Tracks"); + vtxEffVsRecoTrkHist->GetXaxis()->SetTitle("N_{RC}"); + vtxEffVsRecoTrkHist->GetYaxis()->SetTitle("Vertexing Efficiency"); + vtxEffVsRecoTrkHist->GetYaxis()->SetRangeUser(0, 1.2); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Make plots and save to PDF file + + // Update Style + gStyle->SetPadGridX(0); + gStyle->SetPadGridY(0); + gStyle->SetOptStat(0); + + std::cout<<"Making plots..."<cd(1); + auto func = new TF1("func","x",-1,40); + func->SetLineStyle(2); func->SetLineColor(1); + + hgVr->Draw(); + func->Draw("SAMEL"); + + //Vertexing Efficiency + TCanvas *c2 = new TCanvas("c2","Vertexing Efficiency",800,600); + int nevents(0); + for(int i=0; i<=heff->GetNbinsX(); i++) nevents = nevents + heff->GetBinContent(i); + heff->Scale(100./nevents); + + c2->cd(1); + heff->Draw("pe"); + + // MC vx versus vy Vertex + TCanvas *c3 = new TCanvas("c3","MC Vertices",800,600); + c3->cd(1); + hg1->Draw(); + + // MC vr versus vz Vertex + TCanvas *c4 = new TCanvas("c4","MC Vertices",800,600); + c4->cd(1); + hg2->Draw(); + + // Reconstructed Vertex vx versus vy + TCanvas *c5 = new TCanvas("c5","Reconstructed Vertices",800,600); + c5->cd(1); + hr1->Draw(); + + // Reconstructed Vertex vr versus vz + TCanvas *c6 = new TCanvas("c6","Reconstructed Vertices",800,600); + c6->cd(1); + hr2->Draw(); + + //Vertex Resolution vs MC Tracks + TCanvas *c7 = new TCanvas("c7","VtxResX vs MCTrks",800,600); + c7->cd(1); + hres1g->Draw("COLZ"); + + TCanvas *c8 = new TCanvas("c8","VtxResY vs MCTrks",800,600); + c8->cd(1); + hres2g->Draw("COLZ"); + + TCanvas *c9 = new TCanvas("c9","VtxResZ vs MCTrks",800,600); + c9->cd(1); + hres3g->Draw("COLZ"); + + //Vertex Resolution vs RC Tracks + TCanvas *c10 = new TCanvas("c10","VtxResX vs RCTrks",800,600); + c10->cd(1); + hres1r->Draw("COLZ"); + + TCanvas *c11 = new TCanvas("c11","VtxResY vs RCTrks",800,600); + c11->cd(1); + hres2r->Draw("COLZ"); + + TCanvas *c12 = new TCanvas("c12","VtxResZ vs RCTrks",800,600); + c12->cd(1); + hres3r->Draw("COLZ"); + + // Res Sigma v/s MC tracks + TCanvas *c13 = new TCanvas("c13","Vertex Resolution vs MC Tracks",800,600); + c13->cd(1); + + TH1D *resXsigmaM = (TH1D*)gDirectory->Get("vtxResXvsGenTrk_2"); + TH1D *resYsigmaM = (TH1D*)gDirectory->Get("vtxResYvsGenTrk_2"); + TH1D *resZsigmaM = (TH1D*)gDirectory->Get("vtxResZvsGenTrk_2"); + + resXsigmaM->SetMarkerStyle(20); resYsigmaM->SetMarkerStyle(21); resZsigmaM->SetMarkerStyle(22); + resXsigmaM->SetMarkerSize(1.2); resYsigmaM->SetMarkerSize(1.2); resZsigmaM->SetMarkerSize(1.2); + resXsigmaM->SetMarkerColor(kBlue); resYsigmaM->SetMarkerColor(kRed); resZsigmaM->SetMarkerColor(kBlack); + resXsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks"); resYsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks"); + resZsigmaM->SetTitle("Vertex Resolution Sigma vs MC Tracks"); + resZsigmaM->GetXaxis()->SetTitle("N_{MC}"); + resXsigmaM->GetYaxis()->SetTitle("#sigma (mm)"); resYsigmaM->GetYaxis()->SetTitle("#sigma (mm)"); resZsigmaM->GetYaxis()->SetTitle("#sigma (mm)"); + resXsigmaM->GetYaxis()->SetRangeUser(0, 1); resYsigmaM->GetYaxis()->SetRangeUser(0, 1); resZsigmaM->GetYaxis()->SetRangeUser(0, 1); + + resXsigmaM->Draw("P"); + resYsigmaM->Draw("PSAME"); + resZsigmaM->Draw("PSAME"); + + TLegend *legend1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed + legend1->AddEntry(resXsigmaM, "v_{x}", "lep"); + legend1->AddEntry(resYsigmaM, "v_{y}", "lep"); + legend1->AddEntry(resZsigmaM, "v_{z}", "lep"); + legend1->Draw(); + + // Res Mean v/s MC tracks + TCanvas *c14 = new TCanvas("c14","Vertex Resolution Mean vs MC Tracks",800,600); + c14->cd(1); + + TH1D *resXmeanM = (TH1D*)gDirectory->Get("vtxResXvsGenTrk_1"); + TH1D *resYmeanM = (TH1D*)gDirectory->Get("vtxResYvsGenTrk_1"); + TH1D *resZmeanM = (TH1D*)gDirectory->Get("vtxResZvsGenTrk_1"); + + resXmeanM->SetMarkerStyle(20); resYmeanM->SetMarkerStyle(21); resZmeanM->SetMarkerStyle(22); + resXmeanM->SetMarkerSize(1.2); resYmeanM->SetMarkerSize(1.2); resZmeanM->SetMarkerSize(1.2); + resXmeanM->SetMarkerColor(kBlue); resYmeanM->SetMarkerColor(kRed); resZmeanM->SetMarkerColor(kBlack); + resXmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks"); resYmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks"); + resZmeanM->SetTitle("Vertex Resolution Mean vs MC Tracks"); + resZmeanM->GetXaxis()->SetTitle("N_{MC}"); + resXmeanM->GetYaxis()->SetTitle("#mu (mm)"); resYmeanM->GetYaxis()->SetTitle("#mu (mm)"); resZmeanM->GetYaxis()->SetTitle("#mu (mm)"); + resXmeanM->GetYaxis()->SetRangeUser(-1, 1); resYmeanM->GetYaxis()->SetRangeUser(-1, 1); resZmeanM->GetYaxis()->SetRangeUser(-1, 1); + + resXmeanM->Draw("P"); + resYmeanM->Draw("PSAME"); + resZmeanM->Draw("PSAME"); + + TLegend *legend2 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed + legend2->AddEntry(resXmeanM, "v_{x}", "lep"); + legend2->AddEntry(resYmeanM, "v_{y}", "lep"); + legend2->AddEntry(resZmeanM, "v_{z}", "lep"); + legend2->Draw(); + + //Res Sigma v/s RC tracks + TCanvas *c15 = new TCanvas("c15","Vertex Resolution vs RC Tracks",800,600); + c15->cd(1); + + TH1D *resXsigmaR = (TH1D*)gDirectory->Get("vtxResXvsRecoTrk_2"); + TH1D *resYsigmaR = (TH1D*)gDirectory->Get("vtxResYvsRecoTrk_2"); + TH1D *resZsigmaR = (TH1D*)gDirectory->Get("vtxResZvsRecoTrk_2"); + + resXsigmaR->SetMarkerStyle(20); resYsigmaR->SetMarkerStyle(21); resZsigmaR->SetMarkerStyle(22); + resXsigmaR->SetMarkerSize(1.2); resYsigmaR->SetMarkerSize(1.2); resZsigmaR->SetMarkerSize(1.2); + resXsigmaR->SetMarkerColor(kBlue); resYsigmaR->SetMarkerColor(kRed); resZsigmaR->SetMarkerColor(kBlack); + resXsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks"); resYsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks"); + resZsigmaR->SetTitle("Vertex Resolution Sigma vs RC Tracks"); + resZsigmaR->GetXaxis()->SetTitle("N_{RC}"); + resXsigmaR->GetYaxis()->SetTitle("#sigma (mm)"); resYsigmaR->GetYaxis()->SetTitle("#sigma (mm)"); resZsigmaR->GetYaxis()->SetTitle("#sigma (mm)"); + resXsigmaR->GetYaxis()->SetRangeUser(0, 1); resYsigmaR->GetYaxis()->SetRangeUser(0, 1); resZsigmaR->GetYaxis()->SetRangeUser(0, 1); + + resXsigmaR->Draw("P"); + resYsigmaR->Draw("PSAME"); + resZsigmaR->Draw("PSAME"); + + TLegend *legend3 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed + legend3->AddEntry(resXsigmaR, "v_{x}", "lep"); + legend3->AddEntry(resYsigmaR, "v_{y}", "lep"); + legend3->AddEntry(resZsigmaR, "v_{z}", "lep"); + legend3->Draw(); + + //Res Mean v/s RC tracks + TCanvas *c16 = new TCanvas("c16","Vertex Resolution Mean vs RC Tracks",800,600); + c16->cd(1); + + TH1D *resXmeanR = (TH1D*)gDirectory->Get("vtxResXvsRecoTrk_1"); + TH1D *resYmeanR = (TH1D*)gDirectory->Get("vtxResYvsRecoTrk_1"); + TH1D *resZmeanR = (TH1D*)gDirectory->Get("vtxResZvsRecoTrk_1"); + + resXmeanR->SetMarkerStyle(20); resYmeanR->SetMarkerStyle(21); resZmeanR->SetMarkerStyle(22); + resXmeanR->SetMarkerSize(1.2); resYmeanR->SetMarkerSize(1.2); resZmeanR->SetMarkerSize(1.2); + resXmeanR->SetMarkerColor(kBlue); resYmeanR->SetMarkerColor(kRed); resZmeanR->SetMarkerColor(kBlack); + resXmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks"); resYmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks"); + resZmeanR->SetTitle("Vertex Resolution Mean vs RC Tracks"); + resZmeanR->GetXaxis()->SetTitle("N_{RC}"); + resXmeanR->GetYaxis()->SetTitle("#mu (mm)"); resYmeanR->GetYaxis()->SetTitle("#mu (mm)"); resZmeanR->GetYaxis()->SetTitle("#mu (mm)"); + resXmeanR->GetYaxis()->SetRangeUser(-1, 1); resYmeanR->GetYaxis()->SetRangeUser(-1, 1); resZmeanR->GetYaxis()->SetRangeUser(-1, 1); + + resXmeanR->Draw("P"); + resYmeanR->Draw("PSAME"); + resZmeanR->Draw("PSAME"); + + TLegend *legend4 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed + legend4->AddEntry(resXmeanR, "v_{x}", "lep"); + legend4->AddEntry(resYmeanR, "v_{y}", "lep"); + legend4->AddEntry(resZmeanR, "v_{z}", "lep"); + legend4->Draw(); + + //Vertexing Efficiency vs MC Tracks + TCanvas *c17 = new TCanvas("c17","Vertexing Efficiency vs MC Tracks",800,600); + c17->cd(1); + + vtxEffVsGenTrkHist->Draw("P"); + + //Vertexing Efficiency vs RC Tracks + TCanvas *c18 = new TCanvas("c18","Vertexing Efficiency vs RC Tracks",800,600); + c18->cd(1); + + vtxEffVsRecoTrkHist->Draw("P"); + + //-------------------------------------------------------------------------------------------------------------------------------------------- + + // Print plots to pdf file + c1->Print(fmt::format("{}.pdf[", output_prefix).c_str()); + c1->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c2->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c3->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c4->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c5->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c6->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c7->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c8->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c9->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c10->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c11->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c12->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c13->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c14->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c15->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c16->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c17->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c18->Print(fmt::format("{}.pdf", output_prefix).c_str()); + c18->Print(fmt::format("{}.pdf]", output_prefix).c_str()); + +} + +//Defining Fitting Function +double StudentT(double *x, double *par){ + double norm = par[0]; + double mean = par[1]; + double sigma = par[2]; + double nu = par[3]; + + double pi = 3.14; + 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) ); + return st; +}