From bd1cd146c7266b0e1f91f649207fe4ffcd8edc19 Mon Sep 17 00:00:00 2001
From: khushi-singla-21 <142303997+khushi-singla-21@users.noreply.github.com>
Date: Tue, 12 Nov 2024 00:46:31 +0530
Subject: [PATCH] Added benchmark for vertexing studies in DIS events (#101)
---
.../tracking_performances_dis/README.md | 2 +-
.../tracking_performances_dis/Snakefile | 96 ++++-
.../analysis/vtx_dis_analysis.cxx | 307 ++++++++++++++
.../analysis/vtx_dis_plots.cxx | 379 ++++++++++++++++++
.../tracking_performances_dis/config.yml | 8 +
5 files changed, 774 insertions(+), 18 deletions(-)
create mode 100644 benchmarks/tracking_performances_dis/analysis/vtx_dis_analysis.cxx
create mode 100644 benchmarks/tracking_performances_dis/analysis/vtx_dis_plots.cxx
diff --git a/benchmarks/tracking_performances_dis/README.md b/benchmarks/tracking_performances_dis/README.md
index 7df7078e..3d04f31e 100644
--- a/benchmarks/tracking_performances_dis/README.md
+++ b/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)
diff --git a/benchmarks/tracking_performances_dis/Snakefile b/benchmarks/tracking_performances_dis/Snakefile
index 7774f2a1..2d521730 100644
--- a/benchmarks/tracking_performances_dis/Snakefile
+++ b/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]}"
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 <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
+
+}
+
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 <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;
+}
diff --git a/benchmarks/tracking_performances_dis/config.yml b/benchmarks/tracking_performances_dis/config.yml
index 066e8df7..76eeacec 100644
--- a/benchmarks/tracking_performances_dis/config.yml
+++ b/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
+
--
GitLab