diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 86cce78e96ffa83682ea8984fdd000ad0642b8b5..002bb608794dcb43a4a08fd6e8a9d1263aeff90a 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -136,6 +136,7 @@ include:
- local: 'benchmarks/ecal_gaps/config.yml'
- local: 'benchmarks/tracking_detectors/config.yml'
- local: 'benchmarks/tracking_performances/config.yml'
+ - local: 'benchmarks/tracking_performances_dis/config.yml'
- local: 'benchmarks/barrel_ecal/config.yml'
- local: 'benchmarks/barrel_hcal/config.yml'
- local: 'benchmarks/zdc/config.yml'
@@ -158,6 +159,7 @@ deploy_results:
- "collect_results:pid"
- "collect_results:tracking_performance"
- "collect_results:tracking_performance_campaigns"
+ - "collect_results:tracking_performances_dis"
- "collect_results:zdc"
- "collect_results:zdc_lyso"
script:
diff --git a/Snakefile b/Snakefile
index d153297c64c2d0bb19bf544a67d875854045b6d5..738bb5e5289c458a0f957fdf02f3a2047800caed 100644
--- a/Snakefile
+++ b/Snakefile
@@ -6,6 +6,7 @@ include: "benchmarks/barrel_ecal/Snakefile"
include: "benchmarks/ecal_gaps/Snakefile"
include: "benchmarks/material_scan/Snakefile"
include: "benchmarks/tracking_performances/Snakefile"
+include: "benchmarks/tracking_performances_dis/Snakefile"
include: "benchmarks/zdc_lyso/Snakefile"
use_s3 = config["remote_provider"].lower() == "s3"
diff --git a/benchmarks/tracking_performances_dis/README.md b/benchmarks/tracking_performances_dis/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..7df7078ea23a7059f0f08223fda5df72f033bba8
--- /dev/null
+++ b/benchmarks/tracking_performances_dis/README.md
@@ -0,0 +1,10 @@
+Detector Benchmark for tracking performance with DIS events
+============================================================
+
+## Overview
+This benchmark generates and reconstructs DIS events. The generated charged-particle spectra is then compared to the spectra of truth-seeded and real-seeded reconstructed tracks. The hit-based track purity is also plotted.
+
+## Contact
+[Barak Schmookler](baraks@ucr.edu)
+
+
diff --git a/benchmarks/tracking_performances_dis/Snakefile b/benchmarks/tracking_performances_dis/Snakefile
new file mode 100644
index 0000000000000000000000000000000000000000..7774f2a112c6331d49d5d96e64fd46b79b3562a6
--- /dev/null
+++ b/benchmarks/tracking_performances_dis/Snakefile
@@ -0,0 +1,148 @@
+import os
+
+rule compile_analysis:
+ input:
+ "{path}/{filename}.cxx",
+ output:
+ "{path}/{filename}_cxx.d",
+ "{path}/{filename}_cxx.so",
+ "{path}/{filename}_cxx_ACLiC_dict_rdict.pcm",
+ shell:
+ """
+root -l -b -q -e '.L {input}+'
+"""
+
+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"
+
+# Process the generated HepMC files through the simulation
+rule trk_dis_sim:
+ input:
+ warmup="warmup/{DETECTOR_CONFIG}.edm4hep.root",
+ output:
+ "sim_output/tracking_performances_dis/{DETECTOR_CONFIG}/pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_beamEffects_xAngle=-0.025_hiDiv_{INDEX}.edm4hep.root",
+ params:
+ N_EVENTS=200,
+ shell:
+ """
+ddsim \
+ --runType batch \
+ --part.minimalKineticEnergy 1000*GeV \
+ --filter.tracker edep0 \
+ -v WARNING \
+ --numberOfEvents {params.N_EVENTS} \
+ --compactFile $DETECTOR_PATH/{wildcards.DETECTOR_CONFIG}.xml \
+ --inputFiles root://dtn-eic.jlab.org//work/eic2/EPIC/EVGEN/DIS/NC/{wildcards.EBEAM}x{wildcards.PBEAM}/minQ2={wildcards.MINQ2}/pythia8NCDIS_{wildcards.EBEAM}x{wildcards.PBEAM}_minQ2={wildcards.MINQ2}_beamEffects_xAngle=-0.025_hiDiv_vtxfix_{wildcards.INDEX}.hepmc3.tree.root \
+ --outputFile {output}
+"""
+
+# Process the files produced in the previous step through EICRecon
+rule trk_dis_reco:
+ input:
+ "sim_output/tracking_performances_dis/{DETECTOR_CONFIG}/pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_beamEffects_xAngle=-0.025_hiDiv_{INDEX}.edm4hep.root",
+ output:
+ "sim_output/tracking_performances_dis/{DETECTOR_CONFIG}/pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_beamEffects_xAngle=-0.025_hiDiv_{INDEX}.edm4eic.root",
+ shell:
+ """
+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
+"""
+
+# Process the files -- either from the campaign or local running -- through the analysis script
+rule trk_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",
+ 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",
+ wildcard_constraints:
+ PREFIX= ".*",
+ EBEAM="\d+",
+ PBEAM="\d+",
+ MINQ2="\d+",
+ INDEX="\d+",
+ shell:
+ """
+cat > {output.config} <<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"
+}}
+EOF
+root -l -b -q '{input.script}+("{output.config}")'
+"""
+
+#Merge all the files produced in the previous step
+rule trk_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)],
+ 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",
+ wildcard_constraints:
+ PREFIX= ".*",
+ EBEAM="\d+",
+ PBEAM="\d+",
+ MINQ2="\d+",
+ NUM_FILES="\d+",
+ shell:
+ """
+cat > {output.config} <<EOF
+{{
+ "hists_file": "{output.hists}",
+ "detector": "{wildcards.DETECTOR_CONFIG}",
+ "ebeam": {wildcards.EBEAM},
+ "pbeam": {wildcards.PBEAM},
+ "Min_Q2": {wildcards.MINQ2},
+ "nfiles": {wildcards.NUM_FILES},
+ "output_prefix": "$(dirname "{output.hists}")/plots"
+}}
+EOF
+hadd {output.hists} {input}
+"""
+
+#Process the merged file through the plotting script
+rule trk_dis_plots:
+ input:
+ script="benchmarks/tracking_performances_dis/analysis/trk_dis_plots.cxx",
+ script_compiled="benchmarks/tracking_performances_dis/analysis/trk_dis_plots_cxx.so",
+ config="results/tracking_performances_dis/{DETECTOR_CONFIG}/{PREFIX}pythia8NCDIS_{EBEAM}x{PBEAM}_minQ2={MINQ2}_combined_{NUM_FILES}/config.json",
+ output:
+ "results/tracking_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:
+ input:
+ "results/tracking_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:
+ "results/tracking_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/trk_dis_analysis.cxx b/benchmarks/tracking_performances_dis/analysis/trk_dis_analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..52dfdaf70339d40329e3eb08c2e8ec2c3c0fb308
--- /dev/null
+++ b/benchmarks/tracking_performances_dis/analysis/trk_dis_analysis.cxx
@@ -0,0 +1,236 @@
+#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 trk_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());
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Initialize reader
+ TTreeReader tr(mychain);
+
+ // Generated Particle Information
+ TTreeReaderArray<int> gen_status(tr, "MCParticles.generatorStatus");
+ TTreeReaderArray<int> gen_pid(tr, "MCParticles.PDG");
+ TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
+ TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
+ TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
+ TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass"); //Not important here
+ TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
+ TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
+ TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
+ TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
+
+ // Reconstructed real-seeded tracks (charged particles)
+ TTreeReaderArray<float> rec_px(tr, "ReconstructedChargedParticles.momentum.x");
+ TTreeReaderArray<float> rec_py(tr, "ReconstructedChargedParticles.momentum.y");
+ TTreeReaderArray<float> rec_pz(tr, "ReconstructedChargedParticles.momentum.z");
+ TTreeReaderArray<float> rec_mass(tr, "ReconstructedChargedParticles.mass");
+ TTreeReaderArray<int> rec_type(tr, "ReconstructedChargedParticles.type"); //Type 0: successful eta/phi match to generated particle
+ TTreeReaderArray<int> rec_pdg(tr, "ReconstructedChargedParticles.PDG"); //Uses PID lookup table information
+
+ // Reconstructed truth-seeded tracks (charged particles)
+ TTreeReaderArray<float> rec_ts_px(tr, "ReconstructedTruthSeededChargedParticles.momentum.x");
+ TTreeReaderArray<float> rec_ts_py(tr, "ReconstructedTruthSeededChargedParticles.momentum.y");
+ TTreeReaderArray<float> rec_ts_pz(tr, "ReconstructedTruthSeededChargedParticles.momentum.z");
+ TTreeReaderArray<float> rec_ts_mass(tr, "ReconstructedTruthSeededChargedParticles.mass");
+ TTreeReaderArray<int> rec_ts_type(tr, "ReconstructedTruthSeededChargedParticles.type"); //Type 0: successful eta/phi match to generated particle
+ TTreeReaderArray<int> rec_ts_pdg(tr, "ReconstructedTruthSeededChargedParticles.PDG"); //Uses PID lookup table information
+
+ // Hit-based track to MC Particle association weight
+ TTreeReaderArray<float> hit_assoc_weight(tr,"CentralCKFTrackAssociations.weight"); //Real-seeded tracking
+ TTreeReaderArray<float> hit_assoc_weight_ts(tr,"CentralCKFTruthSeededTrackAssociations.weight"); //Truth-seeded tracking
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+ // Define Histograms
+
+ //Eta distribution of generated charged particles
+ TH1 *h1a = new TH1D("h1a","Generated Charged Particles",100,-4,4);
+ h1a->GetXaxis()->SetTitle("#eta_{gen.}");h1a->GetXaxis()->CenterTitle();
+ h1a->SetLineColor(kTeal);h1a->SetLineWidth(2);
+
+ TH1 *h1a1 = new TH1D("h1a1","Generated Charged Particles",100,-4,4); //Minimum momentum cut of Pt > 200 MeV/c
+ h1a1->GetXaxis()->SetTitle("#eta_{gen.}");h1a1->GetXaxis()->CenterTitle();
+ h1a1->SetLineColor(kRed);h1a1->SetLineWidth(2);
+
+ TH1 *h1a2 = new TH1D("h1a2","Generated Charged Particles",100,-4,4); //Minimum momentum cut of Pt > 500 MeV/c
+ h1a2->GetXaxis()->SetTitle("#eta_{gen.}");h1a2->GetXaxis()->CenterTitle();
+ h1a2->SetLineColor(kBlack);h1a2->SetLineWidth(2);
+ h1a2->SetFillColor(kBlack);h1a2->SetFillStyle(3244);
+
+ //Eta distribution of reconstructed real-seeded tracks (charged particles)
+ TH1 *h1b = new TH1D("h1b","Reconstructed Real-seeded tracks",100,-4,4);
+ h1b->GetXaxis()->SetTitle("#eta_{rec.}");h1b->GetXaxis()->CenterTitle();
+ h1b->SetLineColor(kGreen);h1b->SetLineWidth(2);
+
+ TH1 *h1b1 = new TH1D("h1b1","Reconstructed Real-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 200 MeV/c
+ h1b1->GetXaxis()->SetTitle("#eta_{rec.}");h1b1->GetXaxis()->CenterTitle();
+ h1b1->SetLineColor(kBlue);h1b1->SetLineWidth(2);
+
+ TH1 *h1b2 = new TH1D("h1b2","Reconstructed Real-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 500 MeV/c
+ h1b2->GetXaxis()->SetTitle("#eta_{rec.}");h1b2->GetXaxis()->CenterTitle();
+ h1b2->SetLineColor(kRed);h1b2->SetLineWidth(2);
+ h1b2->SetMarkerColor(kRed);h1b2->SetMarkerStyle(kFullCrossX);
+
+ //Eta distribution of reconstructed truth-seeded tracks (charged particles)
+ TH1 *h1c = new TH1D("h1c","Reconstructed Truth-seeded tracks",100,-4,4);
+ h1c->GetXaxis()->SetTitle("#eta_{rec.}");h1c->GetXaxis()->CenterTitle();
+ h1c->SetLineColor(kRed);h1c->SetLineWidth(2);
+
+ TH1 *h1c1 = new TH1D("h1c1","Reconstructed Truth-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 200 MeV/c
+ h1c1->GetXaxis()->SetTitle("#eta_{rec.}");h1c1->GetXaxis()->CenterTitle();
+ h1c1->SetLineColor(kOrange);h1c1->SetLineWidth(2);
+
+ TH1 *h1c2 = new TH1D("h1c2","Reconstructed Truth-seeded tracks",100,-4,4); //Minimum momentum cut of Pt > 500 MeV/c
+ h1c2->GetXaxis()->SetTitle("#eta_{rec.}");h1c2->GetXaxis()->CenterTitle();
+ h1c2->SetLineColor(kMagenta);h1c2->SetLineWidth(2);
+ h1c2->SetMarkerColor(kMagenta);h1c2->SetMarkerStyle(kFullCrossX);
+
+ //Track purity
+ TH1 *h2a = new TH1D("h2a","Real-seeded tracks purity",100,0,1.1);
+ h2a->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h2a->GetXaxis()->CenterTitle();
+ h2a->GetYaxis()->SetTitle("Number of Tracks");h2a->GetYaxis()->CenterTitle();
+ h2a->SetLineColor(kBlack);h2a->SetLineWidth(2);
+
+ TH1 *h2b = new TH1D("h2b","Truth-seeded tracks purity",100,0,1.1);
+ h2b->GetXaxis()->SetTitle("Fraction of track measurements from a given MC Particle");h2b->GetXaxis()->CenterTitle();
+ h2b->GetYaxis()->SetTitle("Number of Tracks");h2b->GetYaxis()->CenterTitle();
+ h2b->SetLineColor(kBlack);h2b->SetLineWidth(2);
+
+ //Define additional variables
+ TLorentzVector gen_vec;
+ TVector3 gen_vertex;
+
+ TLorentzVector rec_vec;
+ TVector3 track_vec; //Reconstructed track momentum vector
+ int counter(0);
+
+ //Loop over events
+ std::cout<<"Analyzing "<<mychain->GetEntries()<<" events!"<<std::endl;
+ while (tr.Next()) {
+
+ if(counter%100==0) std::cout<<"Analyzing event "<<counter<<std::endl;
+ counter++;
+
+ //Loop over generated particles
+ for(size_t igen=0;igen<gen_status.GetSize();igen++){
+
+ auto charge = gen_charge[igen];
+ auto status = gen_status[igen];
+
+ //Require final-state, charged particle (no secondaries)
+ if(status==1 && fabs(charge) > 0.01 ){
+ gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
+ gen_vertex.SetXYZ(gen_vx[igen],gen_vy[igen],gen_vz[igen]);
+
+ //Fill eta histogram
+ h1a->Fill(gen_vec.Eta());
+ if( gen_vec.Pt()>0.2 ) h1a1->Fill(gen_vec.Eta());
+ if( gen_vec.Pt()>0.5 ) h1a2->Fill(gen_vec.Eta());
+ }
+ } //End loop over generated particles
+
+ //Loop over reconstructed real-seeded charged particles (copy of tracks with PID info)
+ size_t rec_mult = rec_type.GetSize();
+
+ for(size_t irec=0;irec<rec_mult;irec++){
+
+ rec_vec.SetXYZM(rec_px[irec],rec_py[irec],rec_pz[irec],rec_mass[irec]);
+
+ //Fill histograms
+ h1b->Fill(rec_vec.Eta());
+ if( rec_vec.Pt() > 0.2 ) h1b1->Fill(rec_vec.Eta());
+ if( rec_vec.Pt() > 0.5 ) h1b2->Fill(rec_vec.Eta());
+
+ } //End loop over reconstructed particles
+
+ //Loop over reconstructed truth-seeded charged particles (copy of tracks with PID info)
+ size_t rec_ts_mult = rec_ts_type.GetSize();
+
+ for(size_t irec=0;irec<rec_ts_mult;irec++){
+
+ rec_vec.SetXYZM(rec_ts_px[irec],rec_ts_py[irec],rec_ts_pz[irec],rec_ts_mass[irec]);
+
+ //Fill histograms
+ h1c->Fill(rec_vec.Eta());
+ if( rec_vec.Pt() > 0.2 ) h1c1->Fill(rec_vec.Eta());
+ if( rec_vec.Pt() > 0.5 ) h1c2->Fill(rec_vec.Eta());
+
+ } //End loop over reconstructed particles
+
+ // Loop over truth-seeded hit-based associations
+ for(size_t iassoc=0;iassoc<hit_assoc_weight.GetSize();iassoc++){
+
+ auto assoc_weight = hit_assoc_weight[iassoc];
+ h2a->Fill(assoc_weight);
+
+ } //End loop over hit associations
+
+ // Loop over truth-seeded hit-based associations
+ for(size_t iassoc=0;iassoc<hit_assoc_weight_ts.GetSize();iassoc++){
+
+ auto assoc_weight = hit_assoc_weight_ts[iassoc];
+ h2b->Fill(assoc_weight);
+
+ } //End loop over hit associations
+
+ } //End loop over events
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ ofile->Write(); // Write histograms to file
+ ofile->Close(); // Close output file
+
+}
\ No newline at end of file
diff --git a/benchmarks/tracking_performances_dis/analysis/trk_dis_plots.cxx b/benchmarks/tracking_performances_dis/analysis/trk_dis_plots.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7a69816136db8970a563384ba1b1fc0448fc1c57
--- /dev/null
+++ b/benchmarks/tracking_performances_dis/analysis/trk_dis_plots.cxx
@@ -0,0 +1,254 @@
+#include <fstream>
+#include <iostream>
+#include <string>
+
+#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TMath.h>
+#include <TStyle.h>
+#include <TLegend.h>
+
+#include "fmt/color.h"
+#include "fmt/core.h"
+
+#include "nlohmann/json.hpp"
+
+void trk_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;
+
+ TH1D* h1a = (TH1D*) file->Get("h1a");
+ TH1D* h1a1 = (TH1D*) file->Get("h1a1");
+ TH1D* h1a2 = (TH1D*) file->Get("h1a2");
+
+ TH1D* h1b = (TH1D*) file->Get("h1b");
+ TH1D* h1b1 = (TH1D*) file->Get("h1b1");
+ TH1D* h1b2 = (TH1D*) file->Get("h1b2");
+
+ TH1D* h1c = (TH1D*) file->Get("h1c");
+ TH1D* h1c1 = (TH1D*) file->Get("h1c1");
+ TH1D* h1c2 = (TH1D*) file->Get("h1c2");
+
+ TH1D* h2a = (TH1D*) file->Get("h2a");
+ TH1D* h2b = (TH1D*) file->Get("h2b");
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Make ratio histograms
+ TH1 *h1rb1 = new TH1D("h1rb1","",100,-4,4); //Real-seeded tracks (Pt > 200 MeV/c cut)
+ TH1 *h1rc1 = new TH1D("h1rc1","",100,-4,4); //Truth-seeded tracks (Pt > 200 MeV/c cut)
+ TH1 *h1rb2 = new TH1D("h1rb2","",100,-4,4); //Real-seeded tracks (Pt > 500 MeV/c cut)
+ TH1 *h1rc2 = new TH1D("h1rc2","",100,-4,4); //Truth-seeded tracks (Pt > 500 MeV/c cut)
+
+ h1rb1 = (TH1*) h1b1->Clone("h1rb1");
+ h1rb1->Divide(h1a1);
+ h1rb1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
+ h1rb1->GetXaxis()->SetTitle("#eta");h1rb1->GetXaxis()->CenterTitle();
+ h1rb1->GetYaxis()->SetTitle("Ratio");h1rb1->GetYaxis()->CenterTitle();
+
+ h1rc1 = (TH1*) h1c1->Clone("h1rc1");
+ h1rc1->Divide(h1a1);
+ h1rc1->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 200 MeV/c");
+ h1rc1->GetXaxis()->SetTitle("#eta");h1rc1->GetXaxis()->CenterTitle();
+ h1rc1->GetYaxis()->SetTitle("Ratio");h1rc1->GetYaxis()->CenterTitle();
+
+ h1rb2 = (TH1*) h1b2->Clone("h1rb2");
+ h1rb2->Divide(h1a2);
+ h1rb2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
+ h1rb2->GetXaxis()->SetTitle("#eta");h1rb2->GetXaxis()->CenterTitle();
+ h1rb2->GetYaxis()->SetTitle("Ratio");h1rb2->GetYaxis()->CenterTitle();
+
+ h1rc2 = (TH1*) h1c2->Clone("h1rc2");
+ h1rc2->Divide(h1a2);
+ h1rc2->SetTitle("Ratio of recontructed to generated particle counts: P_{T} > 500 MeV/c");
+ h1rc2->GetXaxis()->SetTitle("#eta");h1rc2->GetXaxis()->CenterTitle();
+ h1rc2->GetYaxis()->SetTitle("Ratio");h1rc2->GetYaxis()->CenterTitle();
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Make plots and save to PDF file
+
+ // Update Style
+ gStyle->SetPadGridX(0);
+ gStyle->SetPadGridY(0);
+ gStyle->SetOptStat(0);
+
+ std::cout<<"Making plots..."<<std::endl;
+
+ //Generated charged particles
+ TCanvas *c1a = new TCanvas("c1a");
+ h1a->Draw();
+ h1a1->Draw("same");
+ h1a2->Draw("same");
+
+ TLegend *leg1a = new TLegend(0.25,0.6,0.6,0.875);
+ leg1a->SetBorderSize(0);leg1a->SetFillStyle(0);
+ leg1a->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1a->AddEntry("h1a","All generated charged particles","l");
+ leg1a->AddEntry("h1a1","+ P_{T} > 200 MeV/c","l");
+ leg1a->AddEntry("h1a2","+ P_{T} > 500 MeV/c","l");
+ leg1a->Draw();
+
+ //Real-seeded tracks
+ TCanvas *c1b = new TCanvas("c1b");
+ h1b->Draw();
+ h1b1->Draw("same");
+ h1b2->Draw("same");
+
+ TLegend *leg1b = new TLegend(0.25,0.6,0.6,0.875);
+ leg1b->SetBorderSize(0);leg1b->SetFillStyle(0);
+ leg1b->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1b->AddEntry("h1b","All real-seeded tracks","l");
+ leg1b->AddEntry("h1b1","+ P_{T} > 200 MeV/c","l");
+ leg1b->AddEntry("h1b2","+ P_{T} > 500 MeV/c","l");
+ leg1b->Draw();
+
+ //Truth-seeded tracks
+ TCanvas *c1c = new TCanvas("c1c");
+ h1c->Draw();
+ h1c1->Draw("same");
+ h1c2->Draw("same");
+
+ TLegend *leg1c = new TLegend(0.25,0.6,0.6,0.875);
+ leg1c->SetBorderSize(0);leg1c->SetFillStyle(0);
+ leg1c->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1c->AddEntry("h1c","All truth-seeded tracks","l");
+ leg1c->AddEntry("h1c1","+ P_{T} > 200 MeV/c","l");
+ leg1c->AddEntry("h1c2","+ P_{T} > 500 MeV/c","l");
+ leg1c->Draw();
+
+ //Comparison 1
+ TCanvas *c1d = new TCanvas("c1d");
+ auto frame_d1 = c1d->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
+ frame_d1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_d1->GetXaxis()->CenterTitle();
+ h1a1->Draw("same");
+ h1b1->Draw("same");
+ h1c1->Draw("same");
+
+ TLegend *leg1d = new TLegend(0.15,0.675,0.6,0.875);
+ leg1d->SetBorderSize(0);leg1d->SetFillStyle(0);
+ leg1d->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1d->AddEntry("h1a1","Generated charged particles w/ P_{T} > 200 MeV/c","l");
+ leg1d->AddEntry("h1b1","Real-seeded tracks w/ P_{T} > 200 MeV/c","l");
+ leg1d->AddEntry("h1c1","Truth-seeded tracks w/ P_{T} > 200 MeV/c","l");
+ leg1d->Draw();
+
+ //Comparison 2a
+ TCanvas *c1e = new TCanvas("c1e");
+ auto frame_e1 = c1e->DrawFrame(-4,0,4,1.2*h1a1->GetMaximum());
+ frame_e1->GetXaxis()->SetTitle("#eta_{gen} or #eta_{rec}");frame_e1->GetXaxis()->CenterTitle();
+ h1a2->Draw("same");
+ h1b2->Draw("P same");
+
+ TLegend *leg1e = new TLegend(0.15,0.675,0.6,0.875);
+ leg1e->SetBorderSize(0);leg1e->SetFillStyle(0);
+ leg1e->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1e->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
+ leg1e->AddEntry("h1b2","Real-seeded tracks w/ P_{T} > 500 MeV/c","p");
+ leg1e->Draw();
+
+ //Comparison 2b
+ TCanvas *c1e1 = new TCanvas("c1e1");
+ frame_e1->Draw();
+ h1a2->Draw("same");
+ h1c2->Draw("P same");
+
+ TLegend *leg1e1 = new TLegend(0.15,0.675,0.6,0.875);
+ leg1e1->SetBorderSize(0);leg1e1->SetFillStyle(0);
+ leg1e1->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1e1->AddEntry("h1a2","Generated charged particles w/ P_{T} > 500 MeV/c","fl");
+ leg1e1->AddEntry("h1c2","Truth-seeded tracks w/ P_{T} > 500 MeV/c","p");
+ leg1e1->Draw();
+
+ //Comparison 1 ratio
+ TCanvas *c1f = new TCanvas("c1f");
+ h1rb1->Draw();
+ h1rc1->Draw("same");
+
+ TLegend *leg1f = new TLegend(0.575,0.25,0.875,0.45);
+ leg1f->SetBorderSize(0);leg1f->SetFillStyle(0);
+ leg1f->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1f->AddEntry("h1rb1","Real-seeded tracking","l");
+ leg1f->AddEntry("h1rc1","Truth-seeded tracking","l");
+ leg1f->Draw();
+
+ //Comparison 2 ratio
+ TCanvas *c1g = new TCanvas("c1g");
+ h1rb2->Draw();
+ h1rc2->Draw("same");
+
+ TLegend *leg1g = new TLegend(0.575,0.25,0.875,0.45);
+ leg1g->SetBorderSize(0);leg1g->SetFillStyle(0);
+ leg1g->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg1g->AddEntry("h1rb2","Real-seeded tracking","l");
+ leg1g->AddEntry("h1rc2","Truth-seeded tracking","l");
+ leg1g->Draw();
+
+ //Hit-based associations -- real-seeded tracks
+ TCanvas *c2a = new TCanvas("c2a");
+ c2a->SetLogy();
+ h2a->Draw();
+
+ TLegend *leg2a = new TLegend(0.25,0.6,0.6,0.875);
+ leg2a->SetBorderSize(0);leg2a->SetFillStyle(0);
+ leg2a->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg2a->Draw();
+
+ //Hit-based associations -- truth-seeded tracks
+ TCanvas *c2b = new TCanvas("c2b");
+ c2b->SetLogy();
+ h2b->Draw();
+
+ TLegend *leg2b = new TLegend(0.25,0.6,0.6,0.875);
+ leg2b->SetBorderSize(0);leg2a->SetFillStyle(0);
+ leg2b->SetHeader(Form("Pythia8: %dx%d GeV, Q^{2} > %d GeV^{2}",ebeam,pbeam,Q2_min));
+ leg2b->Draw();
+
+ //--------------------------------------------------------------------------------------------------------------------------------------------
+
+ // Print plots to pdf file
+ c1a->Print(fmt::format("{}.pdf[", output_prefix).c_str());
+ c1a->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1b->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1c->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1d->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1e->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1e1->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1f->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c1g->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c2a->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c2b->Print(fmt::format("{}.pdf", output_prefix).c_str());
+ c2b->Print(fmt::format("{}.pdf]", output_prefix).c_str());
+
+}
\ No newline at end of file
diff --git a/benchmarks/tracking_performances_dis/config.yml b/benchmarks/tracking_performances_dis/config.yml
new file mode 100644
index 0000000000000000000000000000000000000000..961df5492dfe5871cea6bd441f6d5b9d23c101f8
--- /dev/null
+++ b/benchmarks/tracking_performances_dis/config.yml
@@ -0,0 +1,35 @@
+sim:tracking_performances_dis:
+ extends: .det_benchmark
+ stage: simulate
+ script:
+ - snakemake --cache --cores 5 results/tracking_performances_dis/epic_craterlake_tracking_only/pythia8NCDIS_18x275_minQ2=1_combined_5/hists.root
+ retry:
+ max: 2
+ when:
+ - runner_system_failure
+
+bench:tracking_performances_dis:
+ extends: .det_benchmark
+ stage: benchmarks
+ needs:
+ - ["sim:tracking_performances_dis"]
+ script:
+ - snakemake --cores 1 trk_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
+
+collect_results:tracking_performances_dis:
+ extends: .det_benchmark
+ stage: collect
+ needs:
+ - "bench:tracking_performances_dis"
+ script:
+ - ls -lrht
+ - mv results{,_save}/ # move results directory out of the way to preserve it
+ - snakemake --cores 1 --delete-all-output trk_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