From 28424bd660ad2562a0992223d805021f82189d7f Mon Sep 17 00:00:00 2001
From: Dmitry Kalinkin <dkalinkin@bnl.gov>
Date: Wed, 6 Sep 2023 20:31:55 +0000
Subject: [PATCH] add dis/jets benchmark
---
benchmarks/dis/analysis/jets.cxx | 558 +++++++++++++++++++++++++++++++
benchmarks/dis/dis.sh | 8 +
benchmarks/dis/env.sh | 1 +
3 files changed, 567 insertions(+)
create mode 100644 benchmarks/dis/analysis/jets.cxx
diff --git a/benchmarks/dis/analysis/jets.cxx b/benchmarks/dis/analysis/jets.cxx
new file mode 100644
index 00000000..66c53c7a
--- /dev/null
+++ b/benchmarks/dis/analysis/jets.cxx
@@ -0,0 +1,558 @@
+#include "common_bench/benchmark.h"
+#include "common_bench/mt.h"
+#include "common_bench/util.h"
+#include "common_bench/plot.h"
+
+#include "ROOT/RDataFrame.hxx"
+#include <TCanvas.h>
+#include <TH1D.h>
+#include <TVector3.h>
+
+#include "fmt/color.h"
+#include "fmt/core.h"
+
+#include "nlohmann/json.hpp"
+
+// Jet Benchmarks
+
+#define ETHRESH 5.0
+
+using ints = ROOT::VecOps::RVec<int>;
+using floats = ROOT::VecOps::RVec<float>;
+using vecs = ROOT::VecOps::RVec<TVector3>;
+
+int jets(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 results_path = config["results_path"];
+ const std::string plot_tag = config["plot_tag"];
+ const std::string test_tag = config["test_tag"];
+ const int ebeam = config["ebeam"];
+ const int pbeam = config["pbeam"];
+
+ fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
+ "Running jet analysis...\n");
+ fmt::print(" - Detector package: {}\n", detector);
+ fmt::print(" - input file: {}\n", rec_file);
+ fmt::print(" - test tag: {}\n", test_tag);
+ fmt::print(" - ebeam: {}\n", ebeam);
+ fmt::print(" - pbeam: {}\n", pbeam);
+
+ using namespace ROOT;
+
+ ROOT::EnableImplicitMT(kNumThreads);
+ ROOT::RDataFrame d("events", rec_file);
+
+ TFile output(fmt::format("{}/jets/{}.hist.root", results_path, plot_tag).c_str(), "RECREATE");
+
+ // Define Lambdas for Selection
+ auto getNJetsThresh = [](ints &type, floats &E){
+ ints mult;
+ int n=0;
+ for(auto i=0U; i<type.size(); i++)
+ {
+ if(type[i] == 0 && E[i] >= ETHRESH) n++;
+ }
+ mult.push_back(n);
+ return mult;
+ }; // Number of Jets Above Threshold
+
+ auto getNJets = [](ints &type, floats &E){
+ ints mult;
+ int n=0;
+ for(auto i=0U; i<type.size(); i++)
+ {
+ if(type[i] == 0) n++;
+ }
+ mult.push_back(n);
+ return mult;
+ };
+
+ auto getJetEThresh = [](ints &type, floats &E){
+ floats nrg;
+ for(auto i=0U; i<E.size(); i++)
+ {
+ if(type[i] == 0 && E[i] >= ETHRESH) nrg.push_back(E[i]);
+ }
+ return nrg;
+ }; // Jet Energy
+
+ auto getJetE = [](ints &type, floats &E){
+ floats nrg;
+ for(auto i=0U; i<E.size(); i++)
+ {
+ if(type[i] == 0) nrg.push_back(E[i]);
+ }
+ return nrg;
+ }; // Jet Energy
+
+ auto getJetVecThresh = [](ints &type, floats &E, floats &x, floats &y, floats &z){
+ vecs v;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ TVector3 tmp(x[i],y[i],z[i]);
+ if(type[i] == 0 && E[i] >= ETHRESH) v.push_back(tmp);
+ }
+ return v;
+ }; // Jet 3-Momentum
+
+ auto getJetVec = [](ints &type, floats &E, floats &x, floats &y, floats &z){
+ vecs v;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ TVector3 tmp(x[i],y[i],z[i]);
+ if(type[i] == 0)
+ {
+ //cout << i << " " << type[i] << " " << x[i] << " " << y[i] << " " << z[i] << " " << tmp.PseudoRapidity() << " " << tmp.Phi() << endl;
+ v.push_back(tmp);
+ }
+ }
+ return v;
+ }; // Jet 3-Momentum
+
+ auto getJetTrackDRThresh = [](ints &type, ints &pid, floats &E, floats &x, floats &y, floats &z){
+ floats dr;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
+ {
+ TVector3 jetMom(x[i],y[i],z[i]);
+ int jetID = pid[i];
+
+ for(auto j=0U; j<x.size(); j++) // First loop over tracks
+ {
+ if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
+ {
+ TVector3 trackA(x[j],y[j],z[j]);
+
+ if(j < x.size() - 1) // Overflow protection
+ {
+ for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
+ {
+ if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
+ {
+ TVector3 trackB(x[k],y[k],z[k]);
+
+ float dEta = trackA.PseudoRapidity() - trackB.PseudoRapidity();
+ float dPhi = TVector2::Phi_mpi_pi(trackA.Phi() - trackB.Phi());
+ float dR = std::sqrt(dEta*dEta + dPhi*dPhi);
+
+ dr.push_back(dR);
+ } // Track B
+ } // Second Loop
+ } // Overflow
+ } // Track A
+ } // First Loop
+ } // Jet Select
+ } // Jet Loop
+ return dr;
+ };
+
+ auto getJetTrackDR = [](ints &type, ints &pid, floats &x, floats &y, floats &z){
+ floats dr;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(type[i] == 0) // Look at jet
+ {
+ TVector3 jetMom(x[i],y[i],z[i]);
+ int jetID = pid[i];
+
+ for(auto j=0U; j<x.size(); j++) // First loop over tracks
+ {
+ if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
+ {
+ TVector3 trackA(x[j],y[j],z[j]);
+
+ if(j < x.size() - 1) // Overflow protection
+ {
+ for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
+ {
+ if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
+ {
+ TVector3 trackB(x[k],y[k],z[k]);
+
+ float dEta = trackA.PseudoRapidity() - trackB.PseudoRapidity();
+ float dPhi = TVector2::Phi_mpi_pi(trackA.Phi() - trackB.Phi());
+ float dR = std::sqrt(dEta*dEta + dPhi*dPhi);
+
+ dr.push_back(dR);
+ } // Track B
+ } // Second Loop
+ } // Overflow
+ } // Track A
+ } // First Loop
+ } // Jet Select
+ } // Jet Loop
+ return dr;
+ }; // Look at distance between particles in jet
+
+ auto getJetDupTrackThresh = [](ints &type, ints &pid, floats &E, floats &x, floats &y, floats &z){
+ ints dup;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
+ {
+ int hasDupTrack = 0;
+
+ TVector3 jetMom(x[i],y[i],z[i]);
+ int jetID = pid[i];
+
+ for(auto j=0U; j<x.size(); j++) // First loop over tracks
+ {
+ if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
+ {
+ TVector3 trackA(x[j],y[j],z[j]);
+
+ if(j < x.size() - 1) // Overflow protection
+ {
+ for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
+ {
+ if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
+ {
+ TVector3 trackB(x[k],y[k],z[k]);
+
+ float dEta = trackA.PseudoRapidity() - trackB.PseudoRapidity();
+ float dPhi = TVector2::Phi_mpi_pi(trackA.Phi() - trackB.Phi());
+ float dR = std::sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < 0.05) hasDupTrack = 1;
+ } // Track B
+ } // Second Loop
+ } // Overflow
+ } // Track A
+ } // First Loop
+
+ dup.push_back(hasDupTrack);
+ } // Jet Select
+ } // Jet Loop
+ return dup;
+ }; // Identify jets containing duplicate tracks
+
+ auto getPt = [](vecs &x){
+ floats pt;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ pt.push_back(x[i].Perp());
+ }
+ return pt;
+ }; // Jet Transverse Momentum
+
+ auto getEta = [](vecs &x){
+ floats eta;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ eta.push_back(x[i].PseudoRapidity());
+ }
+ return eta;
+ }; // Jet Psudorapidity
+
+ auto getPhi = [](vecs &x){
+ floats phi;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ phi.push_back(x[i].Phi());
+ }
+ return phi;
+ }; // Jet Phi
+
+ auto getMatchDeltaR = [](vecs &vecR, vecs &vecG){
+ floats deltaR;
+ for(auto i=0U; i<vecR.size(); i++)
+ {
+ float minDeltaR = 10000.0;
+ for(auto j=0U; j<vecG.size(); j++)
+ {
+ double dEta = vecR[i].PseudoRapidity() - vecG[j].PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(vecR[i].Phi() - vecG[j].Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < minDeltaR) minDeltaR = dR;
+ }
+ deltaR.push_back(minDeltaR);
+ }
+ return deltaR;
+ };
+
+ auto getMatchE = [](vecs &vecR, vecs &vecG, floats &E){
+ floats nrg;
+ for(auto i=0U; i<vecR.size(); i++)
+ {
+ float minDeltaR = 10000.0;
+ float minE = -1.0;
+ for(auto j=0U; j<vecG.size(); j++)
+ {
+ double dEta = vecR[i].PseudoRapidity() - vecG[j].PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(vecR[i].Phi() - vecG[j].Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < minDeltaR)
+ {
+ minDeltaR = dR;
+ minE = E[j];
+ }
+ }
+ nrg.push_back(minE);
+ }
+ return nrg;
+ };
+
+ auto getMatchEta = [](vecs &vecR, vecs &vecG){
+ floats eta;
+ for(auto i=0U; i<vecR.size(); i++)
+ {
+ float minDeltaR = 10000.0;
+ float minEta = -100.0;
+ for(auto j=0U; j<vecG.size(); j++)
+ {
+ double dEta = vecR[i].PseudoRapidity() - vecG[j].PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(vecR[i].Phi() - vecG[j].Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < minDeltaR)
+ {
+ minDeltaR = dR;
+ minEta = vecG[j].PseudoRapidity();
+ }
+ }
+ eta.push_back(minEta);
+ }
+ return eta;
+ };
+
+ auto getMatchPhi = [](vecs &vecR, vecs &vecG){
+ floats phi;
+ for(auto i=0U; i<vecR.size(); i++)
+ {
+ float minDeltaR = 10000.0;
+ float minPhi = -100.0;
+ for(auto j=0U; j<vecG.size(); j++)
+ {
+ double dEta = vecR[i].PseudoRapidity() - vecG[j].PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(vecR[i].Phi() - vecG[j].Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < minDeltaR)
+ {
+ minDeltaR = dR;
+ minPhi = vecG[j].Phi();
+ }
+ }
+ phi.push_back(minPhi);
+ }
+ return phi;
+ };
+
+ auto calcRes = [](floats &x, floats &y){
+ floats res;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ res.push_back((x[i]-y[i])/y[i]);
+ }
+ return res;
+ };
+
+ auto filterDR = [](floats &x, floats &dr){
+ floats e;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(dr[i] < 0.25) e.push_back(x[i]);
+ }
+ return e;
+ };
+
+ auto filterDRI = [](ints &x, floats &dr){
+ ints e;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(dr[i] < 0.25) e.push_back(x[i]);
+ }
+ return e;
+ };
+
+ auto filterDup = [](floats &x, ints &dup){
+ floats e;
+ for(auto i=0U; i<x.size(); i++)
+ {
+ if(dup[i] == 0) e.push_back(x[i]);
+ }
+ return e;
+ };
+
+ // Define RDataFrame Columns
+ auto d1 = d.Define("jetMultThresh",getNJetsThresh,{"ReconstructedJets.type","ReconstructedJets.energy"})
+ .Define("jetEThresh",getJetEThresh,{"ReconstructedJets.type","ReconstructedJets.energy"})
+ .Define("jetTrackDRThresh",getJetTrackDRThresh,{"ReconstructedJets.type","ReconstructedJets.PDG","ReconstructedJets.energy","ReconstructedJets.momentum.x","ReconstructedJets.momentum.y","ReconstructedJets.momentum.z"})
+ .Define("jetDupTrackThresh",getJetDupTrackThresh,{"ReconstructedJets.type","ReconstructedJets.PDG","ReconstructedJets.energy","ReconstructedJets.momentum.x","ReconstructedJets.momentum.y","ReconstructedJets.momentum.z"})
+ .Define("jetVecThresh",getJetVecThresh,{"ReconstructedJets.type","ReconstructedJets.energy","ReconstructedJets.momentum.x","ReconstructedJets.momentum.y","ReconstructedJets.momentum.z"})
+ .Define("jetPtThresh",getPt,{"jetVecThresh"})
+ .Define("jetEtaThresh",getEta,{"jetVecThresh"})
+ .Define("jetPhiThresh",getPhi,{"jetVecThresh"})
+ .Define("genJetMult",getNJets,{"GeneratedJets.type","GeneratedJets.energy"})
+ .Define("genJetE",getJetE,{"GeneratedJets.type","GeneratedJets.energy"})
+ .Define("genJetTrackDR",getJetTrackDR,{"GeneratedJets.type","GeneratedJets.PDG","GeneratedJets.momentum.x","GeneratedJets.momentum.y","GeneratedJets.momentum.z"})
+ .Define("genJetVec",getJetVec,{"GeneratedJets.type","GeneratedJets.energy","GeneratedJets.momentum.x","GeneratedJets.momentum.y","GeneratedJets.momentum.z"})
+ .Define("genJetPt",getPt,{"genJetVec"})
+ .Define("genJetEta",getEta,{"genJetVec"})
+ .Define("genJetPhi",getPhi,{"genJetVec"})
+ .Define("matchJetDeltaR",getMatchDeltaR,{"jetVecThresh","genJetVec"})
+ .Define("matchJetE",getMatchE,{"jetVecThresh","genJetVec","genJetE"})
+ .Define("matchJetEta",getMatchEta,{"jetVecThresh","genJetVec"})
+ .Define("matchJetPhi",getMatchPhi,{"jetVecThresh","genJetVec"})
+ .Define("matchJetResE",calcRes,{"jetEThresh","matchJetE"})
+ .Define("jetEThreshDR",filterDR,{"jetEThresh","matchJetDeltaR"})
+ .Define("matchJetEDR",filterDR,{"matchJetE","matchJetDeltaR"})
+ .Define("jetEtaThreshDR",filterDR,{"jetEtaThresh","matchJetDeltaR"})
+ .Define("jetPhiThreshDR",filterDR,{"jetPhiThresh","matchJetDeltaR"})
+ .Define("matchJetEtaDR",filterDR,{"matchJetEta","matchJetDeltaR"})
+ .Define("matchJetPhiDR",filterDR,{"matchJetPhi","matchJetDeltaR"})
+ .Define("matchJetResEDR",calcRes,{"jetEThreshDR","matchJetEDR"})
+ .Define("jetDupTrackThreshDR",filterDRI,{"jetDupTrackThresh","matchJetDeltaR"})
+ .Define("jetEThreshDRDup",filterDup,{"jetEThreshDR","jetDupTrackThreshDR"})
+ .Define("matchJetEDRDup",filterDup,{"matchJetEDR","jetDupTrackThreshDR"})
+ .Define("jetEtaThreshDRDup",filterDup,{"jetEtaThreshDR","jetDupTrackThreshDR"})
+ .Define("jetPhiThreshDRDup",filterDup,{"jetPhiThreshDR","jetDupTrackThreshDR"})
+ .Define("matchJetEtaDRDup",filterDup,{"matchJetEtaDR","jetDupTrackThreshDR"})
+ .Define("matchJetPhiDRDup",filterDup,{"matchJetPhiDR","jetDupTrackThreshDR"})
+ .Define("matchJetResEDRDup",calcRes,{"jetEThreshDRDup","matchJetEDRDup"});
+
+ // Book Histograms
+ auto h_type = d.Histo1D({"h_type",";Type",3, 0.,3.}, "ReconstructedJets.type");
+ auto h_nJetsThresh = d1.Histo1D({"h_nJetsThresh","Number of Reco Jets (E >= 5 GeV);Num Jets",20,0.,20.}, "jetMultThresh");
+ auto h_jetEnergyThresh = d1.Histo1D({"h_jetEnergyThresh","Reco Jet Energy (E >= 5 GeV);Energy [GeV]",300,0.,300.}, "jetEThresh");
+ auto h_jetTrackDRThresh = d1.Histo1D({"h_jetTrackDRThresh","Reco Jet: Pairwise Distance Between Tracks (E >= 5 GeV);Delta R",1000,0.,10.}, "jetTrackDRThresh");
+ auto h_jetDupTrackThresh = d1.Histo1D({"h_jetDupTrackThresh","Reco Jets with Duplicate Tracks (E >= 5 GeV);0 = No Duplicate 1 = Duplicate",5,0.,5.}, "jetDupTrackThresh");
+ auto h_jetPtThresh = d1.Histo1D({"h_jetPtThresh","Reco Jet Pt (E >= 5 GeV);Jet Pt [GeV]",100,0.,50.}, "jetPtThresh");
+ auto h_jetEtaThresh = d1.Histo1D({"h_jetEtaThresh","Reco Jet Eta (E >= 5 GeV);Eta",100,-5.,5.}, "jetEtaThresh");
+ auto h_jetPhiThresh = d1.Histo1D({"h_jetPhiThresh","Reco Jet Phi (E >= 5 GeV);Phi [Rad]",100,-TMath::Pi(),TMath::Pi()}, "jetPhiThresh");
+ auto h_jetPhiVsEtaThresh = d1.Histo2D({"h_jetPhiVsEtaThresh","Reco Jet Phi Vs Eta (E >= 5 GeV);Jet Eta;Jet Phi",100,-5.,5.,100,-TMath::Pi(),TMath::Pi()}, "jetEtaThresh", "jetPhiThresh");
+ auto h_jetEVsEtaThresh = d1.Histo2D({"h_jetEVsEtaThresh","Reco Jet E Vs Eta (E >= 5 GeV);Jet Eta;Jet Energy [GeV]",100,-5.,5.,300,0.,300.}, "jetEtaThresh", "jetEThresh");
+
+ auto h_nGenJets = d1.Histo1D({"h_nGenJets","Number of Truth Jets;Num Jets",20,0.,20.}, "genJetMult");
+ auto h_genJetEnergy = d1.Histo1D({"h_genJetEnergy","Truth Jet Energy;Energy [GeV]",300,0.,300.}, "genJetE");
+ auto h_genJetTrackDR = d1.Histo1D({"h_genJetTrackDR","Truth Jet: Pairwise Distance Between Charged Particles;Delta R",1000,0.,10.}, "genJetTrackDR");
+ auto h_genJetPt = d1.Histo1D({"h_genJetPt","Truth Jet Pt;Jet Pt [GeV]",100,0.,50.}, {"genJetPt"});
+ auto h_genJetEta = d1.Histo1D({"h_genJetEta","Truth Jet Eta;Eta",100,-5.,5.}, "genJetEta");
+ auto h_genJetPhi = d1.Histo1D({"h_genJetPhi","Truth Jet Phi;Phi [Rad]",100,-TMath::Pi(),TMath::Pi()}, {"genJetPhi"});
+ auto h_genJetPhiVsEta = d1.Histo2D({"h_genJetPhiVsEta","Truth Jet Phi Vs Eta;Jet Eta;Jet Phi",100,-5.,5.,100,-TMath::Pi(),TMath::Pi()}, "genJetEta", "genJetPhi");
+ auto h_genJetEVsEta = d1.Histo2D({"h_genJetEVsEta","Truth Jet E Vs Eta;Jet Eta;Jet Energy [GeV]",100,-5.,5.,300,0.,300.}, "genJetEta", "genJetE");
+
+ auto h_jetMatchDeltaR = d1.Histo1D({"h_jetMatchDeltaR","Distance Between Closest Reco-Truth Jet;Delta R",1000,0.,5.}, "matchJetDeltaR");
+ auto h_jetMatchTruthVsRecoE = d1.Histo2D({"h_jetMatchTruthVsRecoE","Truth Vs Reco Jet Energy;Reco Energy [GeV];Truth Energy [GeV]",300,0.,300.,300,0.,300.}, "jetEThresh","matchJetE");
+ auto h_jetMatchTruthVsRecoEta = d1.Histo2D({"h_jetMatchTruthVsRecoEta","Truth Vs Reco Jet Eta;Reco Eta;Truth Eta",100,-5.,5.,100,-5.,5.}, "jetEtaThresh","matchJetEta");
+ auto h_jetMatchTruthVsRecoPhi = d1.Histo2D({"h_jetMatchTruthVsRecoPhi","Truth Vs Reco Jet Phi;Reco Phi;Truth Phi",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi()}, "jetPhiThresh","matchJetPhi");
+ auto h_jetMatchERes = d1.Histo1D({"h_jetMatchERes","(Reco - Truth)/Truth Jet Energy;Energy Diff",200,-10.,10.}, "matchJetResE");
+ auto h_jetMatchEResVsE = d1.Histo2D({"h_jetMatchEResVsE","(Reco - Truth)/Truth Jet Energy Vs Reco Energy;Energy [GeV];Energy Diff",300,0.,300,2000,-10.,10.}, "jetEThresh","matchJetResE");
+ auto h_jetMatchEResVsEta = d1.Histo2D({"h_jetMatchEResVsEta","(Reco - Truth)/Truth Jet Energy Vs Reco Eta;Reco Jet Eta;Energy Diff",100,-5.,5.,2000,-10.,10.}, "jetEtaThresh", "matchJetResE");
+
+ auto h_jetMatchTruthVsRecoEDR = d1.Histo2D({"h_jetMatchTruthVsRecoEDR","Truth Vs Reco Jet Energy (Delta R < 0.25);Reco Energy [GeV];Truth Energy [GeV]",300,0.,300.,300,0.,300.}, "jetEThreshDR","matchJetEDR");
+ auto h_jetMatchTruthVsRecoEtaDR = d1.Histo2D({"h_jetMatchTruthVsRecoEtaDR","Truth Vs Reco Jet Eta (Delta R < 0.25);Reco Eta;Truth Eta",100,-5.,5.,100,-5.,5.}, "jetEtaThreshDR","matchJetEtaDR");
+ auto h_jetMatchTruthVsRecoPhiDR = d1.Histo2D({"h_jetMatchTruthVsRecoPhiDR","Truth Vs Reco Jet Phi (Delta R < 0.25);Reco Phi;Truth Phi",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi()}, "jetPhiThreshDR","matchJetPhiDR");
+ auto h_jetMatchEResDR = d1.Histo1D({"h_jetMatchEResDR","(Reco - Truth)/Truth Jet Energy (Delta R < 0.25);Energy Diff",200,-10.,10.}, "matchJetResEDR");
+ auto h_jetMatchEResVsEDR = d1.Histo2D({"h_jetMatchEResVsEDR","(Reco - Truth)/Truth Jet Energy Vs Reco Energy (Delta R < 0.25);Energy [GeV];Energy Diff",300,0.,300,2000,-10.,10.}, "jetEThreshDR","matchJetResEDR");
+ auto h_jetMatchEResVsEtaDR = d1.Histo2D({"h_jetMatchEResVsEtaDR","(Reco - Truth)/Truth Jet Energy Vs Reco Eta (Delta R < 0.25);Reco Jet Eta;Energy Diff",100,-5.,5.,2000,-10.,10.}, "jetEtaThreshDR", "matchJetResEDR");
+
+ auto h_jetMatchTruthVsRecoEDRDup = d1.Histo2D({"h_jetMatchTruthVsRecoEDRDup","Truth Vs Reco Jet Energy (Delta R < 0.25 No Duplicate);Reco Energy [GeV];Truth Energy [GeV]",300,0.,300.,300,0.,300.}, "jetEThreshDRDup","matchJetEDRDup");
+ auto h_jetMatchTruthVsRecoEtaDRDup = d1.Histo2D({"h_jetMatchTruthVsRecoEtaDRDup","Truth Vs Reco Jet Eta (Delta R < 0.25 No Duplicate);Reco Eta;Truth Eta",100,-5.,5.,100,-5.,5.}, "jetEtaThreshDRDup","matchJetEtaDRDup");
+ auto h_jetMatchTruthVsRecoPhiDRDup = d1.Histo2D({"h_jetMatchTruthVsRecoPhiDRDup","Truth Vs Reco Jet Phi (Delta R < 0.25 No Duplicate);Reco Phi;Truth Phi",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi()}, "jetPhiThreshDRDup","matchJetPhiDRDup");
+ auto h_jetMatchEResDRDup = d1.Histo1D({"h_jetMatchEResDRDup","(Reco - Truth)/Truth Jet Energy (Delta R < 0.25 No Duplicate);Energy Diff",200,-10.,10.}, "matchJetResEDRDup");
+ auto h_jetMatchEResVsEDRDup = d1.Histo2D({"h_jetMatchEResVsEDRDup","(Reco - Truth)/Truth Jet Energy Vs Reco Energy (Delta R < 0.25 No Duplicate);Energy [GeV];Energy Diff",300,0.,300,2000,-10.,10.}, "jetEThreshDRDup","matchJetResEDRDup");
+ auto h_jetMatchEResVsEtaDRDup = d1.Histo2D({"h_jetMatchEResVsEtaDRDup","(Reco - Truth)/Truth Jet Energy Vs Reco Eta (Delta R < 0.25 No Duplicate);Reco Jet Eta;Energy Diff",100,-5.,5.,2000,-10.,10.}, "jetEtaThreshDRDup", "matchJetResEDRDup");
+
+ // Render histograms
+ { TCanvas c; h_type->Draw(); c.Print(fmt::format("{}/jets/{}_type.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_nJetsThresh->Draw(); c.Print(fmt::format("{}/jets/{}_nJetsThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetEnergyThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetEnergyThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetTrackDRThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetTrackDRThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetDupTrackThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetDupTrackThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetPtThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetPtThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetEtaThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetEtaThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetPhiThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetPhiThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetPhiVsEtaThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetPhiVsEtaThresh.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetEVsEtaThresh->Draw(); c.Print(fmt::format("{}/jets/{}_jetEVsEtaThresh.png", results_path, plot_tag).c_str()); }
+
+ { TCanvas c; h_nGenJets->Draw(); c.Print(fmt::format("{}/jets/{}_nGenJets.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetEnergy->Draw(); c.Print(fmt::format("{}/jets/{}_genJetEnergy.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetTrackDR->Draw(); c.Print(fmt::format("{}/jets/{}_genJetTrackDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetPt->Draw(); c.Print(fmt::format("{}/jets/{}_genJetPt.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetEta->Draw(); c.Print(fmt::format("{}/jets/{}_genJetEta.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetPhi->Draw(); c.Print(fmt::format("{}/jets/{}_genJetPhi.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetPhiVsEta->Draw(); c.Print(fmt::format("{}/jets/{}_genJetPhiVsEta.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_genJetEVsEta->Draw(); c.Print(fmt::format("{}/jets/{}_genJetEVsEta.png", results_path, plot_tag).c_str()); }
+
+ { TCanvas c; h_jetMatchDeltaR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchDeltaR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoE->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoE.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoEta->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoEta.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoPhi->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoPhi.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchERes->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchERes.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsE->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsE.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsEta->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsEta.png", results_path, plot_tag).c_str()); }
+
+ { TCanvas c; h_jetMatchTruthVsRecoEDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoEDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoEtaDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoEtaDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoPhiDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoPhiDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsEDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsEDR.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsEtaDR->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsEtaDR.png", results_path, plot_tag).c_str()); }
+
+ { TCanvas c; h_jetMatchTruthVsRecoEDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoEDRDup.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoEtaDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoEtaDRDup.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchTruthVsRecoPhiDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchTruthVsRecoPhiDRDup.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResDRDup.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsEDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsEDRDup.png", results_path, plot_tag).c_str()); }
+ { TCanvas c; h_jetMatchEResVsEtaDRDup->Draw(); c.Print(fmt::format("{}/jets/{}_jetMatchEResVsEtaDRDup.png", results_path, plot_tag).c_str()); }
+
+ // Write Histograms
+ h_type->Write();
+ h_nJetsThresh->Write();
+ h_jetEnergyThresh->Write();
+ h_jetTrackDRThresh->Write();
+ h_jetDupTrackThresh->Write();
+ h_jetPtThresh->Write();
+ h_jetEtaThresh->Write();
+ h_jetPhiThresh->Write();
+ h_jetPhiVsEtaThresh->Write();
+ h_jetEVsEtaThresh->Write();
+
+ h_nGenJets->Write();
+ h_genJetEnergy->Write();
+ h_genJetTrackDR->Write();
+ h_genJetPt->Write();
+ h_genJetEta->Write();
+ h_genJetPhi->Write();
+ h_genJetPhiVsEta->Write();
+ h_genJetEVsEta->Write();
+
+ h_jetMatchDeltaR->Write();
+ h_jetMatchTruthVsRecoE->Write();
+ h_jetMatchTruthVsRecoEta->Write();
+ h_jetMatchTruthVsRecoPhi->Write();
+ h_jetMatchERes->Write();
+ h_jetMatchEResVsE->Write();
+ h_jetMatchEResVsEta->Write();
+
+ h_jetMatchTruthVsRecoEDR->Write();
+ h_jetMatchTruthVsRecoEtaDR->Write();
+ h_jetMatchTruthVsRecoPhiDR->Write();
+ h_jetMatchEResDR->Write();
+ h_jetMatchEResVsEDR->Write();
+ h_jetMatchEResVsEtaDR->Write();
+
+ h_jetMatchTruthVsRecoEDRDup->Write();
+ h_jetMatchTruthVsRecoEtaDRDup->Write();
+ h_jetMatchTruthVsRecoPhiDRDup->Write();
+ h_jetMatchEResDRDup->Write();
+ h_jetMatchEResVsEDRDup->Write();
+ h_jetMatchEResVsEtaDRDup->Write();
+
+ // Write output and close file
+ output.Write();
+ output.Close();
+
+ return 0;
+}
diff --git a/benchmarks/dis/dis.sh b/benchmarks/dis/dis.sh
index e46b7ae2..f1a551d9 100755
--- a/benchmarks/dis/dis.sh
+++ b/benchmarks/dis/dis.sh
@@ -108,6 +108,8 @@ cat << EOF > ${CONFIG}
"rec_file": "${REC_FILE}",
"detector": "${DETECTOR}",
"output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+ "results_path": "${RESULTS_PATH}",
+ "plot_tag": "${PLOT_TAG}",
"ebeam": ${EBEAM},
"pbeam": ${PBEAM},
"minq2": ${MINQ2},
@@ -133,6 +135,12 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
+root -b -q "benchmarks/dis/analysis/jets.cxx+g(\"${CONFIG}\")"
+if [[ "$?" -ne "0" ]] ; then
+ echo "ERROR running dis_electron script"
+ exit 1
+fi
+
## =============================================================================
## Step 5: finalize
diff --git a/benchmarks/dis/env.sh b/benchmarks/dis/env.sh
index 1afe4f05..50f6bb04 100644
--- a/benchmarks/dis/env.sh
+++ b/benchmarks/dis/env.sh
@@ -49,6 +49,7 @@ echo "TMP_PATH: ${TMP_PATH}"
RESULTS_PATH="results/${BENCHMARK_TAG}/${BEAM_TAG}/minQ2=${MINQ2}"
mkdir -p ${RESULTS_PATH}
mkdir -p "${RESULTS_PATH}/truth_reconstruction"
+mkdir -p "${RESULTS_PATH}/jets"
export RESULTS_PATH=`realpath ${RESULTS_PATH}`
echo "RESULTS_PATH: ${RESULTS_PATH}"
--
GitLab