diff --git a/benchmarks/dis/analysis/jets.cxx b/benchmarks/dis/analysis/jets.cxx
index 66c53c7aa8258179b16edc2a19233d0f4e43d4f8..2f696813a58efc937b975a2c6c090ea3172074a1 100644
--- a/benchmarks/dis/analysis/jets.cxx
+++ b/benchmarks/dis/analysis/jets.cxx
@@ -3,9 +3,15 @@
#include "common_bench/util.h"
#include "common_bench/plot.h"
-#include "ROOT/RDataFrame.hxx"
#include <TCanvas.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TGraph.h>
#include <TH1D.h>
+#include <TH2D.h>
+#include <TTreeReader.h>
+#include <TTreeReaderArray.h>
+#include <TLegend.h>
#include <TVector3.h>
#include "fmt/color.h"
@@ -14,12 +20,9 @@
#include "nlohmann/json.hpp"
// Jet Benchmarks
+// Author: B. Page (bpage@bnl.gov)
-#define ETHRESH 5.0
-
-using ints = ROOT::VecOps::RVec<int>;
-using floats = ROOT::VecOps::RVec<float>;
-using vecs = ROOT::VecOps::RVec<TVector3>;
+// To run: snakemake -c1 results/epic_craterlake/dis/10on100/minQ2=1/jets
int jets(const std::string& config_name)
{
@@ -44,515 +47,1268 @@ int jets(const std::string& config_name)
fmt::print(" - ebeam: {}\n", ebeam);
fmt::print(" - pbeam: {}\n", pbeam);
- using namespace ROOT;
+ const bool PRINT = true;
- ROOT::EnableImplicitMT(kNumThreads);
- ROOT::RDataFrame d("events", rec_file);
+ // Input
+ TChain *mychain = new TChain("events");
+ mychain->Add(rec_file.c_str());
- TFile output(fmt::format("{}/jets/{}.hist.root", results_path, plot_tag).c_str(), "RECREATE");
+ const int seabornRed = TColor::GetColor(213, 94, 0);
- // 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;
- };
+// Seaborn Green: #009E73 -> (0, 158, 115)
+const int seabornGreen = TColor::GetColor(0, 158, 115);
- 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
+// Seaborn Blue: #56B4E9 -> (86, 180, 233)
+const int seabornBlue = TColor::GetColor(100, 149, 237);
- 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
+ // Output
+ //TFile *ofile = TFile::Open("test_24-05-0.hist.root","RECREATE");
- 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
+ // TTreeReader
+ TTreeReader tree_reader(mychain);
- 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
+ // Set Delta R Cut
+ float DELTARCUT = 0.05;
+
+ // Reco Jets
+ TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedJets.type"};
+ TTreeReaderArray<float> recoNRG = {tree_reader, "ReconstructedChargedJets.energy"};
+ TTreeReaderArray<int> recoPDG = {tree_reader, "ReconstructedChargedJets.PDG"};
+ TTreeReaderArray<float> recoMomX = {tree_reader, "ReconstructedChargedJets.momentum.x"};
+ TTreeReaderArray<float> recoMomY = {tree_reader, "ReconstructedChargedJets.momentum.y"};
+ TTreeReaderArray<float> recoMomZ = {tree_reader, "ReconstructedChargedJets.momentum.z"};
+ TTreeReaderArray<float> recoM = {tree_reader, "ReconstructedChargedJets.mass"};
+ TTreeReaderArray<unsigned int> partsBegin = {tree_reader, "ReconstructedChargedJets.particles_begin"};
+ TTreeReaderArray<unsigned int> partsEnd = {tree_reader, "ReconstructedChargedJets.particles_end"};
+
+ TTreeReaderArray<int> recoClustIndex = {tree_reader, "_ReconstructedChargedJets_clusters.index"};
+ TTreeReaderArray<int> recoTrackIndex = {tree_reader, "_ReconstructedChargedJets_tracks.index"};
+ TTreeReaderArray<int> recoPartIndex = {tree_reader, "_ReconstructedChargedJets_particles.index"};
+
+ // Reconstructed Particles
+ TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
+ TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
+ TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
+ TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
+ TTreeReaderArray<int> recoPartPDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
+ TTreeReaderArray<float> recoPartNRG = {tree_reader, "ReconstructedChargedParticles.energy"};
+
+ // Generated Jets
+ TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
+ TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
+ TTreeReaderArray<int> genPDG = {tree_reader, "GeneratedChargedJets.PDG"};
+ TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
+ TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
+ TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
+ TTreeReaderArray<float> genM = {tree_reader, "GeneratedChargedJets.mass"};
+ TTreeReaderArray<unsigned int> genPartsBegin = {tree_reader, "GeneratedChargedJets.particles_begin"};
+ TTreeReaderArray<unsigned int> genPartsEnd = {tree_reader, "GeneratedChargedJets.particles_end"};
+
+ TTreeReaderArray<int> genPartIndex = {tree_reader, "_GeneratedChargedJets_particles.index"};
+ //TTreeReaderArray<int> genChargedIndex = {tree_reader, "GeneratedChargedParticles_objIdx.index"};
+
+ // MC
+ //TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
+ TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
+ TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
+ TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
+ TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
+ TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
+
+ // Define Histograms
+ TH1D *counter = new TH1D("counter","",10,0.,10.);
+
+
+ // Reco
+ TH1D *numRecoChargedJetsECutHist = new TH1D("numRecoChargedJetsECut","",20,0.,20.);
+ TH1D *recoChargedJetEHist = new TH1D("recoChargedJetE","",300,0.,300.);
+ TH1D *recoChargedJetEtaECutHist = new TH1D("recoChargedJetEtaECut","",60,-3.,3.);
+ TH2D *recoChargedJetEvsEtaHist = new TH2D("recoChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
+ TH2D *recoChargedJetPhiVsEtaECutHist = new TH2D("recoChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numRecoChargedJetsECutNoElecHist = new TH1D("numRecoChargedJetsECutNoElec","",20,0.,20.);
+ TH1D *recoChargedJetENoElecHist = new TH1D("recoChargedJetENoElec","",300,0.,300.);
+ TH1D *recoChargedJetEtaECutNoElecHist = new TH1D("recoChargedJetEtaECutNoElec","",60,-3.,3.);
+ TH2D *recoChargedJetEvsEtaNoElecHist = new TH2D("recoChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
+ TH2D *recoChargedJetPhiVsEtaECutNoElecHist = new TH2D("recoChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numRecoChargedJetPartsHist = new TH1D("numRecoChargedJetParts","",20,0.,20.);
+ TH1D *recoChargedJetPartEHist = new TH1D("recoChargedJetPartE","",500,0.,100.);
+ TH1D *recoChargedJetPartEtaHist = new TH1D("recoChargedJetPartEta","",80,-4.,4.);
+ TH2D *recoChargedJetPartEvsEtaHist = new TH2D("recoChargedJetPartEvsEta","",80,-4.,4.,500,0.,100.);
+ TH2D *recoChargedJetPartPhiVsEtaHist = new TH2D("recoChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numRecoChargedJetPartsNoElecHist = new TH1D("numRecoChargedJetPartsNoElec","",20,0.,20.);
+ TH1D *recoChargedJetPartENoElecHist = new TH1D("recoChargedJetPartENoElec","",500,0.,100.);
+ TH1D *recoChargedJetPartEtaNoElecHist = new TH1D("recoChargedJetPartEtaNoElec","",80,-4.,4.);
+ TH2D *recoChargedJetPartEvsEtaNoElecHist = new TH2D("recoChargedJetPartEvsEtaNoElec","",80,-4.,4.,500,0.,100.);
+ TH2D *recoChargedJetPartPhiVsEtaNoElecHist = new TH2D("recoChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *recoChargedJetPartPairwiseDeltaRHist = new TH1D("recoChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
+
+ // Gen
+ TH1D *numGenChargedJetsECutHist = new TH1D("numGenChargedJetsECut","",20,0.,20.);
+ TH1D *genChargedJetEHist = new TH1D("genChargedJetE","",300,0.,300.);
+ TH1D *genChargedJetEtaECutHist = new TH1D("genChargedJetEtaECut","",60,-3.,3.);
+ TH2D *genChargedJetEvsEtaHist = new TH2D("genChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
+ TH2D *genChargedJetPhiVsEtaECutHist = new TH2D("genChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numGenChargedJetsECutNoElecHist = new TH1D("numGenChargedJetsECutNoElec","",20,0.,20.);
+ TH1D *genChargedJetENoElecHist = new TH1D("genChargedJetENoElec","",300,0.,300.);
+ TH1D *genChargedJetEtaECutNoElecHist = new TH1D("genChargedJetEtaECutNoElec","",60,-3.,3.);
+ TH2D *genChargedJetEvsEtaNoElecHist = new TH2D("genChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
+ TH2D *genChargedJetPhiVsEtaECutNoElecHist = new TH2D("genChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numGenChargedJetPartsHist = new TH1D("numGenChargedJetParts","",20,0.,20.);
+ TH1D *genChargedJetPartEHist = new TH1D("genChargedJetPartE","",500,0.,100.);
+ TH1D *genChargedJetPartEtaHist = new TH1D("genChargedJetPartEta","",80,-4.,4.);
+ TH2D *genChargedJetPartEvsEtaHist = new TH2D("genChargedJetPartEvsEta","",80,-4.,4.,500,0.,100.);
+ TH2D *genChargedJetPartPhiVsEtaHist = new TH2D("genChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *numGenChargedJetPartsNoElecHist = new TH1D("numGenChargedJetPartsNoElec","",20,0.,20.);
+ TH1D *genChargedJetPartENoElecHist = new TH1D("genChargedJetPartENoElec","",500,0.,100.);
+ TH1D *genChargedJetPartEtaNoElecHist = new TH1D("genChargedJetPartEtaNoElec","",80,-4.,4.);
+ TH2D *genChargedJetPartEvsEtaNoElecHist = new TH2D("genChargedJetPartEvsEtaNoElec","",80,-4.,4.,500,0.,100.);
+ TH2D *genChargedJetPartPhiVsEtaNoElecHist = new TH2D("genChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
+
+ TH1D *genChargedJetPartPairwiseDeltaRHist = new TH1D("genChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
+
+ // Matched
+ TH1D *matchJetDeltaRHist = new TH1D("matchJetDeltaR","",5000,0.,5.);
+ TH1D *matchJetDeltaRBackHist = new TH1D("matchJetDeltaRBack","",5000,0.,5.);
+ TH2D *recoVsGenChargedJetEtaHist = new TH2D("recoVsGenChargedJetEta","",80,-4.,4.,80,-4.,4.);
+ TH2D *recoVsGenChargedJetPhiHist = new TH2D("recoVsGenChargedJetPhi","",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
+ TH2D *recoVsGenChargedJetEHist = new TH2D("recoVsGenChargedJetE","",100,0.,100.,100,0.,100.);
+ TH2D *recoVsGenChargedJetENoDRHist = new TH2D("recoVsGenChargedJetENoDRHist","",100,0.,100.,100,0.,100.);
+ TH2D *recoVsGenChargedJetENoDupHist = new TH2D("recoVsGenChargedJetENoDup","",100,0.,100.,100,0.,100.);
+
+ TH2D *jetResVsEtaHist = new TH2D("jetResVsEta","",80,-4.,4.,10000,-10.,10.);
+ TH2D *jetResVsEHist = new TH2D("jetResVsE","",100,0.,100.,10000,-10.,10.);
+ TH2D *jetResVsENegEtaHist = new TH2D("jetResVsENegEta","",20,0.,100.,10000,-10.,10.);
+ TH2D *jetResVsEMidEtaHist = new TH2D("jetResVsEMidEta","",20,0.,100.,10000,-10.,10.);
+ TH2D *jetResVsEPosEtaHist = new TH2D("jetResVsEPosEta","",20,0.,100.,10000,-10.,10.);
- auto getJetTrackDRThresh = [](ints &type, ints &pid, floats &E, floats &x, floats &y, floats &z){
- floats dr;
- for(auto i=0U; i<x.size(); i++)
+ TH2D *jetResVsENegEtaNoDupHist = new TH2D("jetResVsENegEtaNoDup","",20,0.,100.,10000,-10.,10.);
+ TH2D *jetResVsEMidEtaNoDupHist = new TH2D("jetResVsEMidEtaNoDup","",20,0.,100.,10000,-10.,10.);
+ TH2D *jetResVsEPosEtaNoDupHist = new TH2D("jetResVsEPosEtaNoDup","",20,0.,100.,10000,-10.,10.);
+
+
+ // Loop Through Events
+ int NEVENTS = 0;
+ while(tree_reader.Next()) {
+
+ if(NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
+
+ counter->Fill(0);
+
+ //////////////////////////////////////////////////////////////////////////
+ ////////////////////// Analyze Reconstructed Jets //////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ int numRecoChargedJets = 0;
+ int numRecoChargedJetsNoElec = 0;
+ for(unsigned int i=0; i<recoType.GetSize(); i++)
{
- if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
+ TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
+
+ counter->Fill(3);
+
+ // Place eta cut to avoid edges of tracking acceptance
+ if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
+
+ // Place a minimum energy condition for several plots
+ bool ECut = recoNRG[i] > 5.0;
+
+ if(ECut) numRecoChargedJets++;
+
+ recoChargedJetEHist->Fill(recoNRG[i]);
+ if(ECut) recoChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
+ recoChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
+ if(ECut) recoChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
+
+ bool noElectron = true;
+ if(ECut)
{
- TVector3 jetMom(x[i],y[i],z[i]);
- int jetID = pid[i];
+ // Find Jets with Electrons
+ for(unsigned int m=partsBegin[i]; m<partsEnd[i]; m++)
+ {
+ if(recoPartM[recoPartIndex[m]] > 0.00050 && recoPartM[recoPartIndex[m]] < 0.00052)
+ noElectron = false;
+ }
- for(auto j=0U; j<x.size(); j++) // First loop over tracks
+ for(unsigned int j=partsBegin[i]; j<partsEnd[i]; j++)
{
- if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
+ // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
+ // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
+ double mX = recoPartMomX[recoPartIndex[j]];
+ double mY = recoPartMomY[recoPartIndex[j]];
+ double mZ = recoPartMomZ[recoPartIndex[j]];
+ double mM = recoPartM[recoPartIndex[j]];
+ double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
+
+ TVector3 partMom(mX,mY,mZ);
+
+ recoChargedJetPartEHist->Fill(tmpE);
+ recoChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
+ recoChargedJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
+ recoChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
+
+ if(noElectron)
{
- TVector3 trackA(x[j],y[j],z[j]);
+ recoChargedJetPartENoElecHist->Fill(tmpE);
+ recoChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
+ recoChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
+ recoChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
+ }
- if(j < x.size() - 1) // Overflow protection
+ // Pairwise Distance Between Constituents
+ if(j<(partsEnd[i]-1))
+ {
+ for(unsigned int k=j+1; k<partsEnd[i]; k++)
{
- 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
+ double mXB = recoPartMomX[recoPartIndex[k]];
+ double mYB = recoPartMomY[recoPartIndex[k]];
+ double mZB = recoPartMomZ[recoPartIndex[k]];
+
+ TVector3 partMomB(mXB,mYB,mZB);
+
+ double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ recoChargedJetPartPairwiseDeltaRHist->Fill(dR);
+ }
+ }
+ }
+ numRecoChargedJetPartsHist->Fill(partsEnd[i] - partsBegin[i]);
+ if(noElectron) numRecoChargedJetPartsNoElecHist->Fill(partsEnd[i] - partsBegin[i]);
+ }
+
+ // No Electrons
+ if(noElectron)
{
- TVector3 jetMom(x[i],y[i],z[i]);
- int jetID = pid[i];
+ recoChargedJetENoElecHist->Fill(recoNRG[i]);
+ if(ECut) recoChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
+ recoChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
+ if(ECut) recoChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
- 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(ECut) numRecoChargedJetsNoElec++;
+ }
+ }
+ numRecoChargedJetsECutHist->Fill(numRecoChargedJets);
+ numRecoChargedJetsECutNoElecHist->Fill(numRecoChargedJetsNoElec);
- 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++)
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////// Analyze Generator Jets ////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ int numGenChargedJets = 0;
+ int numGenChargedJetsNoElec = 0;
+ for(unsigned int i=0; i<genType.GetSize(); i++)
{
- if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
- {
- int hasDupTrack = 0;
+ TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
+
+ counter->Fill(4);
+
+ // Place eta cut to avoid edges of tracking acceptance
+ if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
+
+ // Place a minimum energy condition for several plots
+ bool ECut = genNRG[i] > 5.0;
+
+ if(ECut) numGenChargedJets++;
+
+ genChargedJetEHist->Fill(genNRG[i]);
+ if(ECut) genChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
+ genChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
+ if(ECut) genChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
- TVector3 jetMom(x[i],y[i],z[i]);
- int jetID = pid[i];
+ bool noElectron = true;
+ if(ECut)
+ {
+ // Find Jets with Electrons
+ for(unsigned int m=genPartsBegin[i]; m<genPartsEnd[i]; m++)
+ {
+ if(mcM[genPartIndex[m]] > 0.00050 && mcM[genPartIndex[m]] < 0.00052)
+ noElectron = false;
+ }
- for(auto j=0U; j<x.size(); j++) // First loop over tracks
+ for(unsigned int j=genPartsBegin[i]; j<genPartsEnd[i]; j++)
{
- if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
+ // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
+ // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
+ double mX = mcMomX[genPartIndex[j]];
+ double mY = mcMomY[genPartIndex[j]];
+ double mZ = mcMomZ[genPartIndex[j]];
+ double mM = mcM[genPartIndex[j]];
+ double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
+
+ TVector3 partMom(mX,mY,mZ);
+
+ genChargedJetPartEHist->Fill(tmpE);
+ genChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
+ genChargedJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
+ genChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
+
+ if(noElectron)
{
- TVector3 trackA(x[j],y[j],z[j]);
+ genChargedJetPartENoElecHist->Fill(tmpE);
+ genChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
+ genChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
+ genChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
+ }
- if(j < x.size() - 1) // Overflow protection
+ // Pairwise Distance Between Constituents
+ if(j<(genPartsEnd[i]-1))
+ {
+ for(unsigned int k=j+1; k<genPartsEnd[i]; k++)
{
- 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
+ double mXB = mcMomX[genPartIndex[k]];
+ double mYB = mcMomY[genPartIndex[k]];
+ double mZB = mcMomZ[genPartIndex[k]];
- auto getEta = [](vecs &x){
- floats eta;
- for(auto i=0U; i<x.size(); i++)
- {
- eta.push_back(x[i].PseudoRapidity());
- }
- return eta;
- }; // Jet Psudorapidity
+ TVector3 partMomB(mXB,mYB,mZB);
- auto getPhi = [](vecs &x){
- floats phi;
- for(auto i=0U; i<x.size(); i++)
- {
- phi.push_back(x[i].Phi());
+ double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ genChargedJetPartPairwiseDeltaRHist->Fill(dR);
+ }
+ }
+ }
+ numGenChargedJetPartsHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
+ if(noElectron) numGenChargedJetPartsNoElecHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
+ }
+
+ // No Electrons
+ if(noElectron)
+ {
+ genChargedJetENoElecHist->Fill(genNRG[i]);
+ if(ECut) genChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
+ genChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
+ if(ECut) genChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
+
+ if(ECut) numGenChargedJetsNoElec++;
+ }
}
- return phi;
- }; // Jet Phi
+ numGenChargedJetsECutHist->Fill(numGenChargedJets);
+ numGenChargedJetsECutNoElecHist->Fill(numGenChargedJetsNoElec);
- auto getMatchDeltaR = [](vecs &vecR, vecs &vecG){
- floats deltaR;
- for(auto i=0U; i<vecR.size(); i++)
+
+ //////////////////////////////////////////////////////////////////////////
+ ///////////////////////////// Matched Jets /////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ for(unsigned int i=0; i<genType.GetSize(); i++)
{
- float minDeltaR = 10000.0;
- for(auto j=0U; j<vecG.size(); j++)
+ TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
+
+ // Place eta cut to avoid edges of tracking acceptance
+ //if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
+
+ // Place a minimum energy condition
+ //if(genNRG[i] < 5.0) continue;
+
+ // Don't Look at Electron Jets
+ bool hasElectron = false;
+ // Find Jets with Electrons
+ for(unsigned int m=genPartsBegin[i]; m<genPartsEnd[i]; m++)
{
- 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++)
+ if(mcM[genPartIndex[m]] > 0.00050 && mcM[genPartIndex[m]] < 0.00052)
+ hasElectron = true;
+ }
+ //if(hasElectron) continue;
+
+ // Find Matching Reconstructed Jet
+ double minDeltaR = 999.;
+ int minIndex = -1;
+ for(unsigned int j=0; j<recoType.GetSize(); j++)
{
- double dEta = vecR[i].PseudoRapidity() - vecG[j].PseudoRapidity();
- double dPhi = TVector2::Phi_mpi_pi(vecR[i].Phi() - vecG[j].Phi());
+ TVector3 recoMom(recoMomX[j],recoMomY[j],recoMomZ[j]);
+
+ double dEta = jetMom.PseudoRapidity() - recoMom.PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - recoMom.Phi());
double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
-
- if(dR < minDeltaR)
+
+ if(dR < minDeltaR)
{
minDeltaR = dR;
- minE = E[j];
+ minIndex = 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++)
+ }
+
+ // Do Backwards Match
+ double minDeltaRBack = 999.;
+ double minIndexBack = -1;
+ if(minIndex > -1)
{
- 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)
+ TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
+ for(unsigned int j=0; j<genType.GetSize(); j++)
{
- minDeltaR = dR;
- minEta = vecG[j].PseudoRapidity();
+ TVector3 genMom(genMomX[j],genMomY[j],genMomZ[j]);
+
+ double dEta = recoMatchMom.PseudoRapidity() - genMom.PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(recoMatchMom.Phi() - genMom.Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < minDeltaRBack)
+ {
+ minDeltaRBack = dR;
+ minIndexBack = j;
+ }
}
- }
- 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++)
+ // Look at Best Match
+ if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1 && !hasElectron) matchJetDeltaRHist->Fill(minDeltaR);
+ if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1) matchJetDeltaRBackHist->Fill(minDeltaR);
+ if(minIndex > -1 && genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && !hasElectron)
{
- 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)
+ TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
+
+ recoVsGenChargedJetENoDRHist->Fill(genNRG[i],recoNRG[minIndex]);
+
+ if(minDeltaR < DELTARCUT)
{
- minDeltaR = dR;
- minPhi = vecG[j].Phi();
+ recoVsGenChargedJetEtaHist->Fill(jetMom.PseudoRapidity(),recoMatchMom.PseudoRapidity());
+ recoVsGenChargedJetPhiHist->Fill(jetMom.Phi(),recoMatchMom.Phi());
+ recoVsGenChargedJetEHist->Fill(genNRG[i],recoNRG[minIndex]);
+
+ double jetERes = (recoNRG[minIndex] - genNRG[i])/genNRG[i];
+
+ jetResVsEtaHist->Fill(jetMom.PseudoRapidity(),jetERes);
+ jetResVsEHist->Fill(genNRG[i],jetERes);
+ if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
+ jetResVsENegEtaHist->Fill(genNRG[i],jetERes);
+ if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
+ jetResVsEMidEtaHist->Fill(genNRG[i],jetERes);
+ if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
+ jetResVsEPosEtaHist->Fill(genNRG[i],jetERes);
+
+ // Check for Duplicate Tracks
+ bool noDuplicate = true;
+ for(unsigned int j=partsBegin[minIndex]; j<partsEnd[minIndex]; j++)
+ {
+ double mX = recoPartMomX[recoPartIndex[j]];
+ double mY = recoPartMomY[recoPartIndex[j]];
+ double mZ = recoPartMomZ[recoPartIndex[j]];
+ double mM = recoPartM[recoPartIndex[j]];
+ double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
+
+ TVector3 partMom(mX,mY,mZ);
+
+ // Pairwise Distance Between Constituents
+ if(j<(partsEnd[minIndex]-1))
+ {
+ for(unsigned int k=j+1; k<partsEnd[minIndex]; k++)
+ {
+ double mXB = recoPartMomX[recoPartIndex[k]];
+ double mYB = recoPartMomY[recoPartIndex[k]];
+ double mZB = recoPartMomZ[recoPartIndex[k]];
+
+ TVector3 partMomB(mXB,mYB,mZB);
+
+ double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
+ double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
+ double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
+
+ if(dR < 0.02) noDuplicate = false;
+ }
+ }
+ }
+
+ if(noDuplicate)
+ {
+ recoVsGenChargedJetENoDupHist->Fill(genNRG[i],recoNRG[minIndex]);
+
+ if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
+ jetResVsENegEtaNoDupHist->Fill(genNRG[i],jetERes);
+ if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
+ jetResVsEMidEtaNoDupHist->Fill(genNRG[i],jetERes);
+ if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
+ jetResVsEPosEtaNoDupHist->Fill(genNRG[i],jetERes);
+ }
}
- }
- 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;
- };
+ NEVENTS++;
+ }
- 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;
- };
+
+ gStyle->SetOptStat(0);
+ //////////////////////// Reconstructed Jets Plots ////////////////////////
+ // Reco Number
+ TCanvas *c1 = new TCanvas("c1","Number Reco Jets",800,600);
+ c1->Clear();
+ c1->Divide(1,1);
- 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;
- };
+ c1->cd(1);
+ numRecoChargedJetsECutHist->Draw("HIST");
+ numRecoChargedJetsECutNoElecHist->SetLineColor(seabornRed);
+ numRecoChargedJetsECutNoElecHist->Draw("HISTSAME");
+ numRecoChargedJetsECutHist->SetLineWidth(2); // Set line width to 2 (adjust as needed)
+numRecoChargedJetsECutNoElecHist->SetLineWidth(2);
+ numRecoChargedJetsECutHist->SetTitle("Reconstructed Jets per Event (|eta| < 2.5 && E > 5);Number");
- 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();
+TLegend *legend1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+legend1->AddEntry(numRecoChargedJetsECutHist, "With Electrons", "l");
+legend1->AddEntry(numRecoChargedJetsECutNoElecHist, "No Electrons", "l");
+legend1->Draw();
+
+
+
+ gPad->SetLogy();
+ if(PRINT) c1->Print((results_path+"/jets/numberRecoJets.png").c_str()); // Number of reconstructed jets per event with energy > 5 GeV and Abs(eta) < 2.5
+ delete c1;
+
+ // Reco Energy
+ TCanvas *c2 = new TCanvas("c2","Reco Jet Energy",800,600);
+ c2->Clear();
+ c2->Divide(1,1);
+
+ c2->cd(1);
+ recoChargedJetEHist->Draw("HIST");
+ recoChargedJetENoElecHist->SetLineColor(seabornRed);
+ recoChargedJetENoElecHist->Draw("HISTSAME");
+
+ recoChargedJetEHist->SetLineWidth(2);
+ recoChargedJetENoElecHist->SetLineWidth(2);
+ recoChargedJetEHist->SetTitle("Reconstructed Jet Energy (|eta| < 2.5);Energy [GeV]");
+
+TLegend *legend2 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+legend2->AddEntry(recoChargedJetEHist, "With Electrons", "l");
+legend2->AddEntry(recoChargedJetENoElecHist, "No Electrons", "l");
+legend2->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c2->Print((results_path+"/jets/recoJetEnergy.png").c_str()); // Energy spectrum of reconstructed jets with Abs(eta) < 2.5
+
+ delete c2;
+ // Reco Eta
+ TCanvas *c3 = new TCanvas("c3","Reco Jet Eta",800,600);
+ c3->Clear();
+ c3->Divide(1,1);
+
+ c3->cd(1);
+ recoChargedJetEtaECutHist->Draw("HIST");
+ recoChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
+ recoChargedJetEtaECutNoElecHist->Draw("HISTSAME");
+
+ recoChargedJetEtaECutHist->SetLineWidth(2);
+ recoChargedJetEtaECutNoElecHist->SetLineWidth(2);
+ recoChargedJetEtaECutHist->SetTitle("Reconstructed Jet Eta (E > 5);Eta");
+
+//add legend
+TLegend *legend3 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+legend3->AddEntry(recoChargedJetEtaECutHist, "With Electrons", "l");
+legend3->AddEntry(recoChargedJetEtaECutNoElecHist, "No Electrons", "l");
+legend3->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c3->Print((results_path+"/jets/recoJetEta.png").c_str()); // Eta spectrum of reconstructed jets with energy > 5 GeV
+ delete c3;
+ // Reco E Vs Eta
+ TCanvas *c4 = new TCanvas("c4","Reco Jet E Vs Eta",800,600);
+ c4->Clear();
+ c4->Divide(1,1);
+
+ c4->cd(1);
+ recoChargedJetEvsEtaHist->Draw("COLZ");
+ recoChargedJetEvsEtaHist->SetTitle("Reconstructed Jet Energy Vs Eta;Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c4->Print((results_path+"/jets/recoJetEnergyvsEta.png").c_str()); // Energy vs eta of reconstructed jets
+
+ // Reco Phi Vs Eta
+ TCanvas *c5 = new TCanvas("c5","Reco Jet Phi Vs Eta",800,600);
+ c5->Clear();
+ c5->Divide(1,1);
+
+ c5->cd(1);
+ recoChargedJetPhiVsEtaECutHist->Draw("COLZ");
+ recoChargedJetPhiVsEtaECutHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5);Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c5->Print((results_path+"/jets/recoJetPhivsEta.png").c_str()); // Phi vs eta of reconstructed jets
+
+ // Num Particles Per Reco Jet
+ TCanvas *c6 = new TCanvas("c6","Number Constituents Per Reco Jet",800,600);
+ c6->Clear();
+ c6->Divide(1,1);
+
+ c6->cd(1);
+ numRecoChargedJetPartsHist->Draw("HIST");
+ numRecoChargedJetPartsNoElecHist->SetLineColor(seabornRed);
+ numRecoChargedJetPartsNoElecHist->Draw("HISTSAME");
+
+ numRecoChargedJetPartsHist->SetLineWidth(2);
+ numRecoChargedJetPartsNoElecHist->SetLineWidth(2);
+ numRecoChargedJetPartsHist->SetTitle("Number of Constituents Per Reco Jet;Number of Constituents");
+
+TLegend *legend6 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+legend6->AddEntry(numRecoChargedJetPartsHist, "With Electrons", "l");
+legend6->AddEntry(numRecoChargedJetPartsNoElecHist, "No Electrons", "l");
+legend6->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c6->Print((results_path+"/jets/numConstituentsPerRecoJet.png").c_str()); // Number of constituents in reconstructed jets
+
+ // Reco Part Energy
+ TCanvas *c7 = new TCanvas("c7","Reco Jet Constituent Energy",800,600);
+ c7->Clear();
+ c7->Divide(1,1);
+
+ c7->cd(1);
+ recoChargedJetPartEHist->Draw("HIST");
+ recoChargedJetPartENoElecHist->SetLineColor(seabornRed);
+ recoChargedJetPartENoElecHist->Draw("HISTSAME");
+
+ recoChargedJetPartEHist->SetLineWidth(2);
+ recoChargedJetPartENoElecHist->SetLineWidth(2);
+ recoChargedJetPartEHist->SetTitle("Reconstructed Jet Constituent Energy;Energy [GeV]");
+
+ TLegend *legend7 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend7->AddEntry(recoChargedJetPartEHist, "With Electrons", "l");
+ legend7->AddEntry(recoChargedJetPartENoElecHist, "No Electrons", "l");
+ legend7->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c7->Print((results_path+"/jets/recoJetConstituentEnergy.png").c_str()); // Energy of reconstructed jet constituents
+
+ // Reco Part Eta
+ TCanvas *c8 = new TCanvas("c8","Reco Jet Constituent Eta",800,600);
+ c8->Clear();
+ c8->Divide(1,1);
+
+ c8->cd(1);
+ recoChargedJetPartEtaHist->Draw("HIST");
+ recoChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
+ recoChargedJetPartEtaNoElecHist->Draw("HISTSAME");
+
+ recoChargedJetPartEtaHist->SetLineWidth(2);
+ recoChargedJetPartEtaNoElecHist->SetLineWidth(2);
+
+ recoChargedJetPartEtaHist->SetTitle("Reconstructed Jet Constituent Eta;Eta");
+
+ TLegend *legend8 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend8->AddEntry(recoChargedJetPartEtaHist, "With Electrons", "l");
+ legend8->AddEntry(recoChargedJetPartEtaNoElecHist, "No Electrons", "l");
+ legend8->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c8->Print((results_path+"/jets/recoJetConstituentEta.png").c_str()); // Eta of reconstructed jet constituents
+
+ // Reco Part E Vs Eta
+ TCanvas *c9 = new TCanvas("c9","Reco Jet Constituent E Vs Eta",800,600);
+ c9->Clear();
+ c9->Divide(1,1);
+
+ c9->cd(1);
+ recoChargedJetPartEvsEtaHist->Draw("COLZ");
+ recoChargedJetPartEvsEtaHist->SetTitle("Reconstructed Jet Constituent Energy Vs Eta;Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c9->Print((results_path+"/jets/recoJetConstituentEnergyVsEta.png").c_str()); // Energy vs eta of reconstructed jet constituents
+
+ // Reco Part Phi Vs Eta
+ TCanvas *c10 = new TCanvas("c10","Reco Jet Constituent Phi Vs Eta",800,600);
+ c10->Clear();
+ c10->Divide(1,1);
+
+ c10->cd(1);
+ recoChargedJetPartPhiVsEtaHist->Draw("COLZ");
+ recoChargedJetPartPhiVsEtaHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta;Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c10->Print((results_path+"/jets/recoJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of reconstructed jet constituents
+
+ // Reco Constituent Pairwise delta R
+ TCanvas *c11 = new TCanvas("c11","Reco Jet Constituent Pairwise Delta R",800,600);
+ c11->Clear();
+ c11->Divide(1,1);
+
+ c11->cd(1);
+ recoChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
+ recoChargedJetPartPairwiseDeltaRHist->SetTitle("Pairwise Constituent Delta R;Delta R");
+ recoChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
+ gPad->SetLogy();
+ if(PRINT) c11->Print((results_path+"/jets/recoJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in reconstructed jets
+
+ // Reco E Vs Eta No Electron Jets
+ TCanvas *c12 = new TCanvas("c12","Reco Jet E Vs Eta (No Electrons)",800,600);
+ c12->Clear();
+ c12->Divide(1,1);
+
+ c12->cd(1);
+ recoChargedJetEvsEtaNoElecHist->Draw("COLZ");
+ recoChargedJetEvsEtaNoElecHist->SetTitle("Reconstructed Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c12->Print((results_path+"/jets/recoJetEnergyVsEtaNoElectron.png").c_str()); // Reconstructed jet energy - no jets containing electrons included
+
+ // Reco Phi Vs Eta No Electron Jets
+ TCanvas *c13 = new TCanvas("c13","Reco Jet Phi Vs Eta (No Electrons)",800,600);
+ c13->Clear();
+ c13->Divide(1,1);
+
+ c13->cd(1);
+ recoChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
+ recoChargedJetPhiVsEtaECutNoElecHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c13->Print((results_path+"/jets/recoJetPhiVsEtaNoElectron.png").c_str()); // Reconstructed Jet phi vs eta - no jets containing electrons included
+
+ // Reco Part E Vs Eta No Electron Jets
+ TCanvas *c14 = new TCanvas("c14","Reco Jet Constituent E Vs Eta (No Electrons)",800,600);
+ c14->Clear();
+ c14->Divide(1,1);
+
+ c14->cd(1);
+ recoChargedJetPartEvsEtaNoElecHist->Draw("COLZ");
+ recoChargedJetPartEvsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c14->Print((results_path+"/jets/recoJetConstituentEnergyVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent energy vs eta - no jets containing electrons included
+
+ // Reco Part Phi Vs Eta No Electron Jets
+ TCanvas *c15 = new TCanvas("c15","Reco Jet Constituent Phi Vs Eta (No Electrons)",800,600);
+ c15->Clear();
+ c15->Divide(1,1);
+
+ c15->cd(1);
+ recoChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
+ recoChargedJetPartPhiVsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c15->Print((results_path+"/jets/recoJetConstituentPhiVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent phi vs eta - no jets containing electrons included
+
+
+ //////////////////////// Generated Jets Plots ////////////////////////
+ // Gen Number
+ TCanvas *c16 = new TCanvas("c16","Number Gen Jets",800,600);
+ c16->Clear();
+ c16->Divide(1,1);
+
+ c16->cd(1);
+ numGenChargedJetsECutHist->Draw("HIST");
+ numGenChargedJetsECutNoElecHist->SetLineColor(seabornRed);
+ numGenChargedJetsECutNoElecHist->Draw("HISTSAME");
+
+ numGenChargedJetsECutHist->SetLineWidth(2);
+ numGenChargedJetsECutNoElecHist->SetLineWidth(2);
+
+ numGenChargedJetsECutHist->SetTitle("Generator Jets per Event (|eta| < 2.5 && E > 5);Number");
+
+ TLegend *legend16 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend16->AddEntry(numGenChargedJetsECutHist, "With Electrons", "l");
+ legend16->AddEntry(numGenChargedJetsECutNoElecHist, "No Electrons", "l");
+ legend16->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c16->Print((results_path+"/jets/numberGenJets.png").c_str()); // Number of generator jets per event with energy > 5 GeV and Abs(eta) < 2.5
+
+ // Gen Energy
+ TCanvas *c17 = new TCanvas("c17","Gen Jet Energy",800,600);
+ c17->Clear();
+ c17->Divide(1,1);
+
+ c17->cd(1);
+ genChargedJetEHist->Draw("HIST");
+ genChargedJetENoElecHist->SetLineColor(seabornRed);
+ genChargedJetENoElecHist->Draw("HISTSAME");
+
+ genChargedJetEHist->SetLineWidth(2);
+ genChargedJetENoElecHist->SetLineWidth(2);
+
+ genChargedJetEHist->SetTitle("Generator Jet Energy (|eta| < 2.5);Energy [GeV]");
+ TLegend *legend17 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend17->AddEntry(genChargedJetEHist, "With Electrons", "l");
+ legend17->AddEntry(genChargedJetENoElecHist, "No Electrons", "l");
+ legend17->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c17->Print((results_path+"/jets/genJetEnergy.png").c_str()); // Energy spectrum of generated jets with Abs(eta) < 2.5
+
+ // Gen Eta
+ TCanvas *c18 = new TCanvas("c18","Gen Jet Eta",800,600);
+ c18->Clear();
+ c18->Divide(1,1);
+
+ c18->cd(1);
+ genChargedJetEtaECutHist->Draw("HIST");
+ genChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
+ genChargedJetEtaECutNoElecHist->Draw("HISTSAME");
+
+ genChargedJetEtaECutHist->SetLineWidth(2);
+ genChargedJetEtaECutNoElecHist->SetLineWidth(2);
+
+ genChargedJetEtaECutHist->SetTitle("Generator Jet Eta (E > 5);Eta");
+
+ TLegend *legend18 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend18->AddEntry(genChargedJetEtaECutHist, "With Electrons", "l");
+ legend18->AddEntry(genChargedJetEtaECutNoElecHist, "No Electrons", "l");
+ legend18->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c18->Print((results_path+"/jets/genJetEta.png").c_str()); // Eta spectrum of generator jets with energy > 5 GeV
+
+ // Gen E Vs Eta
+ TCanvas *c19 = new TCanvas("c19","Gen Jet E Vs Eta",800,600);
+ c19->Clear();
+ c19->Divide(1,1);
+
+ c19->cd(1);
+ genChargedJetEvsEtaHist->Draw("COLZ");
+ genChargedJetEvsEtaHist->SetTitle("Generator Jet Energy Vs Eta;Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c19->Print((results_path+"/jets/genJetEnergyvsEta.png").c_str()); // Energy vs eta of generator jets
+
+ // Gen Phi Vs Eta
+ TCanvas *c20 = new TCanvas("c20","Gen Jet Phi Vs Eta",800,600);
+ c20->Clear();
+ c20->Divide(1,1);
+
+ c20->cd(1);
+ genChargedJetPhiVsEtaECutHist->Draw("COLZ");
+ genChargedJetPhiVsEtaECutHist->SetTitle("Generator Jet Phi Vs Eta (E > 5);Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c20->Print((results_path+"/jets/genJetPhiVsEta.png").c_str()); // Phi vs eta of generator jets
+
+ // Num Particles Per Gen Jet
+ TCanvas *c21 = new TCanvas("c21","Number Constituents Per Gen Jet",800,600);
+ c21->Clear();
+ c21->Divide(1,1);
+
+ c21->cd(1);
+ numGenChargedJetPartsHist->Draw("HIST");
+ numGenChargedJetPartsNoElecHist->SetLineColor(seabornRed);
+ numGenChargedJetPartsNoElecHist->Draw("HISTSAME");
+
+ numGenChargedJetPartsHist->SetLineWidth(2);
+ numGenChargedJetPartsNoElecHist->SetLineWidth(2);
+
+ numGenChargedJetPartsHist->SetTitle("Number of Constituents Per Gen Jet;Number of Constituents");
+
+ TLegend *legend21 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend21->AddEntry(numGenChargedJetPartsHist, "With Electrons", "l");
+ legend21->AddEntry(numGenChargedJetPartsNoElecHist, "No Electrons", "l");
+ legend21->Draw();
+ gPad->SetLogy();
+ if(PRINT) c21->Print((results_path+"/jets/numConstituentsPerGenJet.png").c_str()); // Number of constituents in generator jets
+
+ // Gen Part Energy
+ TCanvas *c22 = new TCanvas("c22","Gen Jet Constituent Energy",800,600);
+ c22->Clear();
+ c22->Divide(1,1);
+
+ c22->cd(1);
+ genChargedJetPartEHist->Draw("HIST");
+ genChargedJetPartENoElecHist->SetLineColor(seabornRed);
+ genChargedJetPartENoElecHist->Draw("HISTSAME");
+
+ genChargedJetPartEHist->SetLineWidth(2);
+ genChargedJetPartENoElecHist->SetLineWidth(2);
+
+ genChargedJetPartEHist->SetTitle("Generator Jet Constituent Energy;Energy [GeV]");
+
+ TLegend *legend22 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend22->AddEntry(genChargedJetPartEHist, "With Electrons", "l");
+ legend22->AddEntry(genChargedJetPartENoElecHist, "No Electrons", "l");
+ legend22->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c22->Print((results_path+"/jets/genJetConstituentEnergy.png").c_str()); // Energy of generator jet constituents
+
+ // Gen Part Eta
+ TCanvas *c23 = new TCanvas("c23","Gen Jet Constituent Eta",800,600);
+ c23->Clear();
+ c23->Divide(1,1);
+
+ c23->cd(1);
+ genChargedJetPartEtaHist->Draw("HIST");
+ genChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
+ genChargedJetPartEtaNoElecHist->Draw("HISTSAME");
+
+ genChargedJetPartEtaHist->SetLineWidth(2);
+ genChargedJetPartEtaNoElecHist->SetLineWidth(2);
+
+ genChargedJetPartEtaHist->SetTitle("Generator Jet Constituent Eta;Eta");
+
+ TLegend *legend23 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
+ legend23->AddEntry(genChargedJetPartEtaHist, "With Electrons", "l");
+ legend23->AddEntry(genChargedJetPartEtaNoElecHist, "No Electrons", "l");
+ legend23->Draw();
+
+ gPad->SetLogy();
+ if(PRINT) c23->Print((results_path+"/jets/genJetConstituentEta.png").c_str()); // Eta of generator jet constituents
+
+ // Gen Part E Vs Eta
+ TCanvas *c24 = new TCanvas("c24","Gen Jet Constituent E Vs Eta",800,600);
+ c24->Clear();
+ c24->Divide(1,1);
+
+ c24->cd(1);
+ genChargedJetPartEvsEtaHist->Draw("COLZ");
+ genChargedJetPartEvsEtaHist->SetTitle("Generator Jet Constituent Energy Vs Eta;Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c24->Print((results_path+"/jets/genJetConstituentEnergyVsEta.png").c_str()); // Energy vs eta of generator jet constituents
+
+ // Gen Part Phi Vs Eta
+ TCanvas *c25 = new TCanvas("c25","Gen Jet Constituent Phi Vs Eta",800,600);
+ c25->Clear();
+ c25->Divide(1,1);
+
+ c25->cd(1);
+ genChargedJetPartPhiVsEtaHist->Draw("COLZ");
+ genChargedJetPartPhiVsEtaHist->SetTitle("Generator Jet Constituent Phi Vs Eta;Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c25->Print((results_path+"/jets/genJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of generator jet constituents
+
+ // Gen Constituent Pairwise delta R
+ TCanvas *c26 = new TCanvas("c26","Gen Jet Constituent Pairwise Delta R",800,600);
+ c26->Clear();
+ c26->Divide(1,1);
+
+ c26->cd(1);
+ genChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
+ genChargedJetPartPairwiseDeltaRHist->SetTitle("Generator Jet Pairwise Constituent Delta R;Delta R");
+ genChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
+ gPad->SetLogy();
+ if(PRINT) c26->Print((results_path+"/jets/genJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in generator jets
+
+ // Gen E Vs Eta No Electron Jets
+ TCanvas *c27 = new TCanvas("c27","Gen Jet E Vs Eta (No Electrons)",800,600);
+ c27->Clear();
+ c27->Divide(1,1);
+
+ c27->cd(1);
+ genChargedJetEvsEtaNoElecHist->Draw("COLZ");
+ genChargedJetEvsEtaNoElecHist->SetTitle("Generator Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c27->Print((results_path+"/jets/genJetEnergyVsEtaNoElectron.png").c_str()); // Generator jet energy vs eta - no jets containing electrons included
+
+ // Gen Phi Vs Eta No Electron Jets
+ TCanvas *c28 = new TCanvas("c28","Gen Jet Phi Vs Eta (No Electrons)",800,600);
+ c28->Clear();
+ c28->Divide(1,1);
+
+ c28->cd(1);
+ genChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
+ genChargedJetPhiVsEtaECutNoElecHist->SetTitle("Generator Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
+ gPad->SetLogz();
+ if(PRINT) c28->Print((results_path+"/jets/genJetPhiVsEtaNoElectron.png").c_str()); // Generator Jet phi vs eta - no jets containing electrons included
+
+ // Gen Part E Vs Eta No Electron Jets
+ TCanvas *c29 = new TCanvas("c29","Gen Jet Constituent E Vs Eta (No Electrons)",800,600);
+ c29->Clear();
+ c29->Divide(1,1);
+
+ c29->cd(1);
+ genChargedJetPartEvsEtaNoElecHist->Draw("COLZ");
+ genChargedJetPartEvsEtaNoElecHist->SetTitle("Generator Jet Constituent Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
+ gPad->SetLogz();
+ if(PRINT) c29->Print((results_path+"/jets/genJetConstituentEnergyVsEtaNoElectron.png").c_str()); // Generator jet constituent energy vs eta - no jets containing electrons included
+
+ // Gen Part Phi Vs Eta No Electron Jets
+ TCanvas *c30 = new TCanvas("c30","Gen Jet Constituent Phi Vs Eta (No Electrons)",800,600);
+ c30->Clear();
+ c30->Divide(1,1);
+
+ c30->cd(1);
+ genChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
+ genChargedJetPartPhiVsEtaNoElecHist->SetTitle("Generator Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
+ gPad->SetLogz();
+ //c30->Print((results_path+"/jets/recoJetEvsEta.png").c_str());
+ if(PRINT) c30->Print((results_path+"/jets/genJetConstituentPhiVsEtaNoElectron.png").c_str()); // Generator jet constituent phi vs eta - no jets containing electrons included
+
+
+ //////////////////////// Matched Jets Plots ////////////////////////
+ // Matched Delta R
+ TCanvas *c31 = new TCanvas("c31","Gen - Reco Delta R",800,600);
+ c31->Clear();
+ c31->Divide(1,1);
+
+ c31->cd(1);
+ matchJetDeltaRHist->Draw("HIST");
+ matchJetDeltaRBackHist->SetLineColor(seabornRed);
+ //matchJetDeltaRBackHist->Draw("HISTSAME");
+ matchJetDeltaRHist->SetTitle("Matched Gen - Reco Jet Delta R;Delta R");
+ gPad->SetLogy();
+ if(PRINT) c31->Print((results_path+"/jets/genRecoJetDeltaR.png").c_str()); // Distance between closest generated and reconstructed jet pair
+
+ // Matched Reco Vs Gen Eta
+ TCanvas *c32 = new TCanvas("c32","Reco Vs Gen Eta",800,600);
+ c32->Clear();
+ c32->Divide(1,1);
+
+ c32->cd(1);
+ recoVsGenChargedJetEtaHist->Draw("COLZ");
+ recoVsGenChargedJetEtaHist->SetTitle("Reconstructed Vs Generator Jet Eta;Gen Eta;Reco Eta");
+ gPad->SetLogz();
+ if(PRINT) c32->Print((results_path+"/jets/matchedRecoVsGenJetEta.png").c_str()); // Matched Reconstructed Vs Generator Jet eta
+
+ // Matched Reco Vs Gen Phi
+ TCanvas *c33 = new TCanvas("c33","Reco Vs Gen Phi",800,600);
+ c33->Clear();
+ c33->Divide(1,1);
+
+ c33->cd(1);
+ recoVsGenChargedJetPhiHist->Draw("COLZ");
+ recoVsGenChargedJetPhiHist->SetTitle("Reconstructed Vs Generator Jet Phi;Gen Phi;Reco Phi");
+ gPad->SetLogz();
+ if(PRINT) c33->Print((results_path+"/jets/matchedRecoVsGenJetPhi.png").c_str()); // Matched reconstructed vs generator jet phi
+
+ // Matched Reco Vs Gen Energy
+ TCanvas *c34 = new TCanvas("c34","Reco Vs Gen Energy",800,600);
+ c34->Clear();
+ c34->Divide(1,1);
+
+ TF1 *f1_34 = new TF1("f1_34","1.0*x + 0.0",1,100);
+ TF1 *f2_34 = new TF1("f2_34","2.0*x + 0.0",1,100);
+ TF1 *f3_34 = new TF1("f3_34","3.0*x + 0.0",1,100);
+
+ c34->cd(1);
+ recoVsGenChargedJetEHist->Draw("COLZ");
+ recoVsGenChargedJetEHist->SetTitle("Reconstructed Vs Generator Jet Energy;Gen E;Reco E");
+ f1_34->Draw("SAME");
+ f2_34->Draw("SAME");
+ f3_34->Draw("SAME");
+ gPad->SetLogz();
+ if(PRINT) c34->Print((results_path+"/jets/matchedRecoVsGenJetEnergy.png").c_str()); // Matched reconstructed vs generator jet energy
+
+ // Jet Res Vs Gen Eta
+ TCanvas *c35 = new TCanvas("c35","Jet Res Vs Gen Eta",800,600);
+ c35->Clear();
+ c35->Divide(1,1);
+
+ c35->cd(1);
+ jetResVsEtaHist->Draw("COLZ");
+ jetResVsEtaHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Eta;Gen Eta;Res");
+ gPad->SetLogz();
+ if(PRINT) c35->Print((results_path+"/jets/matchedJetResolutionVsEta.png").c_str()); // Matched jet resolution vs generator jet eta
+
+ // Jet Res Vs Gen E
+ TCanvas *c36 = new TCanvas("c36","Jet Res Vs Gen E",800,600);
+ c36->Clear();
+ c36->Divide(1,1);
+
+ c36->cd(1);
+ jetResVsEHist->Draw("COLZ");
+ jetResVsEHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy;Gen E;Res");
+ gPad->SetLogz();
+ if(PRINT) c36->Print((results_path+"/jets/matchedJetResolutionVsEnergy.png").c_str()); // Matched jet resolution vs generator jet energy
+
+ // Jet Res Vs Gen E Neg Eta
+ TCanvas *c37 = new TCanvas("c37","Jet Res Vs Gen E (-2.5 < eta < -1.0)",800,600);
+ c37->Clear();
+ c37->Divide(1,1);
+
+ c37->cd(1);
+ jetResVsENegEtaNoDupHist->Draw("COLZ");
+ jetResVsENegEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-2.5 < eta < -1.0) No Duplicate;Gen E;Res");
+ gPad->SetLogz();
+ if(PRINT) c37->Print((results_path+"/jets/matchedJetResolutionVsEnergyNegEta.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0
+
+ // Jet Res Vs Gen E Mid Eta
+ TCanvas *c38 = new TCanvas("c38","Jet Res Vs Gen E (-1.0 < eta < 1.0)",800,600);
+ c38->Clear();
+ c38->Divide(1,1);
+
+ c38->cd(1);
+ jetResVsEMidEtaNoDupHist->Draw("COLZ");
+ jetResVsEMidEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-1.0 < eta < 1.0) No Duplicate;Gen E;Res");
+ gPad->SetLogz();
+ if(PRINT) c38->Print((results_path+"/jets/matchedJetResolutionVsEnergyMidEta.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0
+ delete c38;
+ // Jet Res Vs Gen E Pos Eta
+ TCanvas *c39 = new TCanvas("c39","Jet Res Vs Gen E (1.0 < eta < 2.5)",800,600);
+ c39->Clear();
+ c39->Divide(1,1);
+
+ c39->cd(1);
+ jetResVsEPosEtaNoDupHist->Draw("COLZ");
+ jetResVsEPosEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (1.0 < eta < 2.5) No Duplicate;Gen E;Res");
+ gPad->SetLogz();
+ if(PRINT) c39->Print((results_path+"/jets/matchedJetResolutionVsEnergyPosEta.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5
+ delete c39;
+
+
+ // Generate Resolution Plots
+ const int BINS = 20;
+ double binCent[BINS];
+ double jesVsENeg[BINS];
+ double jesVsEMid[BINS];
+ double jesVsEPos[BINS];
+ double jerVsENeg[BINS];
+ double jerVsEMid[BINS];
+ double jerVsEPos[BINS];
+
+ std::fill(std::begin(binCent), std::end(binCent), -999.);
+ std::fill(std::begin(jesVsENeg), std::end(jesVsENeg), -999.);
+ std::fill(std::begin(jesVsEMid), std::end(jesVsEMid), -999.);
+ std::fill(std::begin(jesVsEPos), std::end(jesVsEPos), -999.);
+ std::fill(std::begin(jerVsENeg), std::end(jerVsENeg), -999.);
+ std::fill(std::begin(jerVsEMid), std::end(jerVsEMid), -999.);
+ std::fill(std::begin(jerVsEPos), std::end(jerVsEPos), -999.);
+
+ TH1D *pxA = jetResVsENegEtaNoDupHist->ProjectionX("pxA",1,10000);
+ for(int i=0; i<BINS; i++)
+ {
+ binCent[i] = pxA->GetBinCenter(i+1);
+ }
+
+ TCanvas *c40 = new TCanvas("c40","Negative Rapidity Fit Results",800,600);
+ c40->Clear();
+ c40->Divide(5,4);
+
+ TH1D *hA[20];
+ for(int i=1; i<21; i++)
+ {
+ hA[i-1] = (TH1D *)jetResVsENegEtaNoDupHist->ProjectionY(Form("projYA_%d",i),i,i);
+
+ TF1 *myGausA = new TF1("myGausA","gaus",-0.5,0.5);
+ myGausA->SetParameters(hA[i-1]->GetMaximum(),0.0,0.01);
+
+ c40->cd(i);
+ //hA[i-1]->Draw("HIST");
+ hA[i-1]->Fit("myGausA","B","",-0.5,0.5);
+ hA[i-1]->GetXaxis()->SetRangeUser(-1,1);
+ gPad->SetLogy();
+
+ if(hA[i-1]->GetEntries() > 2)
+ {
+ auto fA = hA[i-1]->GetFunction("myGausA");
+
+ if(fA->GetParError(2)/fA->GetParameter(2) < 0.5)
+ {
+ jesVsENeg[i-1] = fA->GetParameter(1);
+ jerVsENeg[i-1] = fA->GetParameter(2);
+ }
+ //cout << fA->GetParameter(0) << " " << fA->GetParameter(1) << " " << fA->GetParameter(2) << endl;
+ //cout << fA->GetParError(0) << " " << fA->GetParError(1) << " " << fA->GetParError(2) << endl;
+ }
+ }
+ if(PRINT) c40->Print((results_path+"/jets/matchedJetResolutionVsEnergyNegEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0 fits
+ delete c40;
+
+ TCanvas *c41 = new TCanvas("c41","Mid Rapidity Fit Results",800,600);
+ c41->Clear();
+ c41->Divide(5,4);
+
+ TH1D *hB[20];
+ for(int i=1; i<21; i++)
+ {
+ hB[i-1] = (TH1D *)jetResVsEMidEtaNoDupHist->ProjectionY(Form("projYB_%d",i),i,i);
+
+ TF1 *myGausB = new TF1("myGausB","gaus",-0.5,0.5);
+ myGausB->SetParameters(hB[i-1]->GetMaximum(),0.0,0.01);
+
+ c41->cd(i);
+ //hA[i-1]->Draw("HIST");
+ hB[i-1]->Fit("myGausB","B","",-0.5,0.5);
+ hB[i-1]->GetXaxis()->SetRangeUser(-1,1);
+ gPad->SetLogy();
+
+ if(hB[i-1]->GetEntries() > 2)
+ {
+ auto fB = hB[i-1]->GetFunction("myGausB");
+
+ if(fB->GetParError(2)/fB->GetParameter(2) < 0.5)
+ {
+ jesVsEMid[i-1] = fB->GetParameter(1);
+ jerVsEMid[i-1] = fB->GetParameter(2);
+ }
+ //cout << fB->GetParameter(0) << " " << fB->GetParameter(1) << " " << fB->GetParameter(2) << endl;
+ //cout << fB->GetParError(0) << " " << fB->GetParError(1) << " " << fB->GetParError(2) << endl;
+ }
+ }
+ if(PRINT) c41->Print((results_path+"/jets/matchedJetResolutionVsEnergyMidEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0 fits
+ delete c41;
+
+ TCanvas *c42 = new TCanvas("c42","Positive Rapidity Fit Results",800,600);
+ c42->Clear();
+ c42->Divide(5,4);
+
+ TH1D *hC[20];
+ for(int i=1; i<21; i++)
+ {
+ hC[i-1] = (TH1D *)jetResVsEPosEtaNoDupHist->ProjectionY(Form("projYC_%d",i),i,i);
+
+ TF1 *myGausC = new TF1("myGausC","gaus",-0.5,0.5);
+ myGausC->SetParameters(hC[i-1]->GetMaximum(),0.0,0.01);
+
+ c42->cd(i);
+ //hA[i-1]->Draw("HIST");
+ hC[i-1]->Fit("myGausC","B","",-0.5,0.5);
+ hC[i-1]->GetXaxis()->SetRangeUser(-1,1);
+ gPad->SetLogy();
+
+ if(hC[i-1]->GetEntries() > 2)
+ {
+ auto fC = hC[i-1]->GetFunction("myGausC");
+
+ if(fC->GetParError(2)/fC->GetParameter(2) < 0.5)
+ {
+ jesVsEPos[i-1] = fC->GetParameter(1);
+ jerVsEPos[i-1] = fC->GetParameter(2);
+ }
+ //cout << fC->GetParameter(0) << " " << fC->GetParameter(1) << " " << fC->GetParameter(2) << endl;
+ //cout << fC->GetParError(0) << " " << fC->GetParError(1) << " " << fC->GetParError(2) << endl;
+ }
+ }
+ if(PRINT) c42->Print((results_path+"/jets/matchedJetResolutionVsEnergyPosEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5 fits
+ delete c42;
+ TCanvas *c43 = new TCanvas("c43","Positive JES/JER",800,600);
+ c43->Clear();
+ c43->Divide(1,1);
+
+ TGraph *gJESvsENeg = new TGraph(BINS,binCent,jesVsENeg);
+ TGraph *gJERvsENeg = new TGraph(BINS,binCent,jerVsENeg);
+
+ TGraph *gJESvsEMid = new TGraph(BINS,binCent,jesVsEMid);
+ TGraph *gJERvsEMid = new TGraph(BINS,binCent,jerVsEMid);
+
+ TGraph *gJESvsEPos = new TGraph(BINS,binCent,jesVsEPos);
+ TGraph *gJERvsEPos = new TGraph(BINS,binCent,jerVsEPos);
+
+ TH2D *test43 = new TH2D("test43","Jet Energy Scale / Resolution Vs Eta;True Eta;JES/JER",1,0.,100.,1,-0.2,0.2);
+ test43->Draw();
+
+ c43->cd(1);
+ gJERvsENeg->Draw("*");
+ gJERvsENeg->SetMarkerStyle(21);
+ gJERvsENeg->SetMarkerSize(1);
+ gJERvsENeg->SetMarkerColor(seabornBlue);
+
+ gJESvsENeg->Draw("*");
+ gJESvsENeg->SetMarkerStyle(26);
+ gJESvsENeg->SetMarkerSize(1);
+ gJESvsENeg->SetMarkerColor(seabornBlue);
+
+ gJERvsEMid->Draw("*");
+ gJERvsEMid->SetMarkerStyle(21);
+ gJERvsEMid->SetMarkerSize(1);
+ gJERvsEMid->SetMarkerColor(seabornRed);
+
+ gJESvsEMid->Draw("*");
+ gJESvsEMid->SetMarkerStyle(26);
+ gJESvsEMid->SetMarkerSize(1);
+ gJESvsEMid->SetMarkerColor(seabornRed);
+
+ gJERvsEPos->Draw("*");
+ gJERvsEPos->SetMarkerStyle(21);
+ gJERvsEPos->SetMarkerSize(1);
+ gJERvsEPos->SetMarkerColor(seabornGreen);
+
+ gJESvsEPos->Draw("*");
+ gJESvsEPos->SetMarkerStyle(26);
+ gJESvsEPos->SetMarkerSize(1);
+ gJESvsEPos->SetMarkerColor(seabornGreen);
+
+ TLegend *legend = new TLegend(0.7,0.7,0.9,0.9);
+legend->AddEntry(gJERvsENeg, "JER, (-2.5 < #eta < -1)", "p");
+legend->AddEntry(gJESvsENeg, "JES, (-2.5 < #eta < -1)","p");
+legend->AddEntry(gJERvsEMid, "JER, (-1 < #eta < 1)", "p");
+legend->AddEntry(gJESvsEMid, "JES, (-1 < #eta < 1)", "p");
+legend->AddEntry(gJERvsEPos, "JER, (1 < #eta < 2.5) ", "p");
+legend->AddEntry(gJESvsEPos, "JES, (1 < #eta < 2.5)", "p");
+legend->Draw();
+
+ if(PRINT) c43->Print((results_path+"/jets/matchedJetScaleResolutionSummary.png").c_str()); // Matched jet JER/JES summary
+ delete c43;
+
+delete mychain;
return 0;
}
diff --git a/benchmarks/dis/config.yml b/benchmarks/dis/config.yml
index bb68fa08d39260872d6446e8bb21fd6fef6cd7dc..7afb5b45e6034e2fe0f8a52d700c15d8347b40fd 100644
--- a/benchmarks/dis/config.yml
+++ b/benchmarks/dis/config.yml
@@ -21,7 +21,12 @@ dis:simulate:
MINQ2: [1, 10, 100, 1000]
timeout: 2 hours
script:
- - snakemake --cores 1 results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/dis_${EBEAM}x${PBEAM}_minQ2=${MINQ2}dis_electrons.json results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/kinematics_correlations results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/truth_reconstruction
+ - |
+ snakemake --cores 1 \
+ results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/dis_${EBEAM}x${PBEAM}_minQ2=${MINQ2}dis_electrons.json \
+ results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/jets/ \
+ results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/kinematics_correlations/ \
+ results/epic_craterlake/dis/${EBEAM}on${PBEAM}/minQ2=${MINQ2}/truth_reconstruction/
retry:
max: 2
when: