Skip to content
Snippets Groups Projects
Unverified Commit 0d19f33d authored by Dmitry Kalinkin's avatar Dmitry Kalinkin Committed by GitHub
Browse files

dis/analysis/jets.cxx: update from the latest code (#13) 

* dis/analysis/jets.cxx: update from the latest code

* dis/config.yml: add back jets bechmarks

* update benchmark

* jets.cxx: add missing includes

* jets.cxx: fix string concatenation

* place plots to jets/

* update instruction
parent 97e60162
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
benchmarks/dis/analysis/jets.cxx
+ 1223
− 467
View file @ 0d19f33d
......@@ -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;
// Input
TChain *mychain = new TChain("events");
mychain->Add(rec_file.c_str());
const int seabornRed = TColor::GetColor(213, 94, 0);
// Seaborn Green: #009E73 -> (0, 158, 115)
const int seabornGreen = TColor::GetColor(0, 158, 115);
// Seaborn Blue: #56B4E9 -> (86, 180, 233)
const int seabornBlue = TColor::GetColor(100, 149, 237);
// Output
//TFile *ofile = TFile::Open("test_24-05-0.hist.root","RECREATE");
// TTreeReader
TTreeReader tree_reader(mychain);
// 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.);
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++)
{
TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
ROOT::EnableImplicitMT(kNumThreads);
ROOT::RDataFrame d("events", rec_file);
counter->Fill(3);
TFile output(fmt::format("{}/jets/{}.hist.root", results_path, plot_tag).c_str(), "RECREATE");
// Place eta cut to avoid edges of tracking acceptance
if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
// 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
// Place a minimum energy condition for several plots
bool ECut = recoNRG[i] > 5.0;
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;
};
if(ECut) numRecoChargedJets++;
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
recoChargedJetEHist->Fill(recoNRG[i]);
if(ECut) recoChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
recoChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
if(ECut) recoChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
auto getJetE = [](ints &type, floats &E){
floats nrg;
for(auto i=0U; i<E.size(); i++)
bool noElectron = true;
if(ECut)
{
// Find Jets with Electrons
for(unsigned int m=partsBegin[i]; m<partsEnd[i]; m++)
{
if(type[i] == 0) nrg.push_back(E[i]);
if(recoPartM[recoPartIndex[m]] > 0.00050 && recoPartM[recoPartIndex[m]] < 0.00052)
noElectron = false;
}
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++)
for(unsigned int j=partsBegin[i]; j<partsEnd[i]; j++)
{
TVector3 tmp(x[i],y[i],z[i]);
if(type[i] == 0 && E[i] >= ETHRESH) v.push_back(tmp);
// 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)
{
recoChargedJetPartENoElecHist->Fill(tmpE);
recoChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
recoChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
recoChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
}
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++)
// Pairwise Distance Between Constituents
if(j<(partsEnd[i]-1))
{
TVector3 tmp(x[i],y[i],z[i]);
if(type[i] == 0)
for(unsigned int k=j+1; k<partsEnd[i]; k++)
{
//cout << i << " " << type[i] << " " << x[i] << " " << y[i] << " " << z[i] << " " << tmp.PseudoRapidity() << " " << tmp.Phi() << endl;
v.push_back(tmp);
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);
}
}
}
return v;
}; // Jet 3-Momentum
numRecoChargedJetPartsHist->Fill(partsEnd[i] - partsBegin[i]);
if(noElectron) numRecoChargedJetPartsNoElecHist->Fill(partsEnd[i] - partsBegin[i]);
}
auto getJetTrackDRThresh = [](ints &type, ints &pid, floats &E, floats &x, floats &y, floats &z){
floats dr;
for(auto i=0U; i<x.size(); i++)
// No Electrons
if(noElectron)
{
if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
recoChargedJetENoElecHist->Fill(recoNRG[i]);
if(ECut) recoChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
recoChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
if(ECut) recoChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
if(ECut) numRecoChargedJetsNoElec++;
}
}
numRecoChargedJetsECutHist->Fill(numRecoChargedJets);
numRecoChargedJetsECutNoElecHist->Fill(numRecoChargedJetsNoElec);
//////////////////////////////////////////////////////////////////////////
//////////////////////// Analyze Generator Jets ////////////////////////
//////////////////////////////////////////////////////////////////////////
int numGenChargedJets = 0;
int numGenChargedJetsNoElec = 0;
for(unsigned int i=0; i<genType.GetSize(); i++)
{
TVector3 jetMom(x[i],y[i],z[i]);
int jetID = pid[i];
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());
for(auto j=0U; j<x.size(); j++) // First loop over tracks
bool noElectron = true;
if(ECut)
{
if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
// Find Jets with Electrons
for(unsigned int m=genPartsBegin[i]; m<genPartsEnd[i]; m++)
{
TVector3 trackA(x[j],y[j],z[j]);
if(mcM[genPartIndex[m]] > 0.00050 && mcM[genPartIndex[m]] < 0.00052)
noElectron = false;
}
if(j < x.size() - 1) // Overflow protection
for(unsigned int j=genPartsBegin[i]; j<genPartsEnd[i]; j++)
{
// 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)
{
for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
genChargedJetPartENoElecHist->Fill(tmpE);
genChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
genChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
genChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
}
// Pairwise Distance Between Constituents
if(j<(genPartsEnd[i]-1))
{
if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
for(unsigned int k=j+1; k<genPartsEnd[i]; k++)
{
TVector3 trackB(x[k],y[k],z[k]);
double mXB = mcMomX[genPartIndex[k]];
double mYB = mcMomY[genPartIndex[k]];
double mZB = mcMomZ[genPartIndex[k]];
TVector3 partMomB(mXB,mYB,mZB);
float dEta = trackA.PseudoRapidity() - trackB.PseudoRapidity();
float dPhi = TVector2::Phi_mpi_pi(trackA.Phi() - trackB.Phi());
float dR = std::sqrt(dEta*dEta + dPhi*dPhi);
double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
dr.push_back(dR);
} // Track B
} // Second Loop
} // Overflow
} // Track A
} // First Loop
} // Jet Select
} // Jet Loop
return dr;
};
genChargedJetPartPairwiseDeltaRHist->Fill(dR);
}
}
}
numGenChargedJetPartsHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
if(noElectron) numGenChargedJetPartsNoElecHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
}
auto getJetTrackDR = [](ints &type, ints &pid, floats &x, floats &y, floats &z){
floats dr;
for(auto i=0U; i<x.size(); i++)
// No Electrons
if(noElectron)
{
if(type[i] == 0) // Look at jet
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++;
}
}
numGenChargedJetsECutHist->Fill(numGenChargedJets);
numGenChargedJetsECutNoElecHist->Fill(numGenChargedJetsNoElec);
//////////////////////////////////////////////////////////////////////////
///////////////////////////// Matched Jets /////////////////////////////
//////////////////////////////////////////////////////////////////////////
for(unsigned int i=0; i<genType.GetSize(); i++)
{
TVector3 jetMom(x[i],y[i],z[i]);
int jetID = pid[i];
TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
for(auto j=0U; j<x.size(); j++) // First loop over tracks
// 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++)
{
if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
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++)
{
TVector3 trackA(x[j],y[j],z[j]);
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(j < x.size() - 1) // Overflow protection
if(dR < minDeltaR)
{
for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
minDeltaR = dR;
minIndex = j;
}
}
// Do Backwards Match
double minDeltaRBack = 999.;
double minIndexBack = -1;
if(minIndex > -1)
{
if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
for(unsigned int j=0; j<genType.GetSize(); j++)
{
TVector3 trackB(x[k],y[k],z[k]);
TVector3 genMom(genMomX[j],genMomY[j],genMomZ[j]);
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
double dEta = recoMatchMom.PseudoRapidity() - genMom.PseudoRapidity();
double dPhi = TVector2::Phi_mpi_pi(recoMatchMom.Phi() - genMom.Phi());
double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
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(dR < minDeltaRBack)
{
if(type[i] == 0 && E[i] >= ETHRESH) // Look at jet with energy above thresh
minDeltaRBack = dR;
minIndexBack = 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)
{
int hasDupTrack = 0;
TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
TVector3 jetMom(x[i],y[i],z[i]);
int jetID = pid[i];
recoVsGenChargedJetENoDRHist->Fill(genNRG[i],recoNRG[minIndex]);
for(auto j=0U; j<x.size(); j++) // First loop over tracks
if(minDeltaR < DELTARCUT)
{
if(type[j] == 1 && (pid[j] == jetID)) // Look at Tracks associated with jet
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++)
{
TVector3 trackA(x[j],y[j],z[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);
if(j < x.size() - 1) // Overflow protection
{
for(auto k=j+1; k<x.size(); k++) // Second loop over tracks
TVector3 partMom(mX,mY,mZ);
// Pairwise Distance Between Constituents
if(j<(partsEnd[minIndex]-1))
{
if(type[k] == 1 && (pid[k] == jetID)) // Look at tracks associated with jet
for(unsigned int k=j+1; k<partsEnd[minIndex]; k++)
{
TVector3 trackB(x[k],y[k],z[k]);
double mXB = recoPartMomX[recoPartIndex[k]];
double mYB = recoPartMomY[recoPartIndex[k]];
double mZB = recoPartMomZ[recoPartIndex[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);
TVector3 partMomB(mXB,mYB,mZB);
if(dR < 0.05) hasDupTrack = 1;
} // Track B
} // Second Loop
} // Overflow
} // Track A
} // First Loop
double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
dup.push_back(hasDupTrack);
} // Jet Select
} // Jet Loop
return dup;
}; // Identify jets containing duplicate tracks
if(dR < 0.02) noDuplicate = false;
}
}
}
auto getPt = [](vecs &x){
floats pt;
for(auto i=0U; i<x.size(); i++)
if(noDuplicate)
{
pt.push_back(x[i].Perp());
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);
}
}
}
}
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());
NEVENTS++;
}
return eta;
}; // Jet Psudorapidity
auto getPhi = [](vecs &x){
floats phi;
for(auto i=0U; i<x.size(); i++)
gStyle->SetOptStat(0);
//////////////////////// Reconstructed Jets Plots ////////////////////////
// Reco Number
TCanvas *c1 = new TCanvas("c1","Number Reco Jets",800,600);
c1->Clear();
c1->Divide(1,1);
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");
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++)
{
phi.push_back(x[i].Phi());
binCent[i] = pxA->GetBinCenter(i+1);
}
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++)
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++)
{
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);
hA[i-1] = (TH1D *)jetResVsENegEtaNoDupHist->ProjectionY(Form("projYA_%d",i),i,i);
if(dR < minDeltaR) minDeltaR = dR;
}
deltaR.push_back(minDeltaR);
}
return deltaR;
};
TF1 *myGausA = new TF1("myGausA","gaus",-0.5,0.5);
myGausA->SetParameters(hA[i-1]->GetMaximum(),0.0,0.01);
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++)
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)
{
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);
auto fA = hA[i-1]->GetFunction("myGausA");
if(dR < minDeltaR)
if(fA->GetParError(2)/fA->GetParameter(2) < 0.5)
{
minDeltaR = dR;
minE = E[j];
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;
}
nrg.push_back(minE);
}
return nrg;
};
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;
auto getMatchEta = [](vecs &vecR, vecs &vecG){
floats eta;
for(auto i=0U; i<vecR.size(); i++)
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++)
{
float minDeltaR = 10000.0;
float minEta = -100.0;
for(auto j=0U; j<vecG.size(); j++)
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)
{
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);
auto fB = hB[i-1]->GetFunction("myGausB");
if(dR < minDeltaR)
if(fB->GetParError(2)/fB->GetParameter(2) < 0.5)
{
minDeltaR = dR;
minEta = vecG[j].PseudoRapidity();
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;
}
eta.push_back(minEta);
}
return eta;
};
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);
auto getMatchPhi = [](vecs &vecR, vecs &vecG){
floats phi;
for(auto i=0U; i<vecR.size(); i++)
TH1D *hC[20];
for(int i=1; i<21; i++)
{
float minDeltaR = 10000.0;
float minPhi = -100.0;
for(auto j=0U; j<vecG.size(); j++)
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)
{
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);
auto fC = hC[i-1]->GetFunction("myGausC");
if(dR < minDeltaR)
if(fC->GetParError(2)/fC->GetParameter(2) < 0.5)
{
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();
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;
}
benchmarks/dis/config.yml
+ 6
− 1
View file @ 0d19f33d
......@@ -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:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment