Tracking performance dis (#61)

.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:

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"

benchmarks/tracking_performances_dis/README.md

+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)
+
+

benchmarks/tracking_performances_dis/Snakefile

+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]}"
+
+

benchmarks/tracking_performances_dis/analysis/trk_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 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

benchmarks/tracking_performances_dis/analysis/trk_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 <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

benchmarks/tracking_performances_dis/config.yml

+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