Added benchmark for vertexing studies in DIS events (#101)

benchmarks/tracking_performances_dis/README.md

@@ -7,4 +7,4 @@ This benchmark generates and reconstructs DIS events. The generated charged-part
 ## Contact
 [Barak Schmookler](baraks@ucr.edu)
 
-
+[Khushi Singla](ksingla.phy@gmail.com)

benchmarks/tracking_performances_dis/Snakefile

@@ -16,6 +16,8 @@ rule trk_dis_compile:
     input:
         "benchmarks/tracking_performances_dis/analysis/trk_dis_analysis_cxx.so",
         "benchmarks/tracking_performances_dis/analysis/trk_dis_plots_cxx.so"
+	"benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis_cxx.so",
+        "benchmarks/tracking_performances_dis/analysis/vtx_dis_plots_cxx.so"
 
 # Process the generated HepMC files through the simulation
 rule trk_dis_sim:
@@ -49,46 +51,65 @@ rule trk_dis_reco:
 set -m # monitor mode to prevent lingering processes
 exec env DETECTOR_CONFIG={wildcards.DETECTOR_CONFIG} \
   eicrecon {input} -Ppodio:output_file={output} \
-  -Ppodio:output_collections=MCParticles,ReconstructedChargedParticles,ReconstructedTruthSeededChargedParticles,CentralCKFTrackAssociations,CentralCKFTruthSeededTrackAssociations
+  -Ppodio:output_collections=MCParticles,ReconstructedChargedParticles,ReconstructedTruthSeededChargedParticles,CentralCKFTrackAssociations,CentralCKFTruthSeededTrackAssociations,CentralTrackVertices
 """
 
 # Process the files -- either from the campaign or local running -- through the analysis script
-rule trk_dis_analysis:
+rule dis_analysis:
     input:
-        script="benchmarks/tracking_performances_dis/analysis/trk_dis_analysis.cxx",
-        script_compiled="benchmarks/tracking_performances_dis/analysis/trk_dis_analysis_cxx.so",
+        script_trk="benchmarks/tracking_performances_dis/analysis/trk_dis_analysis.cxx",
+        script_trk_compiled="benchmarks/tracking_performances_dis/analysis/trk_dis_analysis_cxx.so",
+        script_vtx="benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx",
+        script_vtx_compiled="benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis_cxx.so",
         data="sim_output/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_beamEffects_xAngle=-0.025_hiDiv_{INDEX}.edm4eic.root",
     output:
-        config="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/config.json",
-        hists="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/hists.root",
+        config_trk="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/config.json",
+        hists_trk="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/hists.root",
+        config_vtx="results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/config.json",
+        hists_vtx="results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_{INDEX}/hists.root",
     wildcard_constraints:
         PREFIX= ".*",
         EBEAM="\d+",
         PBEAM="\d+",
-	MINQ2="\d+",
+        MINQ2="\d+",
         INDEX="\d+",
     shell:
         """
-cat > {output.config} <<EOF
+cat > {output.config_trk} <<EOF
+{{
+  "rec_file": "{input.data}",
+  "detector": "{wildcards.DETECTOR_CONFIG}",
+  "ebeam": {wildcards.EBEAM},
+  "pbeam": {wildcards.PBEAM},
+  "Min_Q2": {wildcards.MINQ2},
+  "output_prefix": "$(dirname "{output.hists_trk}")/hists"
+}}
+EOF
+root -l -b -q '{input.script_trk}+("{output.config_trk}")'
+
+cat > {output.config_vtx} <<EOF
 {{
   "rec_file": "{input.data}",
   "detector": "{wildcards.DETECTOR_CONFIG}",
   "ebeam": {wildcards.EBEAM},
   "pbeam": {wildcards.PBEAM},
   "Min_Q2": {wildcards.MINQ2},
-  "output_prefix": "$(dirname "{output.hists}")/hists"
+  "output_prefix": "$(dirname "{output.hists_vtx}")/hists"
 }}
 EOF
-root -l -b -q '{input.script}+("{output.config}")'
+root -l -b -q '{input.script_vtx}+("{output.config_vtx}")'
 """
 
 #Merge all the files produced in the previous step
-rule trk_dis_combine:
+rule dis_combine:
     input:
         lambda wildcards: [f"results/tracking_performances_dis/{wildcards.DETECTOR_CONFIG}/{wildcards.PREFIX}pythia8NCDIS_{wildcards.EBEAM}x{wildcards.PBEAM}_minQ2={wildcards.MINQ2}_{ix}/hists.root" for ix in range(1,int(wildcards.NUM_FILES)+1)],
+        lambda wildcards: [f"results/vertexing_performances_dis/{wildcards.DETECTOR_CONFIG}/{wildcards.PREFIX}pythia8NCDIS_{wildcards.EBEAM}x{wildcards.PBEAM}_minQ2={wildcards.MINQ2}_{ix}/hists.root" for ix in range(1,int(wildcards.NUM_FILES)+1)],
     output:
-        config="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/config.json",
-        hists="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/hists.root",
+        config_trk="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/config.json",
+        hists_trk="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/hists.root",
+        config_vtx="results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/config.json",
+        hists_vtx="results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/hists.root",
     wildcard_constraints:
         PREFIX= ".*",
         EBEAM="\d+",
@@ -97,18 +118,31 @@ rule trk_dis_combine:
         NUM_FILES="\d+",
     shell:
         """
-cat > {output.config} <<EOF
+cat > {output.config_trk} <<EOF
+{{
+  "hists_file": "{output.hists_trk}",
+  "detector": "{wildcards.DETECTOR_CONFIG}",
+  "ebeam": {wildcards.EBEAM},
+  "pbeam": {wildcards.PBEAM},
+  "Min_Q2": {wildcards.MINQ2},
+  "nfiles": {wildcards.NUM_FILES},
+  "output_prefix": "$(dirname "{output.hists_trk}")/plots"
+}}
+EOF
+hadd {output.hists_trk} {input}
+
+cat > {output.config_vtx} <<EOF
 {{
-  "hists_file": "{output.hists}",
+  "hists_file": "{output.hists_vtx}",
   "detector": "{wildcards.DETECTOR_CONFIG}",
   "ebeam": {wildcards.EBEAM},
   "pbeam": {wildcards.PBEAM},
   "Min_Q2": {wildcards.MINQ2},
   "nfiles": {wildcards.NUM_FILES},
-  "output_prefix": "$(dirname "{output.hists}")/plots"
+  "output_prefix": "$(dirname "{output.hists_vtx}")/plots"
 }}
 EOF
-hadd {output.hists} {input}
+hadd {output.hists_vtx} {input}
 """
 
 #Process the merged file through the plotting script
@@ -130,6 +164,23 @@ rule trk_dis_plots:
 root -l -b -q '{input.script}+("{input.config}")'
 """
 
+rule vtx_dis_plots:
+    input:
+        script="benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx",
+        script_compiled="benchmarks/tracking_performances_dis/analysis/vtx_dis_plots_cxx.so",
+        config="results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/config.json",
+    output:
+        "results/vertexing_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/plots.pdf"
+    wildcard_constraints:
+        PREFIX= ".*",
+        EBEAM="\d+",
+        PBEAM="\d+",
+	MINQ2="\d+",
+        NUM_FILES="\d+",
+    shell:
+        """
+root -l -b -q '{input.script}+("{input.config}")'
+"""
 
 #Examples of invocation
 rule trk_dis_run_locally:
@@ -138,6 +189,12 @@ rule trk_dis_run_locally:
     message:
         "See output in {input[0]}"
 
+rule vtx_dis_run_locally:
+    input:
+        "results/vertexing_performances_dis/" + os.environ["DETECTOR_CONFIG"] + "/pythia8NCDIS_18x275_minQ2=1_combined_5/plots.pdf",
+    message:
+        "See output in {input[0]}"
+
 
 rule trk_dis_run_locally_trk_only:
     input:
@@ -146,3 +203,8 @@ rule trk_dis_run_locally_trk_only:
         "See output in {input[0]}"
 
 
+rule vtx_dis_run_locally_trk_only:
+    input:
+        "results/vertexing_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/plots.pdf",
+    message:
+        "See output in {input[0]}"

benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx

+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <TChain.h>
+#include <TFile.h>
+#include <TTreeReader.h>
+#include <TTreeReaderArray.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <Math/Vector3D.h>
+#include <Math/Vector4D.h>
+#include <Math/RotationY.h>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TLorentzVector.h>
+
+#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<int> recoVtxType = {tree_reader, "CentralTrackVertices.type"};
+    TTreeReaderArray<float> recoVtxX = {tree_reader, "CentralTrackVertices.position.x"};
+    TTreeReaderArray<float> recoVtxY = {tree_reader, "CentralTrackVertices.position.y"};
+    TTreeReaderArray<float> recoVtxZ = {tree_reader, "CentralTrackVertices.position.z"};
+  
+    TTreeReaderArray<unsigned int> assoPartBegin = {tree_reader, "CentralTrackVertices.associatedParticles_begin"};
+    TTreeReaderArray<unsigned int> assoPartEnd = {tree_reader, "CentralTrackVertices.associatedParticles_end"};
+    TTreeReaderArray<int> assoPartIndex = {tree_reader, "_CentralTrackVertices_associatedParticles.index"};
+
+    // MC
+    TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
+    TTreeReaderArray<int> mcPDG = {tree_reader, "MCParticles.PDG"};
+    TTreeReaderArray<float> mcCharge = {tree_reader, "MCParticles.charge"};
+    TTreeReaderArray<float> mcMomX = {tree_reader, "MCParticles.momentum.x"};
+    TTreeReaderArray<float> mcMomY = {tree_reader, "MCParticles.momentum.y"};
+    TTreeReaderArray<float> mcMomZ = {tree_reader, "MCParticles.momentum.z"};
+
+    TTreeReaderArray<double> mcVtxX = {tree_reader, "MCParticles.vertex.x"};
+    TTreeReaderArray<double> mcVtxY = {tree_reader, "MCParticles.vertex.y"};
+    TTreeReaderArray<double> mcVtxZ = {tree_reader, "MCParticles.vertex.z"};
+
+    TTreeReaderArray<unsigned int> mcParentBegin = {tree_reader, "MCParticles.parents_begin"};
+    TTreeReaderArray<unsigned int> mcParentEnd = {tree_reader, "MCParticles.parents_end"};
+    TTreeReaderArray<int> mcParentIndex = {tree_reader, "_MCParticles_parents.index"};
+
+    // Reco
+    TTreeReaderArray<int> 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 "<<mychain->GetEntries()<<" events!"<<std::endl;
+    while (tree_reader.Next()) {
+
+	    if(counter%100==0) std::cout<<"Analyzing event "<<counter<<std::endl;
+	    
+	    counter++;
+
+    	//////////////////////////////////////////////////////////////////////////
+    	////////////////////////  Analyze MC Tracks  /////////////////////////
+    	//////////////////////////////////////////////////////////////////////////
+    
+    	//Finding MC vertex using scattered electron
+    	TVector3 mcEvtVtx(-999., -999., -999.);    
+    	for(unsigned int i=0; i<mcGenStat.GetSize(); i++)
+    	{
+		if(mcGenStat[i] != 1) continue;
+		if(mcPDG[i] != 11) continue;
+
+		bool scatEfound = false;
+		// mcParentBegin and mcParentEnd specify the entries from _MCParticles_parents.index 
+		// _MCParticles_parents.index stores the MCParticle index
+		for(unsigned int j=mcParentBegin[i]; j<mcParentEnd[i]; j++)
+		{
+          	int parentPDG = mcPDG[mcParentIndex[j]];
+          	if(parentPDG == 11) scatEfound = true;
+        	}
+       
+		if(scatEfound == false) continue;
+		//Scattered electron found
+		double vtx_mc_x =  mcVtxX[i];
+		double vtx_mc_y =  mcVtxY[i];
+		double vtx_mc_z =  mcVtxZ[i];
+		mcEvtVtx = TVector3(vtx_mc_x, vtx_mc_y, vtx_mc_z);
+    	}
+    	genVtxYvsXHist->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<mcGenStat.GetSize(); i++)
+    	{
+		if(mcGenStat[i] != 1) continue;
+		if(mcCharge[i] == 0) continue;
+	
+		TVector3 mcPartVtx(mcVtxX[i], mcVtxY[i], mcVtxZ[i]);
+		TVector3 vtx_diff = mcPartVtx - mcEvtVtx;
+		if(vtx_diff.Mag() > 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; i<recoVtxType.GetSize(); i++)
+    	{
+		nVtx++;
+	
+		TVector3 recoVtx(recoVtxX[i], recoVtxY[i], recoVtxZ[i]);
+	
+		//Finding the reconstructed vertex closer to the MC vertex
+		TVector3 vtx_diff = recoVtx - mcEvtVtx;
+		if(vtx_diff.Mag() < diff)
+		{
+	    	diff = vtx_diff.Mag();
+	    	recoEvtVtx = recoVtx;
+	    
+	    	for(unsigned int j=assoPartBegin[i]; j<assoPartEnd[i]; j++)
+	    	{
+              	nAssoPart = j;
+            	}
+		}
+    	}
+    
+    	recoVtxEffHist->Fill(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
+
+}
+  

benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx

+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TStyle.h>
+#include <TLegend.h>
+
+#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..."<<std::endl;
+
+    TH2D* hgVr = (TH2D*) file->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..."<<std::endl;
+  
+    // MC vs RC associated particles
+    TCanvas *c1 = new TCanvas("c1","MC vs Reco Tracks",800,600);
+    c1->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;
+}

benchmarks/tracking_performances_dis/config.yml

@@ -15,8 +15,10 @@ bench:tracking_performances_dis:
     - ["sim:tracking_performances_dis"]
   script:
     - snakemake $SNAKEMAKE_FLAGS --cores 1 trk_dis_run_locally_trk_only
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 vtx_dis_run_locally_trk_only
     # avoid uploading intermediate results
     - find results/tracking_performances_dis/ -mindepth 2 -maxdepth 2 -type d ! -name "*combined*" | xargs rm -rfv
+    - find results/vertexing_performances_dis/ -mindepth 2 -maxdepth 2 -type d ! -name "*combined*" | xargs rm -rfv
 
 collect_results:tracking_performances_dis:
   extends: .det_benchmark
@@ -27,9 +29,15 @@ collect_results:tracking_performances_dis:
     - ls -lrht
     - mv results{,_save}/ # move results directory out of the way to preserve it
     - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output trk_dis_run_locally_trk_only
+    - snakemake $SNAKEMAKE_FLAGS --cores 1 --delete-all-output vtx_dis_run_locally_trk_only
     - mv results{_save,}/
     # convert to png
     - |
       gs -sDEVICE=pngalpha -dUseCropBox -r144 \
         -o 'results/tracking_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/plot_%03d.png' \
         results/tracking_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/plots.pdf
+    - |
+      gs -sDEVICE=pngalpha -dUseCropBox -r144 \
+        -o 'results/vertexing_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/plot_%03d.png' \
+        results/vertexing_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/plots.pdf
+