JetsAndHF/JetValidation/JetValidation.C

0001 // ============================================================================
0002 //! \file    JetValidation.C
0003 //! \authors Brian Page (bpage@bnl.gov),
0004 //!          adapted by Derek Anderson (derek.murphy.anderson@protonmail.com)
0005 // ----------------------------------------------------------------------------
0006 //! \brief Adaption of the Jet Benchmark to run
0007 //!   standalone to generate plots for validation.
0008 //!
0009 //! \usage In eic-shell:
0010 //!     root -b -q JetValidation.C'(<input file list>, \
0011 //!                                 <n files to read>, \
0012 //!                                 <output path>)'
0013 // ============================================================================
0014
0015 #include <edm4eic/EDM4eicVersion.h>
0016 #include <TCanvas.h>
0017 #include <TChain.h>
0018 #include <TFile.h>
0019 #include <TGraph.h>
0020 #include <TH1D.h>
0021 #include <TH2D.h>
0022 #include <TTree.h>
0023 #include <TTreeReader.h>
0024 #include <TTreeReaderArray.h>
0025 #include <TLegend.h>
0026 #include <TVector3.h>
0027 #include <fstream>
0028 #include <string>
0029 #include <vector>
0030
0031 #include "fmt/color.h"
0032 #include "fmt/core.h"
0033
0034 ///! Default input file list
0035 const std::string DefaultInFileList = "filelists/files26060.py8ncdis10x100q100t1000.list";
0036
0037 ///! Default no. of files
0038 const std::size_t DefaultNFiles = 1000;
0039
0040 ///! Default output file path
0041 const std::string DefaultOutPath = ".";
0042
0043 // ----------------------------------------------------------------------------
0044 // Does a branch exist?
0045 // ----------------------------------------------------------------------------
0046 /*! Checks if a branch exists in a TTree.
0047  *!
0048  *! \param[in] tree The tree to check
0049  *! \param[in] branch The branch to check for
0050  */
0051 bool branchExists(TTree* tree, const std::string& branch) {
0052   bool exists = false;
0053   if(tree->GetBranch(branch.c_str())) {
0054     exists = true;
0055   }
0056   return exists;
0057 }
0058
0059 // ----------------------------------------------------------------------------
0060 // Macro body
0061 // ----------------------------------------------------------------------------
0062 /*! Process input files to generate a set of reconstructed,
0063  *! generated jet distributions and save them as PNGs.
0064  *!
0065  *! \param[in]  filelist     Input filelist to use
0066  *! \param[in]  n_files      Number of files to read
0067  *! \param[out] results_path Location to save PNGs to
0068  */
0069 int JetValidation(
0070   const std::string& filelist = DefaultInFileList,
0071   const std::size_t n_files = DefaultNFiles,
0072   const std::string& results_path = DefaultOutPath
0073 ) {
0074
0075   const bool PRINT = true;
0076
0077   std::ifstream mylist(filelist);
0078   if (!mylist.is_open()) {
0079     std::cerr << "PANIC: Couldn't open fileslist '" << filelist << "'!" << std::endl;
0080     return 1;
0081   }
0082
0083   // Load input files
0084   std::string file;
0085   std::size_t i_file = 0;
0086   std::vector<std::string> rec_files;
0087   while (std::getline(mylist, file)) {
0088     rec_files.push_back(file);
0089     ++i_file;
0090     if (i_file == n_files) {
0091       break;
0092     }
0093   }
0094
0095   // Input
0096   TChain *mychain = new TChain("events");
0097   for (const auto& rec_file : rec_files) {
0098     mychain->Add(rec_file.c_str());
0099   }
0100
0101   const int seabornRed = TColor::GetColor(213, 94, 0);
0102
0103   // Seaborn Green: #009E73 -> (0, 158, 115)
0104   const int seabornGreen = TColor::GetColor(0, 158, 115);
0105
0106   // Seaborn Blue: #56B4E9 -> (86, 180, 233)
0107   const int seabornBlue = TColor::GetColor(100, 149, 237);
0108
0109   // TTreeReader
0110   TTreeReader tree_reader(mychain);
0111
0112   // Set Delta R Cut
0113   float DELTARCUT = 0.05;
0114
0115   // Check if area branch exists
0116   //   --> Using jet EDM if it does!
0117   bool useNewEDM = false;
0118 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8,9,0)
0119   {
0120     auto file = TFile::Open(rec_files.front().c_str(), "READ");
0121     auto tree = file->Get<TTree>("events");
0122     bool hasRecoArea = branchExists(tree, "ReconstructedChargedJets.area");
0123     bool hasGenArea = branchExists(tree, "GeneratedChargedJets.area");
0124     useNewEDM = hasRecoArea && hasGenArea;
0125   }
0126 #endif
0127   if (PRINT) {
0128     std::cout << "INFO: Using new EDM? " << useNewEDM << std::endl;
0129   }
0130
0131   // Set area branches to dummy values if not using jet EDM
0132   //   --> These branches won't be used if not using
0133   //       jet EDM
0134   std::string recoAreaBranch = "ReconstructedChargedJets.energy";
0135   std::string genAreaBranch = "GeneratedChargedJets.energy";
0136   if(useNewEDM) {
0137     recoAreaBranch = "ReconstructedChargedJets.area";
0138     genAreaBranch = "GeneratedChargedJets.area";
0139   }
0140
0141   // Use updated constituent branch names if using jet EDM
0142   std::string recoCstsBeginBranch = "ReconstructedChargedJets.particles_begin";
0143   std::string recoCstsEndBranch = "ReconstructedChargedJets.particles_end";
0144   std::string recoCstIndexBranch = "_ReconstructedChargedJets_particles.index";
0145   if (useNewEDM) {
0146     recoCstsBeginBranch = "ReconstructedChargedJets.constituents_begin";
0147     recoCstsEndBranch = "ReconstructedChargedJets.constituents_end";
0148     recoCstIndexBranch = "_ReconstructedChargedJets_constituents.index";
0149   }
0150
0151   std::string genCstsBeginBranch = "GeneratedChargedJets.particles_begin";
0152   std::string genCstsEndBranch = "GeneratedChargedJets.particles_end";
0153   std::string genCstIndexBranch = "_GeneratedChargedJets_particles.index";
0154   if (useNewEDM) {
0155     genCstsBeginBranch = "GeneratedChargedJets.constituents_begin";
0156     genCstsEndBranch = "GeneratedChargedJets.constituents_end";
0157     genCstIndexBranch = "_GeneratedChargedJets_constituents.index";
0158   }
0159
0160   // Reco Jets
0161   TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedJets.type"};
0162   TTreeReaderArray<float> recoNRG = {tree_reader, "ReconstructedChargedJets.energy"};
0163   TTreeReaderArray<float> recoMomX = {tree_reader, "ReconstructedChargedJets.momentum.x"};
0164   TTreeReaderArray<float> recoMomY = {tree_reader, "ReconstructedChargedJets.momentum.y"};
0165   TTreeReaderArray<float> recoMomZ = {tree_reader, "ReconstructedChargedJets.momentum.z"};
0166   TTreeReaderArray<float> recoArea = {tree_reader, recoAreaBranch.c_str()};
0167   TTreeReaderArray<unsigned int> recoCstsBegin = {tree_reader, recoCstsBeginBranch.c_str()};
0168   TTreeReaderArray<unsigned int> recoCstsEnd = {tree_reader, recoCstsEndBranch.c_str()};
0169
0170   TTreeReaderArray<int> recoCstIndex = {tree_reader, recoCstIndexBranch.c_str()};
0171
0172   // Reconstructed Particles
0173   TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0174   TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0175   TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0176   TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
0177   TTreeReaderArray<int> recoPartPDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
0178   TTreeReaderArray<float> recoPartNRG = {tree_reader, "ReconstructedChargedParticles.energy"};
0179
0180   TTreeReaderArray<int> recoPartAssocRec = {tree_reader, "_ReconstructedChargedParticleAssociations_rec.index"}; // Reco <-> MCParticle
0181   TTreeReaderArray<int> recoPartAssocSim = {tree_reader, "_ReconstructedChargedParticleAssociations_sim.index"};
0182   TTreeReaderArray<float> recoPartAssocWeight = {tree_reader, "ReconstructedChargedParticleAssociations.weight"};
0183
0184   // Generated Jets
0185   TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
0186   TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
0187   TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
0188   TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
0189   TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
0190   TTreeReaderArray<float> genArea = {tree_reader, genAreaBranch.c_str()};
0191   TTreeReaderArray<unsigned int> genCstsBegin = {tree_reader, genCstsBeginBranch.c_str()};
0192   TTreeReaderArray<unsigned int> genCstsEnd = {tree_reader, genCstsEndBranch.c_str()};
0193
0194   TTreeReaderArray<int> genPartIndex = {tree_reader, genCstIndexBranch.c_str()};
0195   //TTreeReaderArray<int> genChargedIndex = {tree_reader, "GeneratedChargedParticles_objIdx.index"};
0196
0197   // MC
0198   //TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0199   TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
0200   TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
0201   TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
0202   TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
0203   TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
0204
0205   TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0206   TTreeReaderArray<double> mcMomXPart = {tree_reader, "MCParticles.momentum.x"};
0207   TTreeReaderArray<double> mcMomYPart = {tree_reader, "MCParticles.momentum.y"};
0208   TTreeReaderArray<double> mcMomZPart = {tree_reader, "MCParticles.momentum.z"};
0209   TTreeReaderArray<double> mcMPart = {tree_reader, "MCParticles.mass"};
0210   TTreeReaderArray<int> pdgMCPart = {tree_reader, "MCParticles.PDG"};
0211
0212   // Define Histograms
0213   TH1D *counter = new TH1D("counter","",10,0.,10.);
0214
0215
0216   // Reco
0217   TH1D *numRecoChargedJetsECutHist = new TH1D("numRecoChargedJetsECut","",20,0.,20.);
0218   TH1D *recoChargedJetEHist = new TH1D("recoChargedJetE","",300,0.,300.);
0219   TH1D *recoChargedJetEtaECutHist = new TH1D("recoChargedJetEtaECut","",60,-3.,3.);
0220   TH1D *recoChargedJetAreaECutHist = nullptr;
0221   TH2D *recoChargedJetEvsAreaHist = nullptr;
0222   if (useNewEDM) {
0223     recoChargedJetAreaECutHist = new TH1D("recoChargedJetAreaECut","",250,0.,5.);
0224     recoChargedJetEvsAreaHist = new TH2D("recoChargedJetEvsArea","",250,0.,5.,300,0.,300.);
0225   }
0226   TH2D *recoChargedJetEvsEtaHist = new TH2D("recoChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0227   TH2D *recoChargedJetPhiVsEtaECutHist = new TH2D("recoChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0228
0229   TH1D *numRecoChargedJetsECutNoElecHist = new TH1D("numRecoChargedJetsECutNoElec","",20,0.,20.);
0230   TH1D *recoChargedJetENoElecHist = new TH1D("recoChargedJetENoElec","",300,0.,300.);
0231   TH1D *recoChargedJetEtaECutNoElecHist = new TH1D("recoChargedJetEtaECutNoElec","",60,-3.,3.);
0232   TH1D *recoChargedJetAreaECutNoElecHist = nullptr;
0233   TH2D *recoChargedJetEvsAreaNoElecHist = nullptr;
0234   if (useNewEDM) {
0235     recoChargedJetAreaECutNoElecHist = new TH1D("recoHargedJetAreaECutNoElec","",250,0.,5.);
0236     recoChargedJetEvsAreaNoElecHist = new TH2D("recoChargedJetEvsAreaNoElec","",250,0.,5.,300,0.,300.);
0237   }
0238   TH2D *recoChargedJetEvsEtaNoElecHist = new TH2D("recoChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0239   TH2D *recoChargedJetPhiVsEtaECutNoElecHist = new TH2D("recoChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0240
0241   TH1D *numRecoChargedJetPartsHist = new TH1D("numRecoChargedJetParts","",20,0.,20.);
0242   TH1D *recoChargedJetPartPHist = new TH1D("recoChargedJetPartP","",500,0.,100.);
0243   TH1D *recoChargedJetPartEtaHist = new TH1D("recoChargedJetPartEta","",80,-4.,4.);
0244   TH2D *recoChargedJetPartPvsEtaHist = new TH2D("recoChargedJetPartPvsEta","",80,-4.,4.,500,0.,100.);
0245   TH2D *recoChargedJetPartPhiVsEtaHist = new TH2D("recoChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0246
0247   TH1D *numRecoChargedJetPartsNoElecHist = new TH1D("numRecoChargedJetPartsNoElec","",20,0.,20.);
0248   TH1D *recoChargedJetPartPNoElecHist = new TH1D("recoChargedJetPartPNoElec","",500,0.,100.);
0249   TH1D *recoChargedJetPartEtaNoElecHist = new TH1D("recoChargedJetPartEtaNoElec","",80,-4.,4.);
0250   TH2D *recoChargedJetPartPvsEtaNoElecHist = new TH2D("recoChargedJetPartPvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0251   TH2D *recoChargedJetPartPhiVsEtaNoElecHist = new TH2D("recoChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0252
0253   TH1D *recoChargedJetPartPairwiseDeltaRHist = new TH1D("recoChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0254
0255   // Gen
0256   TH1D *numGenChargedJetsECutHist = new TH1D("numGenChargedJetsECut","",20,0.,20.);
0257   TH1D *genChargedJetEHist = new TH1D("genChargedJetE","",300,0.,300.);
0258   TH1D *genChargedJetEtaECutHist = new TH1D("genChargedJetEtaECut","",60,-3.,3.);
0259   TH1D *genChargedJetAreaECutHist = nullptr;
0260   TH2D *genChargedJetEvsAreaHist = nullptr;
0261   if (useNewEDM) {
0262     genChargedJetAreaECutHist = new TH1D("genChargedJetAreaECut","",250,0.,5.);
0263     genChargedJetEvsAreaHist = new TH2D("genChargedJetEvsAreaHist","",250,0.,5.,300,0.,300.);
0264   }
0265   TH2D *genChargedJetEvsEtaHist = new TH2D("genChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0266   TH2D *genChargedJetPhiVsEtaECutHist = new TH2D("genChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0267
0268   TH1D *numGenChargedJetsECutNoElecHist = new TH1D("numGenChargedJetsECutNoElec","",20,0.,20.);
0269   TH1D *genChargedJetENoElecHist = new TH1D("genChargedJetENoElec","",300,0.,300.);
0270   TH1D *genChargedJetEtaECutNoElecHist = new TH1D("genChargedJetEtaECutNoElec","",60,-3.,3.);
0271   TH1D *genChargedJetAreaECutNoElecHist = nullptr;
0272   TH2D *genChargedJetEvsAreaNoElecHist = nullptr;
0273   if (useNewEDM) {
0274     genChargedJetAreaECutNoElecHist = new TH1D("genChargedJetAreaECutNoElec","",250,0.,5.);
0275     genChargedJetEvsAreaNoElecHist = new TH2D("genChargedJetEvsAreaNoElec","",250,0.,5.,300,0.,300.);
0276   }
0277   TH2D *genChargedJetEvsEtaNoElecHist = new TH2D("genChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0278   TH2D *genChargedJetPhiVsEtaECutNoElecHist = new TH2D("genChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0279
0280   TH1D *numGenChargedJetPartsHist = new TH1D("numGenChargedJetParts","",20,0.,20.);
0281   TH1D *genChargedJetPartPHist = new TH1D("genChargedJetPartP","",500,0.,100.);
0282   TH1D *genChargedJetPartEtaHist = new TH1D("genChargedJetPartEta","",80,-4.,4.);
0283   TH2D *genChargedJetPartPvsEtaHist = new TH2D("genChargedJetPartPvsEta","",80,-4.,4.,500,0.,100.);
0284   TH2D *genChargedJetPartPhiVsEtaHist = new TH2D("genChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0285
0286   TH1D *numGenChargedJetPartsNoElecHist = new TH1D("numGenChargedJetPartsNoElec","",20,0.,20.);
0287   TH1D *genChargedJetPartPNoElecHist = new TH1D("genChargedJetPartPNoElec","",500,0.,100.);
0288   TH1D *genChargedJetPartEtaNoElecHist = new TH1D("genChargedJetPartEtaNoElec","",80,-4.,4.);
0289   TH2D *genChargedJetPartPvsEtaNoElecHist = new TH2D("genChargedJetPartPvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0290   TH2D *genChargedJetPartPhiVsEtaNoElecHist = new TH2D("genChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0291
0292   TH1D *genChargedJetPartPairwiseDeltaRHist = new TH1D("genChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0293
0294   // Matched
0295   TH1D *matchJetDeltaRHist = new TH1D("matchJetDeltaR","",5000,0.,5.);
0296   TH1D *matchJetDeltaRBackHist = new TH1D("matchJetDeltaRBack","",5000,0.,5.);
0297   TH2D *recoVsGenChargedJetEtaHist = new TH2D("recoVsGenChargedJetEta","",80,-4.,4.,80,-4.,4.);
0298   TH2D *recoVsGenChargedJetPhiHist = new TH2D("recoVsGenChargedJetPhi","",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
0299   TH2D *recoVsGenChargedJetAreaHist = nullptr;
0300   if (useNewEDM) {
0301     recoVsGenChargedJetAreaHist = new TH2D("recoVsGenChargedJetArea","",250,0.,5.,250,0.,5.);
0302   }
0303   TH2D *recoVsGenChargedJetEHist = new TH2D("recoVsGenChargedJetE","",100,0.,100.,100,0.,100.);
0304   TH2D *recoVsGenChargedJetENoDRHist = new TH2D("recoVsGenChargedJetENoDRHist","",100,0.,100.,100,0.,100.);
0305   TH2D *recoVsGenChargedJetENoDupHist = new TH2D("recoVsGenChargedJetENoDup","",100,0.,100.,100,0.,100.);
0306
0307   TH2D *jetResVsEtaHist = new TH2D("jetResVsEta","",80,-4.,4.,10000,-10.,10.);
0308   TH2D *jetResVsEHist = new TH2D("jetResVsE","",100,0.,100.,10000,-10.,10.);
0309   TH2D *jetResVsENegEtaHist = new TH2D("jetResVsENegEta","",20,0.,100.,10000,-10.,10.);
0310   TH2D *jetResVsEMidEtaHist = new TH2D("jetResVsEMidEta","",20,0.,100.,10000,-10.,10.);
0311   TH2D *jetResVsEPosEtaHist = new TH2D("jetResVsEPosEta","",20,0.,100.,10000,-10.,10.);
0312
0313   TH2D *jetResVsENegEtaNoDupHist = new TH2D("jetResVsENegEtaNoDup","",20,0.,100.,10000,-10.,10.);
0314   TH2D *jetResVsEMidEtaNoDupHist = new TH2D("jetResVsEMidEtaNoDup","",20,0.,100.,10000,-10.,10.);
0315   TH2D *jetResVsEPosEtaNoDupHist = new TH2D("jetResVsEPosEtaNoDup","",20,0.,100.,10000,-10.,10.);
0316
0317
0318   // Loop Through Events
0319   int NEVENTS = 0;
0320   while(tree_reader.Next()) {
0321
0322     if(NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
0323
0324     counter->Fill(0);
0325
0326     //////////////////////////////////////////////////////////////////////////
0327     //////////////////////  Analyze Reconstructed Jets  //////////////////////
0328     //////////////////////////////////////////////////////////////////////////
0329     int numRecoChargedJets = 0;
0330     int numRecoChargedJetsNoElec = 0;
0331     for(unsigned int i=0; i<recoType.GetSize(); i++)
0332       {
0333     TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
0334
0335     counter->Fill(3);
0336
0337     // Place eta cut to avoid edges of tracking acceptance
0338     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0339
0340     // Place a minimum energy condition for several plots
0341     bool ECut = recoNRG[i] > 5.0;
0342
0343     if(ECut) numRecoChargedJets++;
0344
0345     recoChargedJetEHist->Fill(recoNRG[i]);
0346     if(ECut) recoChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0347     recoChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0348         if(useNewEDM) {
0349           if(ECut) recoChargedJetAreaECutHist->Fill(recoArea[i]);
0350           recoChargedJetEvsAreaHist->Fill(recoArea[i],recoNRG[i]);
0351         }
0352     if(ECut) recoChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0353
0354     // Find Jets with Electrons
0355     bool noElectron = true;
0356     for(unsigned int m=recoCstsBegin[i]; m<recoCstsEnd[i]; m++) // Loop over jet constituents
0357       {
0358         int elecIndex = -1;
0359         double elecIndexWeight = -1.0;
0360         int chargePartIndex = recoCstIndex[m]; // ReconstructedChargedParticle Index for m'th Jet Component
0361         for(unsigned int n=0; n<recoPartAssocRec.GetSize(); n++) // Loop Over All ReconstructedChargedParticleAssociations
0362           {
0363         if(recoPartAssocRec[n] == chargePartIndex) // Select Entry Matching the ReconstructedChargedParticle Index
0364           {
0365             if(recoPartAssocWeight[n] > elecIndexWeight) // Find Particle with Greatest Weight = Contributed Most Hits to Track
0366               {
0367             elecIndex = recoPartAssocSim[n]; // Get Index of MCParticle Associated with ReconstructedChargedParticle
0368             elecIndexWeight = recoPartAssocWeight[n];
0369               }
0370           }
0371           }
0372
0373         if(pdgMCPart[elecIndex] == 11) // Test if Matched Particle is an Electron
0374           noElectron = false;
0375       }
0376
0377     if(ECut)
0378       {
0379         for(unsigned int j=recoCstsBegin[i]; j<recoCstsEnd[i]; j++)
0380           {
0381         // recoCstsBegin and recoCstsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
0382         // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
0383         double mX = recoPartMomX[recoCstIndex[j]];
0384         double mY = recoPartMomY[recoCstIndex[j]];
0385         double mZ = recoPartMomZ[recoCstIndex[j]];
0386         double mM = recoPartM[recoCstIndex[j]];
0387         //double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0388
0389         TVector3 partMom(mX,mY,mZ);
0390
0391         recoChargedJetPartPHist->Fill(partMom.Mag());
0392         recoChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0393         recoChargedJetPartPvsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0394         recoChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0395
0396         if(noElectron)
0397           {
0398             recoChargedJetPartPNoElecHist->Fill(partMom.Mag());
0399             recoChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0400             recoChargedJetPartPvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0401             recoChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0402           }
0403
0404         // Pairwise Distance Between Constituents
0405         if(j<(recoCstsEnd[i]-1))
0406           {
0407             for(unsigned int k=j+1; k<recoCstsEnd[i]; k++)
0408               {
0409             double mXB = recoPartMomX[recoCstIndex[k]];
0410             double mYB = recoPartMomY[recoCstIndex[k]];
0411             double mZB = recoPartMomZ[recoCstIndex[k]];
0412
0413             TVector3 partMomB(mXB,mYB,mZB);
0414
0415             double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0416             double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0417             double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0418
0419             recoChargedJetPartPairwiseDeltaRHist->Fill(dR);
0420               }
0421           }
0422           }
0423         numRecoChargedJetPartsHist->Fill(recoCstsEnd[i] - recoCstsBegin[i]);
0424         if(noElectron) numRecoChargedJetPartsNoElecHist->Fill(recoCstsEnd[i] - recoCstsBegin[i]);
0425       }
0426
0427     // No Electrons
0428     if(noElectron)
0429       {
0430         recoChargedJetENoElecHist->Fill(recoNRG[i]);
0431         if(ECut) recoChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0432         recoChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0433             if(useNewEDM) {
0434               if(ECut) recoChargedJetAreaECutNoElecHist->Fill(recoArea[i]);
0435               recoChargedJetEvsAreaNoElecHist->Fill(recoArea[i],recoNRG[i]);
0436             }
0437         if(ECut) recoChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0438
0439         if(ECut) numRecoChargedJetsNoElec++;
0440       }
0441       }
0442     numRecoChargedJetsECutHist->Fill(numRecoChargedJets);
0443     numRecoChargedJetsECutNoElecHist->Fill(numRecoChargedJetsNoElec);
0444
0445     //////////////////////////////////////////////////////////////////////////
0446     ////////////////////////  Analyze Generator Jets  ////////////////////////
0447     //////////////////////////////////////////////////////////////////////////
0448     int numGenChargedJets = 0;
0449     int numGenChargedJetsNoElec = 0;
0450     for(unsigned int i=0; i<genType.GetSize(); i++)
0451       {
0452     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0453
0454     counter->Fill(4);
0455
0456     // Place eta cut to avoid edges of tracking acceptance
0457     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0458
0459     // Place a minimum energy condition for several plots
0460     bool ECut = genNRG[i] > 5.0;
0461
0462     if(ECut) numGenChargedJets++;
0463
0464     genChargedJetEHist->Fill(genNRG[i]);
0465     if(ECut) genChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0466     genChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0467         if(useNewEDM) {
0468           if(ECut) genChargedJetAreaECutHist->Fill(genArea[i]);
0469           genChargedJetEvsAreaHist->Fill(genArea[i],genNRG[i]);
0470         }
0471     if(ECut) genChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0472
0473     // Find Jets with Electrons
0474     bool noElectron = true;
0475     for(unsigned int m=genCstsBegin[i]; m<genCstsEnd[i]; m++)
0476       {
0477         if(pdg[genPartIndex[m]] == 11)
0478           noElectron = false;
0479       }
0480
0481     if(ECut)
0482       {
0483         for(unsigned int j=genCstsBegin[i]; j<genCstsEnd[i]; j++)
0484           {
0485         // genCstsBegin and genCstsEnd specify the entries from _GeneratedChargedJets_particles.index that make up the jet
0486         // _GeneratedChargedJets_particles.index stores the GeneratedChargedParticles index of the jet constituent
0487         double mX = mcMomX[genPartIndex[j]];
0488         double mY = mcMomY[genPartIndex[j]];
0489         double mZ = mcMomZ[genPartIndex[j]];
0490         double mM = mcM[genPartIndex[j]];
0491         //double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0492
0493         TVector3 partMom(mX,mY,mZ);
0494
0495         genChargedJetPartPHist->Fill(partMom.Mag());
0496         genChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0497         genChargedJetPartPvsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0498         genChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0499
0500         if(noElectron)
0501           {
0502             genChargedJetPartPNoElecHist->Fill(partMom.Mag());
0503             genChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0504             genChargedJetPartPvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Mag());
0505             genChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0506           }
0507
0508         // Pairwise Distance Between Constituents
0509         if(j<(genCstsEnd[i]-1))
0510           {
0511             for(unsigned int k=j+1; k<genCstsEnd[i]; k++)
0512               {
0513             double mXB = mcMomX[genPartIndex[k]];
0514             double mYB = mcMomY[genPartIndex[k]];
0515             double mZB = mcMomZ[genPartIndex[k]];
0516
0517             TVector3 partMomB(mXB,mYB,mZB);
0518
0519             double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0520             double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0521             double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0522
0523             genChargedJetPartPairwiseDeltaRHist->Fill(dR);
0524               }
0525           }
0526           }
0527         numGenChargedJetPartsHist->Fill(genCstsEnd[i] - genCstsBegin[i]);
0528         if(noElectron) numGenChargedJetPartsNoElecHist->Fill(genCstsEnd[i] - genCstsBegin[i]);
0529       }
0530
0531     // No Electrons
0532     if(noElectron)
0533       {
0534         genChargedJetENoElecHist->Fill(genNRG[i]);
0535         if(ECut) genChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0536         genChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0537             if (useNewEDM) {
0538               if(ECut) genChargedJetAreaECutNoElecHist->Fill(genArea[i]);
0539               genChargedJetEvsAreaHist->Fill(genArea[i],genNRG[i]);
0540             }
0541         if(ECut) genChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0542
0543         if(ECut) numGenChargedJetsNoElec++;
0544       }
0545       }
0546     numGenChargedJetsECutHist->Fill(numGenChargedJets);
0547     numGenChargedJetsECutNoElecHist->Fill(numGenChargedJetsNoElec);
0548
0549
0550     //////////////////////////////////////////////////////////////////////////
0551     /////////////////////////////  Matched Jets  /////////////////////////////
0552     //////////////////////////////////////////////////////////////////////////
0553     for(unsigned int i=0; i<genType.GetSize(); i++)
0554       {
0555     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0556
0557     // Place eta cut to avoid edges of tracking acceptance
0558     //if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0559
0560     // Place a minimum energy condition
0561     //if(genNRG[i] < 5.0) continue;
0562
0563     // Don't Look at Electron Jets
0564     bool hasElectron = false;
0565     // Find Jets with Electrons
0566     for(unsigned int m=genCstsBegin[i]; m<genCstsEnd[i]; m++)
0567       {
0568         if(pdg[genPartIndex[m]] == 11)
0569           hasElectron = true;
0570       }
0571     //if(hasElectron) continue;
0572
0573     // Find Matching Reconstructed Jet
0574     double minDeltaR = 999.;
0575     int minIndex = -1;
0576     for(unsigned int j=0; j<recoType.GetSize(); j++)
0577       {
0578         TVector3 recoMom(recoMomX[j],recoMomY[j],recoMomZ[j]);
0579
0580         double dEta = jetMom.PseudoRapidity() - recoMom.PseudoRapidity();
0581         double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - recoMom.Phi());
0582         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0583
0584         if(dR < minDeltaR)
0585           {
0586         minDeltaR = dR;
0587         minIndex = j;
0588           }
0589       }
0590
0591     // Do Backwards Match
0592     double minDeltaRBack = 999.;
0593     double minIndexBack = -1;
0594     if(minIndex > -1)
0595       {
0596         TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0597         for(unsigned int j=0; j<genType.GetSize(); j++)
0598           {
0599         TVector3 genMom(genMomX[j],genMomY[j],genMomZ[j]);
0600
0601         double dEta = recoMatchMom.PseudoRapidity() - genMom.PseudoRapidity();
0602         double dPhi = TVector2::Phi_mpi_pi(recoMatchMom.Phi() - genMom.Phi());
0603         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0604
0605         if(dR < minDeltaRBack)
0606           {
0607             minDeltaRBack = dR;
0608             minIndexBack = j;
0609           }
0610           }
0611       }
0612
0613     // Look at Best Match
0614     if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1 && !hasElectron) matchJetDeltaRHist->Fill(minDeltaR);
0615     if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1) matchJetDeltaRBackHist->Fill(minDeltaR);
0616     if(minIndex > -1 && genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && !hasElectron)
0617       {
0618         TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0619
0620         recoVsGenChargedJetENoDRHist->Fill(genNRG[i],recoNRG[minIndex]);
0621
0622         if(minDeltaR < DELTARCUT)
0623           {
0624         recoVsGenChargedJetEtaHist->Fill(jetMom.PseudoRapidity(),recoMatchMom.PseudoRapidity());
0625         recoVsGenChargedJetPhiHist->Fill(jetMom.Phi(),recoMatchMom.Phi());
0626                 if(useNewEDM) {
0627                   recoVsGenChargedJetAreaHist->Fill(genArea[i],recoArea[minIndex]);
0628                 }
0629         recoVsGenChargedJetEHist->Fill(genNRG[i],recoNRG[minIndex]);
0630
0631         double jetERes = (recoNRG[minIndex] - genNRG[i])/genNRG[i];
0632
0633         jetResVsEtaHist->Fill(jetMom.PseudoRapidity(),jetERes);
0634         jetResVsEHist->Fill(genNRG[i],jetERes);
0635         if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0636           jetResVsENegEtaHist->Fill(genNRG[i],jetERes);
0637         if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0638           jetResVsEMidEtaHist->Fill(genNRG[i],jetERes);
0639         if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0640           jetResVsEPosEtaHist->Fill(genNRG[i],jetERes);
0641
0642         // Check for Duplicate Tracks
0643         bool noDuplicate = true;
0644         for(unsigned int j=recoCstsBegin[minIndex]; j<recoCstsEnd[minIndex]; j++)
0645           {
0646             double mX = recoPartMomX[recoCstIndex[j]];
0647             double mY = recoPartMomY[recoCstIndex[j]];
0648             double mZ = recoPartMomZ[recoCstIndex[j]];
0649             double mM = recoPartM[recoCstIndex[j]];
0650             double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0651
0652             TVector3 partMom(mX,mY,mZ);
0653
0654             // Pairwise Distance Between Constituents
0655             if(j<(recoCstsEnd[minIndex]-1))
0656               {
0657             for(unsigned int k=j+1; k<recoCstsEnd[minIndex]; k++)
0658               {
0659                 double mXB = recoPartMomX[recoCstIndex[k]];
0660                 double mYB = recoPartMomY[recoCstIndex[k]];
0661                 double mZB = recoPartMomZ[recoCstIndex[k]];
0662
0663                 TVector3 partMomB(mXB,mYB,mZB);
0664
0665                 double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0666                 double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0667                 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0668
0669                 if(dR < 0.02) noDuplicate = false;
0670               }
0671               }
0672           }
0673
0674         if(noDuplicate)
0675           {
0676             recoVsGenChargedJetENoDupHist->Fill(genNRG[i],recoNRG[minIndex]);
0677
0678             if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0679               jetResVsENegEtaNoDupHist->Fill(genNRG[i],jetERes);
0680             if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0681               jetResVsEMidEtaNoDupHist->Fill(genNRG[i],jetERes);
0682             if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0683               jetResVsEPosEtaNoDupHist->Fill(genNRG[i],jetERes);
0684           }
0685           }
0686       }
0687       }
0688
0689     NEVENTS++;
0690   }
0691
0692
0693   gStyle->SetOptStat(0);
0694   ////////////////////////  Reconstructed Jets Plots  ////////////////////////
0695   // Reco Number
0696   TCanvas *c1 = new TCanvas("c1","Number Reco Jets",800,600);
0697   c1->Clear();
0698   c1->Divide(1,1);
0699
0700   c1->cd(1);
0701   numRecoChargedJetsECutHist->Draw("HIST");
0702   numRecoChargedJetsECutNoElecHist->SetLineColor(seabornRed);
0703   numRecoChargedJetsECutNoElecHist->Draw("HISTSAME");
0704   numRecoChargedJetsECutHist->SetLineWidth(2); // Set line width to 2 (adjust as needed)
0705 numRecoChargedJetsECutNoElecHist->SetLineWidth(2);
0706   numRecoChargedJetsECutHist->SetTitle("Reconstructed Jets per Event (|eta| < 2.5 && E > 5);Number");
0707
0708 TLegend *legend1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0709 legend1->AddEntry(numRecoChargedJetsECutHist, "With Electrons", "l");
0710 legend1->AddEntry(numRecoChargedJetsECutNoElecHist, "No Electrons", "l");
0711 legend1->Draw();
0712
0713
0714
0715   gPad->SetLogy();
0716   if(PRINT) c1->Print((results_path+"/numberRecoJets.png").c_str()); // Number of reconstructed jets per event with energy > 5 GeV and Abs(eta) < 2.5
0717    delete c1;
0718
0719   // Reco Energy
0720   TCanvas *c2 = new TCanvas("c2","Reco Jet Energy",800,600);
0721   c2->Clear();
0722   c2->Divide(1,1);
0723
0724   c2->cd(1);
0725   recoChargedJetEHist->Draw("HIST");
0726   recoChargedJetENoElecHist->SetLineColor(seabornRed);
0727   recoChargedJetENoElecHist->Draw("HISTSAME");
0728
0729   recoChargedJetEHist->SetLineWidth(2);
0730   recoChargedJetENoElecHist->SetLineWidth(2);
0731   recoChargedJetEHist->SetTitle("Reconstructed Jet Energy (|eta| < 2.5);Energy [GeV]");
0732
0733 TLegend *legend2 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0734 legend2->AddEntry(recoChargedJetEHist, "With Electrons", "l");
0735 legend2->AddEntry(recoChargedJetENoElecHist, "No Electrons", "l");
0736 legend2->Draw();
0737
0738   gPad->SetLogy();
0739   if(PRINT) c2->Print((results_path+"/recoJetEnergy.png").c_str()); // Energy spectrum of reconstructed jets with Abs(eta) < 2.5
0740
0741     delete c2;
0742   // Reco Eta
0743   TCanvas *c3 = new TCanvas("c3","Reco Jet Eta",800,600);
0744   c3->Clear();
0745   c3->Divide(1,1);
0746
0747   c3->cd(1);
0748   recoChargedJetEtaECutHist->Draw("HIST");
0749   recoChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
0750   recoChargedJetEtaECutNoElecHist->Draw("HISTSAME");
0751
0752   recoChargedJetEtaECutHist->SetLineWidth(2);
0753   recoChargedJetEtaECutNoElecHist->SetLineWidth(2);
0754   recoChargedJetEtaECutHist->SetTitle("Reconstructed Jet Eta (E > 5);Eta");
0755
0756 //add legend
0757 TLegend *legend3 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0758 legend3->AddEntry(recoChargedJetEtaECutHist, "With Electrons", "l");
0759 legend3->AddEntry(recoChargedJetEtaECutNoElecHist, "No Electrons", "l");
0760 legend3->Draw();
0761
0762   gPad->SetLogy();
0763   if(PRINT) c3->Print((results_path+"/recoJetEta.png").c_str()); // Eta spectrum of reconstructed jets with energy > 5 GeV
0764     delete c3;
0765
0766   // Reco Area
0767   if(useNewEDM) {
0768     TCanvas *c3_1 = new TCanvas("c3_1","Reco Jet Area",800,600);
0769     c3_1->Clear();
0770     c3_1->Divide(1,1);
0771
0772     c3_1->cd(1);
0773     recoChargedJetAreaECutHist->Draw("HIST");
0774     recoChargedJetAreaECutNoElecHist->SetLineColor(seabornRed);
0775     recoChargedJetAreaECutNoElecHist->Draw("HISTSAME");
0776
0777     recoChargedJetAreaECutHist->SetLineWidth(2);
0778     recoChargedJetAreaECutNoElecHist->SetLineWidth(2);
0779     recoChargedJetAreaECutHist->SetTitle("Reconstructed Jet Area (E > 5);Area");
0780
0781     //add legend
0782     TLegend *legend3_1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0783     legend3_1->AddEntry(recoChargedJetAreaECutHist, "With Electrons", "l");
0784     legend3_1->AddEntry(recoChargedJetAreaECutNoElecHist, "No Electrons", "l");
0785     legend3_1->Draw();
0786
0787     gPad->SetLogy();
0788     if(PRINT) c3_1->Print((results_path+"/recoJetArea.png").c_str()); // Area spectrum of reconstructed jets with energy > 5 GeV
0789       delete c3_1;
0790    }
0791
0792   // Reco E Vs Eta
0793   TCanvas *c4 = new TCanvas("c4","Reco Jet E Vs Eta",800,600);
0794   c4->Clear();
0795   c4->Divide(1,1);
0796
0797   c4->cd(1);
0798   recoChargedJetEvsEtaHist->Draw("COLZ");
0799   recoChargedJetEvsEtaHist->SetTitle("Reconstructed Jet Energy Vs Eta;Eta;Energy [GeV]");
0800   gPad->SetLogz();
0801   if(PRINT) c4->Print((results_path+"/recoJetEnergyvsEta.png").c_str()); // Energy vs eta of reconstructed jets
0802
0803   // Reco E Vs Area
0804   if (useNewEDM) {
0805     TCanvas *c4_1 = new TCanvas("c4_1","Reco Jet E Vs Area",800,600);
0806     c4_1->Clear();
0807     c4_1->Divide(1,1);
0808
0809     c4_1->cd(1);
0810     recoChargedJetEvsAreaHist->Draw("COLZ");
0811     recoChargedJetEvsAreaHist->SetTitle("Reconstructed Jet Energy Vs Area;Area;Energy [GeV]");
0812     gPad->SetLogz();
0813     if(PRINT) c4_1->Print((results_path+"/recoJetEnergyvsArea.png").c_str()); // Energy vs area of reconstructed jets
0814   }
0815
0816   // Reco Phi Vs Eta
0817   TCanvas *c5 = new TCanvas("c5","Reco Jet Phi Vs Eta",800,600);
0818   c5->Clear();
0819   c5->Divide(1,1);
0820
0821   c5->cd(1);
0822   recoChargedJetPhiVsEtaECutHist->Draw("COLZ");
0823   recoChargedJetPhiVsEtaECutHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5);Eta;Phi");
0824   gPad->SetLogz();
0825   if(PRINT) c5->Print((results_path+"/recoJetPhivsEta.png").c_str()); // Phi vs eta of reconstructed jets
0826
0827   // Num Particles Per Reco Jet
0828   TCanvas *c6 = new TCanvas("c6","Number Constituents Per Reco Jet",800,600);
0829   c6->Clear();
0830   c6->Divide(1,1);
0831
0832   c6->cd(1);
0833   numRecoChargedJetPartsHist->Draw("HIST");
0834   numRecoChargedJetPartsNoElecHist->SetLineColor(seabornRed);
0835   numRecoChargedJetPartsNoElecHist->Draw("HISTSAME");
0836
0837   numRecoChargedJetPartsHist->SetLineWidth(2);
0838   numRecoChargedJetPartsNoElecHist->SetLineWidth(2);
0839   numRecoChargedJetPartsHist->SetTitle("Number of Constituents Per Reco Jet;Number of Constituents");
0840
0841 TLegend *legend6 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0842 legend6->AddEntry(numRecoChargedJetPartsHist, "With Electrons", "l");
0843 legend6->AddEntry(numRecoChargedJetPartsNoElecHist, "No Electrons", "l");
0844 legend6->Draw();
0845
0846   gPad->SetLogy();
0847   if(PRINT) c6->Print((results_path+"/numConstituentsPerRecoJet.png").c_str()); // Number of constituents in reconstructed jets
0848
0849   // Reco Part Energy
0850   TCanvas *c7 = new TCanvas("c7","Reco Jet Constituent Momentum",800,600);
0851   c7->Clear();
0852   c7->Divide(1,1);
0853
0854   c7->cd(1);
0855   recoChargedJetPartPHist->Draw("HIST");
0856   recoChargedJetPartPNoElecHist->SetLineColor(seabornRed);
0857   recoChargedJetPartPNoElecHist->Draw("HISTSAME");
0858
0859   recoChargedJetPartPHist->SetLineWidth(2);
0860   recoChargedJetPartPNoElecHist->SetLineWidth(2);
0861   recoChargedJetPartPHist->SetTitle("Reconstructed Jet Constituent Momentum;Momentum [GeV/c]");
0862
0863   TLegend *legend7 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0864   legend7->AddEntry(recoChargedJetPartPHist, "With Electrons", "l");
0865   legend7->AddEntry(recoChargedJetPartPNoElecHist, "No Electrons", "l");
0866   legend7->Draw();
0867
0868   gPad->SetLogy();
0869   if(PRINT) c7->Print((results_path+"/recoJetConstituentMomentum.png").c_str()); // Momentum of reconstructed jet constituents
0870
0871   // Reco Part Eta
0872   TCanvas *c8 = new TCanvas("c8","Reco Jet Constituent Eta",800,600);
0873   c8->Clear();
0874   c8->Divide(1,1);
0875
0876   c8->cd(1);
0877   recoChargedJetPartEtaHist->Draw("HIST");
0878   recoChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
0879   recoChargedJetPartEtaNoElecHist->Draw("HISTSAME");
0880
0881   recoChargedJetPartEtaHist->SetLineWidth(2);
0882   recoChargedJetPartEtaNoElecHist->SetLineWidth(2);
0883
0884   recoChargedJetPartEtaHist->SetTitle("Reconstructed Jet Constituent Eta;Eta");
0885
0886   TLegend *legend8 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0887   legend8->AddEntry(recoChargedJetPartEtaHist, "With Electrons", "l");
0888   legend8->AddEntry(recoChargedJetPartEtaNoElecHist, "No Electrons", "l");
0889   legend8->Draw();
0890
0891   gPad->SetLogy();
0892   if(PRINT) c8->Print((results_path+"/recoJetConstituentEta.png").c_str()); // Eta of reconstructed jet constituents
0893
0894   // Reco Part P Vs Eta
0895   TCanvas *c9 = new TCanvas("c9","Reco Jet Constituent Momentum Vs Eta",800,600);
0896   c9->Clear();
0897   c9->Divide(1,1);
0898
0899   c9->cd(1);
0900   recoChargedJetPartPvsEtaHist->Draw("COLZ");
0901   recoChargedJetPartPvsEtaHist->SetTitle("Reconstructed Jet Constituent Momentum Vs Eta;Eta;Momentum [GeV/c]");
0902   gPad->SetLogz();
0903   if(PRINT) c9->Print((results_path+"/recoJetConstituentMomentumVsEta.png").c_str()); // Momentum vs eta of reconstructed jet constituents
0904
0905   // Reco Part Phi Vs Eta
0906   TCanvas *c10 = new TCanvas("c10","Reco Jet Constituent Phi Vs Eta",800,600);
0907   c10->Clear();
0908   c10->Divide(1,1);
0909
0910   c10->cd(1);
0911   recoChargedJetPartPhiVsEtaHist->Draw("COLZ");
0912   recoChargedJetPartPhiVsEtaHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta;Eta;Phi");
0913   gPad->SetLogz();
0914   if(PRINT) c10->Print((results_path+"/recoJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of reconstructed jet constituents
0915
0916   // Reco Constituent Pairwise delta R
0917   TCanvas *c11 = new TCanvas("c11","Reco Jet Constituent Pairwise Delta R",800,600);
0918   c11->Clear();
0919   c11->Divide(1,1);
0920
0921   c11->cd(1);
0922   recoChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
0923   recoChargedJetPartPairwiseDeltaRHist->SetTitle("Pairwise Constituent Delta R;Delta R");
0924   recoChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
0925   gPad->SetLogy();
0926   if(PRINT) c11->Print((results_path+"/recoJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in reconstructed jets
0927
0928   // Reco E Vs Eta No Electron Jets
0929   TCanvas *c12 = new TCanvas("c12","Reco Jet E Vs Eta (No Electrons)",800,600);
0930   c12->Clear();
0931   c12->Divide(1,1);
0932
0933   c12->cd(1);
0934   recoChargedJetEvsEtaNoElecHist->Draw("COLZ");
0935   recoChargedJetEvsEtaNoElecHist->SetTitle("Reconstructed Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
0936   gPad->SetLogz();
0937   if(PRINT) c12->Print((results_path+"/recoJetEnergyVsEtaNoElectron.png").c_str()); // Reconstructed jet energy - no jets containing electrons included
0938
0939   // Reco Phi Vs Eta No Electron Jets
0940   TCanvas *c13 = new TCanvas("c13","Reco Jet Phi Vs Eta (No Electrons)",800,600);
0941   c13->Clear();
0942   c13->Divide(1,1);
0943
0944   c13->cd(1);
0945   recoChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
0946   recoChargedJetPhiVsEtaECutNoElecHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
0947   gPad->SetLogz();
0948   if(PRINT) c13->Print((results_path+"/recoJetPhiVsEtaNoElectron.png").c_str()); // Reconstructed Jet phi vs eta - no jets containing electrons included
0949
0950   // Reco Part P Vs Eta No Electron Jets
0951   TCanvas *c14 = new TCanvas("c14","Reco Jet Constituent Momentum Vs Eta (No Electrons)",800,600);
0952   c14->Clear();
0953   c14->Divide(1,1);
0954
0955   c14->cd(1);
0956   recoChargedJetPartPvsEtaNoElecHist->Draw("COLZ");
0957   recoChargedJetPartPvsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Momentum Vs Eta (No Electrons);Eta;Momentum [GeV/c]");
0958   gPad->SetLogz();
0959   if(PRINT) c14->Print((results_path+"/recoJetConstituentMomentumVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent momentum vs eta - no jets containing electrons included
0960
0961   // Reco Part Phi Vs Eta No Electron Jets
0962   TCanvas *c15 = new TCanvas("c15","Reco Jet Constituent Phi Vs Eta (No Electrons)",800,600);
0963   c15->Clear();
0964   c15->Divide(1,1);
0965
0966   c15->cd(1);
0967   recoChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
0968   recoChargedJetPartPhiVsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
0969   gPad->SetLogz();
0970   if(PRINT) c15->Print((results_path+"/recoJetConstituentPhiVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent phi vs eta - no jets containing electrons included
0971
0972
0973   ////////////////////////  Generated Jets Plots  ////////////////////////
0974   // Gen Number
0975   TCanvas *c16 = new TCanvas("c16","Number Gen Jets",800,600);
0976   c16->Clear();
0977   c16->Divide(1,1);
0978
0979   c16->cd(1);
0980   numGenChargedJetsECutHist->Draw("HIST");
0981   numGenChargedJetsECutNoElecHist->SetLineColor(seabornRed);
0982   numGenChargedJetsECutNoElecHist->Draw("HISTSAME");
0983
0984   numGenChargedJetsECutHist->SetLineWidth(2);
0985   numGenChargedJetsECutNoElecHist->SetLineWidth(2);
0986
0987   numGenChargedJetsECutHist->SetTitle("Generator Jets per Event (|eta| < 2.5 && E > 5);Number");
0988
0989   TLegend *legend16 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0990   legend16->AddEntry(numGenChargedJetsECutHist, "With Electrons", "l");
0991   legend16->AddEntry(numGenChargedJetsECutNoElecHist, "No Electrons", "l");
0992   legend16->Draw();
0993
0994   gPad->SetLogy();
0995   if(PRINT) c16->Print((results_path+"/numberGenJets.png").c_str()); // Number of generator jets per event with energy > 5 GeV and Abs(eta) < 2.5
0996
0997   // Gen Energy
0998   TCanvas *c17 = new TCanvas("c17","Gen Jet Energy",800,600);
0999   c17->Clear();
1000   c17->Divide(1,1);
1001
1002   c17->cd(1);
1003   genChargedJetEHist->Draw("HIST");
1004   genChargedJetENoElecHist->SetLineColor(seabornRed);
1005   genChargedJetENoElecHist->Draw("HISTSAME");
1006
1007   genChargedJetEHist->SetLineWidth(2);
1008   genChargedJetENoElecHist->SetLineWidth(2);
1009
1010   genChargedJetEHist->SetTitle("Generator Jet Energy (|eta| < 2.5);Energy [GeV]");
1011   TLegend *legend17 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1012   legend17->AddEntry(genChargedJetEHist, "With Electrons", "l");
1013   legend17->AddEntry(genChargedJetENoElecHist, "No Electrons", "l");
1014   legend17->Draw();
1015
1016   gPad->SetLogy();
1017   if(PRINT) c17->Print((results_path+"/genJetEnergy.png").c_str()); // Energy spectrum of generated jets with Abs(eta) < 2.5
1018
1019   // Gen Eta
1020   TCanvas *c18 = new TCanvas("c18","Gen Jet Eta",800,600);
1021   c18->Clear();
1022   c18->Divide(1,1);
1023
1024   c18->cd(1);
1025   genChargedJetEtaECutHist->Draw("HIST");
1026   genChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
1027   genChargedJetEtaECutNoElecHist->Draw("HISTSAME");
1028
1029   genChargedJetEtaECutHist->SetLineWidth(2);
1030   genChargedJetEtaECutNoElecHist->SetLineWidth(2);
1031
1032   genChargedJetEtaECutHist->SetTitle("Generator Jet Eta (E > 5);Eta");
1033
1034   TLegend *legend18 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1035   legend18->AddEntry(genChargedJetEtaECutHist, "With Electrons", "l");
1036   legend18->AddEntry(genChargedJetEtaECutNoElecHist, "No Electrons", "l");
1037   legend18->Draw();
1038
1039   gPad->SetLogy();
1040   if(PRINT) c18->Print((results_path+"/genJetEta.png").c_str()); // Eta spectrum of generator jets with energy > 5 GeV
1041
1042   // Gen Area
1043   if(useNewEDM) {
1044     TCanvas *c18_1 = new TCanvas("c18_1","Gen Jet Area",800,600);
1045     c18_1->Clear();
1046     c18_1->Divide(1,1);
1047
1048     c18_1->cd(1);
1049     genChargedJetAreaECutHist->Draw("HIST");
1050     genChargedJetAreaECutNoElecHist->SetLineColor(seabornRed);
1051     genChargedJetAreaECutNoElecHist->Draw("HISTSAME");
1052
1053     genChargedJetAreaECutHist->SetLineWidth(2);
1054     genChargedJetAreaECutNoElecHist->SetLineWidth(2);
1055
1056     genChargedJetAreaECutHist->SetTitle("Generator Jet Area (E > 5);Area");
1057
1058     TLegend *legend18_1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1059     legend18_1->AddEntry(genChargedJetAreaECutHist, "With Electrons", "l");
1060     legend18_1->AddEntry(genChargedJetAreaECutNoElecHist, "No Electrons", "l");
1061     legend18_1->Draw();
1062
1063     gPad->SetLogy();
1064     if(PRINT) c18_1->Print((results_path+"/genJetArea.png").c_str()); // Area spectrum of generator jets with energy > 5 GeV
1065   }
1066
1067   // Gen E Vs Eta
1068   TCanvas *c19 = new TCanvas("c19","Gen Jet E Vs Eta",800,600);
1069   c19->Clear();
1070   c19->Divide(1,1);
1071
1072   c19->cd(1);
1073   genChargedJetEvsEtaHist->Draw("COLZ");
1074   genChargedJetEvsEtaHist->SetTitle("Generator Jet Energy Vs Eta;Eta;Energy [GeV]");
1075   gPad->SetLogz();
1076   if(PRINT) c19->Print((results_path+"/genJetEnergyvsEta.png").c_str()); // Energy vs eta of generator jets
1077
1078   // Gen E Vs Area
1079   if(useNewEDM) {
1080     TCanvas *c19_1 = new TCanvas("c19_1","Gen Jet E Vs Area",800,600);
1081     c19_1->Clear();
1082     c19_1->Divide(1,1);
1083
1084     c19_1->cd(1);
1085     genChargedJetEvsAreaHist->Draw("COLZ");
1086     genChargedJetEvsAreaHist->SetTitle("Generator Jet Energy Vs Area;Area;Energy [GeV]");
1087     gPad->SetLogz();
1088     if(PRINT) c19_1->Print((results_path+"/genJetEnergyvsArea.png").c_str()); // Energy vs area of generator jets
1089   }
1090
1091   // Gen Phi Vs Eta
1092   TCanvas *c20 = new TCanvas("c20","Gen Jet Phi Vs Eta",800,600);
1093   c20->Clear();
1094   c20->Divide(1,1);
1095
1096   c20->cd(1);
1097   genChargedJetPhiVsEtaECutHist->Draw("COLZ");
1098   genChargedJetPhiVsEtaECutHist->SetTitle("Generator Jet Phi Vs Eta (E > 5);Eta;Phi");
1099   gPad->SetLogz();
1100   if(PRINT) c20->Print((results_path+"/genJetPhiVsEta.png").c_str()); // Phi vs eta of generator jets
1101
1102   // Num Particles Per Gen Jet
1103   TCanvas *c21 = new TCanvas("c21","Number Constituents Per Gen Jet",800,600);
1104   c21->Clear();
1105   c21->Divide(1,1);
1106
1107   c21->cd(1);
1108   numGenChargedJetPartsHist->Draw("HIST");
1109   numGenChargedJetPartsNoElecHist->SetLineColor(seabornRed);
1110   numGenChargedJetPartsNoElecHist->Draw("HISTSAME");
1111
1112   numGenChargedJetPartsHist->SetLineWidth(2);
1113   numGenChargedJetPartsNoElecHist->SetLineWidth(2);
1114
1115   numGenChargedJetPartsHist->SetTitle("Number of Constituents Per Gen Jet;Number of Constituents");
1116
1117   TLegend *legend21 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1118   legend21->AddEntry(numGenChargedJetPartsHist, "With Electrons", "l");
1119   legend21->AddEntry(numGenChargedJetPartsNoElecHist, "No Electrons", "l");
1120   legend21->Draw();
1121   gPad->SetLogy();
1122   if(PRINT) c21->Print((results_path+"/numConstituentsPerGenJet.png").c_str()); // Number of constituents in generator jets
1123
1124   // Gen Part Momentum
1125   TCanvas *c22 = new TCanvas("c22","Gen Jet Constituent Momentum",800,600);
1126   c22->Clear();
1127   c22->Divide(1,1);
1128
1129   c22->cd(1);
1130   genChargedJetPartPHist->Draw("HIST");
1131   genChargedJetPartPNoElecHist->SetLineColor(seabornRed);
1132   genChargedJetPartPNoElecHist->Draw("HISTSAME");
1133
1134   genChargedJetPartPHist->SetLineWidth(2);
1135   genChargedJetPartPNoElecHist->SetLineWidth(2);
1136
1137   genChargedJetPartPHist->SetTitle("Generator Jet Constituent Momentum;Momentum [GeV/c]");
1138
1139   TLegend *legend22 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1140   legend22->AddEntry(genChargedJetPartPHist, "With Electrons", "l");
1141   legend22->AddEntry(genChargedJetPartPNoElecHist, "No Electrons", "l");
1142   legend22->Draw();
1143
1144   gPad->SetLogy();
1145   if(PRINT) c22->Print((results_path+"/genJetConstituentMomentum.png").c_str()); // Momentum of generator jet constituents
1146
1147   // Gen Part Eta
1148   TCanvas *c23 = new TCanvas("c23","Gen Jet Constituent Eta",800,600);
1149   c23->Clear();
1150   c23->Divide(1,1);
1151
1152   c23->cd(1);
1153   genChargedJetPartEtaHist->Draw("HIST");
1154   genChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
1155   genChargedJetPartEtaNoElecHist->Draw("HISTSAME");
1156
1157   genChargedJetPartEtaHist->SetLineWidth(2);
1158   genChargedJetPartEtaNoElecHist->SetLineWidth(2);
1159
1160   genChargedJetPartEtaHist->SetTitle("Generator Jet Constituent Eta;Eta");
1161
1162   TLegend *legend23 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
1163   legend23->AddEntry(genChargedJetPartEtaHist, "With Electrons", "l");
1164   legend23->AddEntry(genChargedJetPartEtaNoElecHist, "No Electrons", "l");
1165   legend23->Draw();
1166
1167   gPad->SetLogy();
1168   if(PRINT) c23->Print((results_path+"/genJetConstituentEta.png").c_str()); // Eta of generator jet constituents
1169
1170   // Gen Part P Vs Eta
1171   TCanvas *c24 = new TCanvas("c24","Gen Jet Constituent Momentum Vs Eta",800,600);
1172   c24->Clear();
1173   c24->Divide(1,1);
1174
1175   c24->cd(1);
1176   genChargedJetPartPvsEtaHist->Draw("COLZ");
1177   genChargedJetPartPvsEtaHist->SetTitle("Generator Jet Constituent Momentum Vs Eta;Eta;Momentum [GeV/c]");
1178   gPad->SetLogz();
1179   if(PRINT) c24->Print((results_path+"/genJetConstituentMomentumVsEta.png").c_str()); // Momentum vs eta of generator jet constituents
1180
1181   // Gen Part Phi Vs Eta
1182   TCanvas *c25 = new TCanvas("c25","Gen Jet Constituent Phi Vs Eta",800,600);
1183   c25->Clear();
1184   c25->Divide(1,1);
1185
1186   c25->cd(1);
1187   genChargedJetPartPhiVsEtaHist->Draw("COLZ");
1188   genChargedJetPartPhiVsEtaHist->SetTitle("Generator Jet Constituent Phi Vs Eta;Eta;Phi");
1189   gPad->SetLogz();
1190   if(PRINT) c25->Print((results_path+"/genJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of generator jet constituents
1191
1192   // Gen Constituent Pairwise delta R
1193   TCanvas *c26 = new TCanvas("c26","Gen Jet Constituent Pairwise Delta R",800,600);
1194   c26->Clear();
1195   c26->Divide(1,1);
1196
1197   c26->cd(1);
1198   genChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
1199   genChargedJetPartPairwiseDeltaRHist->SetTitle("Generator Jet Pairwise Constituent Delta R;Delta R");
1200   genChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
1201   gPad->SetLogy();
1202   if(PRINT) c26->Print((results_path+"/genJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in generator jets
1203
1204   // Gen E Vs Eta No Electron Jets
1205   TCanvas *c27 = new TCanvas("c27","Gen Jet E Vs Eta (No Electrons)",800,600);
1206   c27->Clear();
1207   c27->Divide(1,1);
1208
1209   c27->cd(1);
1210   genChargedJetEvsEtaNoElecHist->Draw("COLZ");
1211   genChargedJetEvsEtaNoElecHist->SetTitle("Generator Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
1212   gPad->SetLogz();
1213   if(PRINT) c27->Print((results_path+"/genJetEnergyVsEtaNoElectron.png").c_str()); // Generator jet energy vs eta - no jets containing electrons included
1214
1215   // Gen Phi Vs Eta No Electron Jets
1216   TCanvas *c28 = new TCanvas("c28","Gen Jet Phi Vs Eta (No Electrons)",800,600);
1217   c28->Clear();
1218   c28->Divide(1,1);
1219
1220   c28->cd(1);
1221   genChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
1222   genChargedJetPhiVsEtaECutNoElecHist->SetTitle("Generator Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
1223   gPad->SetLogz();
1224   if(PRINT) c28->Print((results_path+"/genJetPhiVsEtaNoElectron.png").c_str()); // Generator Jet phi vs eta - no jets containing electrons included
1225
1226   // Gen Part P Vs Eta No Electron Jets
1227   TCanvas *c29 = new TCanvas("c29","Gen Jet Constituent Momentum Vs Eta (No Electrons)",800,600);
1228   c29->Clear();
1229   c29->Divide(1,1);
1230
1231   c29->cd(1);
1232   genChargedJetPartPvsEtaNoElecHist->Draw("COLZ");
1233   genChargedJetPartPvsEtaNoElecHist->SetTitle("Generator Jet Constituent Momentum Vs Eta (No Electrons);Eta;Momentum [GeV/c]");
1234   gPad->SetLogz();
1235   if(PRINT) c29->Print((results_path+"/genJetConstituentMomentumVsEtaNoElectron.png").c_str()); // Generator jet constituent momentum vs eta - no jets containing electrons included
1236
1237   // Gen Part Phi Vs Eta No Electron Jets
1238   TCanvas *c30 = new TCanvas("c30","Gen Jet Constituent Phi Vs Eta (No Electrons)",800,600);
1239   c30->Clear();
1240   c30->Divide(1,1);
1241
1242   c30->cd(1);
1243   genChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
1244   genChargedJetPartPhiVsEtaNoElecHist->SetTitle("Generator Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
1245   gPad->SetLogz();
1246   //c30->Print((results_path+"/recoJetEvsEta.png").c_str());
1247   if(PRINT) c30->Print((results_path+"/genJetConstituentPhiVsEtaNoElectron.png").c_str()); // Generator jet constituent phi vs eta - no jets containing electrons included
1248
1249
1250   ////////////////////////  Matched Jets Plots  ////////////////////////
1251   // Matched Delta R
1252   TCanvas *c31 = new TCanvas("c31","Gen - Reco Delta R",800,600);
1253   c31->Clear();
1254   c31->Divide(1,1);
1255
1256   c31->cd(1);
1257   matchJetDeltaRHist->Draw("HIST");
1258   matchJetDeltaRBackHist->SetLineColor(seabornRed);
1259   //matchJetDeltaRBackHist->Draw("HISTSAME");
1260   matchJetDeltaRHist->SetTitle("Matched Gen - Reco Jet Delta R;Delta R");
1261   gPad->SetLogy();
1262   if(PRINT) c31->Print((results_path+"/genRecoJetDeltaR.png").c_str()); // Distance between closest generated and reconstructed jet pair
1263
1264   // Matched Reco Vs Gen Eta
1265   TCanvas *c32 = new TCanvas("c32","Reco Vs Gen Eta",800,600);
1266   c32->Clear();
1267   c32->Divide(1,1);
1268
1269   c32->cd(1);
1270   recoVsGenChargedJetEtaHist->Draw("COLZ");
1271   recoVsGenChargedJetEtaHist->SetTitle("Reconstructed Vs Generator Jet Eta;Gen Eta;Reco Eta");
1272   gPad->SetLogz();
1273   if(PRINT) c32->Print((results_path+"/matchedRecoVsGenJetEta.png").c_str()); // Matched Reconstructed Vs Generator Jet eta
1274
1275   // Matched Reco Vs Gen Phi
1276   TCanvas *c33 = new TCanvas("c33","Reco Vs Gen Phi",800,600);
1277   c33->Clear();
1278   c33->Divide(1,1);
1279
1280   c33->cd(1);
1281   recoVsGenChargedJetPhiHist->Draw("COLZ");
1282   recoVsGenChargedJetPhiHist->SetTitle("Reconstructed Vs Generator Jet Phi;Gen Phi;Reco Phi");
1283   gPad->SetLogz();
1284   if(PRINT) c33->Print((results_path+"/matchedRecoVsGenJetPhi.png").c_str()); // Matched reconstructed vs generator jet phi
1285
1286   // Matched Reco Vs Gen Area
1287   if(useNewEDM) {
1288     TCanvas *c33_1 = new TCanvas("c33_1","Reco Vs Gen Area",800,600);
1289     c33_1->Clear();
1290     c33_1->Divide(1,1);
1291
1292     c33_1->cd(1);
1293     recoVsGenChargedJetAreaHist->Draw("COLZ");
1294     recoVsGenChargedJetAreaHist->SetTitle("Reconstructed Vs Generator Jet Area;Gen Area;Reco Area");
1295     gPad->SetLogz();
1296     if(PRINT) c33_1->Print((results_path+"/matchedRecoVsGenJetArea.png").c_str()); // Matched reconstructed vs generator jet area
1297   }
1298
1299   // Matched Reco Vs Gen Energy
1300   TCanvas *c34 = new TCanvas("c34","Reco Vs Gen Energy",800,600);
1301   c34->Clear();
1302   c34->Divide(1,1);
1303
1304   TF1 *f1_34 = new TF1("f1_34","1.0*x + 0.0",1,100);
1305   TF1 *f2_34 = new TF1("f2_34","2.0*x + 0.0",1,100);
1306   TF1 *f3_34 = new TF1("f3_34","3.0*x + 0.0",1,100);
1307
1308   c34->cd(1);
1309   recoVsGenChargedJetEHist->Draw("COLZ");
1310   recoVsGenChargedJetEHist->SetTitle("Reconstructed Vs Generator Jet Energy;Gen E;Reco E");
1311   f1_34->Draw("SAME");
1312   f2_34->Draw("SAME");
1313   f3_34->Draw("SAME");
1314   gPad->SetLogz();
1315   if(PRINT) c34->Print((results_path+"/matchedRecoVsGenJetEnergy.png").c_str()); // Matched reconstructed vs generator jet energy
1316
1317   // Jet Res Vs Gen Eta
1318   TCanvas *c35 = new TCanvas("c35","Jet Res Vs Gen Eta",800,600);
1319   c35->Clear();
1320   c35->Divide(1,1);
1321
1322   c35->cd(1);
1323   jetResVsEtaHist->Draw("COLZ");
1324   jetResVsEtaHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Eta;Gen Eta;Res");
1325   gPad->SetLogz();
1326   if(PRINT) c35->Print((results_path+"/matchedJetResolutionVsEta.png").c_str()); // Matched jet resolution vs generator jet eta
1327
1328   // Jet Res Vs Gen E
1329   TCanvas *c36 = new TCanvas("c36","Jet Res Vs Gen E",800,600);
1330   c36->Clear();
1331   c36->Divide(1,1);
1332
1333   c36->cd(1);
1334   jetResVsEHist->Draw("COLZ");
1335   jetResVsEHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy;Gen E;Res");
1336   gPad->SetLogz();
1337   if(PRINT) c36->Print((results_path+"/matchedJetResolutionVsEnergy.png").c_str()); // Matched jet resolution vs generator jet energy
1338
1339   // Jet Res Vs Gen E Neg Eta
1340   TCanvas *c37 = new TCanvas("c37","Jet Res Vs Gen E (-2.5 < eta < -1.0)",800,600);
1341   c37->Clear();
1342   c37->Divide(1,1);
1343
1344   c37->cd(1);
1345   jetResVsENegEtaNoDupHist->Draw("COLZ");
1346   jetResVsENegEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-2.5 < eta < -1.0) No Duplicate;Gen E;Res");
1347   gPad->SetLogz();
1348   if(PRINT) c37->Print((results_path+"/matchedJetResolutionVsEnergyNegEta.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0
1349
1350   // Jet Res Vs Gen E Mid Eta
1351   TCanvas *c38 = new TCanvas("c38","Jet Res Vs Gen E (-1.0 < eta < 1.0)",800,600);
1352   c38->Clear();
1353   c38->Divide(1,1);
1354
1355   c38->cd(1);
1356   jetResVsEMidEtaNoDupHist->Draw("COLZ");
1357   jetResVsEMidEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-1.0 < eta < 1.0) No Duplicate;Gen E;Res");
1358   gPad->SetLogz();
1359   if(PRINT) c38->Print((results_path+"/matchedJetResolutionVsEnergyMidEta.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0
1360     delete c38;
1361   // Jet Res Vs Gen E Pos Eta
1362   TCanvas *c39 = new TCanvas("c39","Jet Res Vs Gen E (1.0 < eta < 2.5)",800,600);
1363   c39->Clear();
1364   c39->Divide(1,1);
1365
1366   c39->cd(1);
1367   jetResVsEPosEtaNoDupHist->Draw("COLZ");
1368   jetResVsEPosEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (1.0 < eta < 2.5) No Duplicate;Gen E;Res");
1369   gPad->SetLogz();
1370   if(PRINT) c39->Print((results_path+"/matchedJetResolutionVsEnergyPosEta.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5
1371   delete c39;
1372
1373
1374   // Generate Resolution Plots
1375   const int BINS = 20;
1376   double binCent[BINS];
1377   double jesVsENeg[BINS];
1378   double jesVsEMid[BINS];
1379   double jesVsEPos[BINS];
1380   double jerVsENeg[BINS];
1381   double jerVsEMid[BINS];
1382   double jerVsEPos[BINS];
1383
1384   std::fill(std::begin(binCent), std::end(binCent), -999.);
1385   std::fill(std::begin(jesVsENeg), std::end(jesVsENeg), -999.);
1386   std::fill(std::begin(jesVsEMid), std::end(jesVsEMid), -999.);
1387   std::fill(std::begin(jesVsEPos), std::end(jesVsEPos), -999.);
1388   std::fill(std::begin(jerVsENeg), std::end(jerVsENeg), -999.);
1389   std::fill(std::begin(jerVsEMid), std::end(jerVsEMid), -999.);
1390   std::fill(std::begin(jerVsEPos), std::end(jerVsEPos), -999.);
1391
1392   TH1D *pxA = jetResVsENegEtaNoDupHist->ProjectionX("pxA",1,10000);
1393   for(int i=0; i<BINS; i++)
1394     {
1395       binCent[i] = pxA->GetBinCenter(i+1);
1396     }
1397
1398   TCanvas *c40 = new TCanvas("c40","Negative Rapidity Fit Results",800,600);
1399   c40->Clear();
1400   c40->Divide(5,4);
1401
1402   TH1D *hA[20];
1403   for(int i=1; i<21; i++)
1404     {
1405       hA[i-1] = (TH1D *)jetResVsENegEtaNoDupHist->ProjectionY(Form("projYA_%d",i),i,i);
1406
1407       TF1 *myGausA = new TF1("myGausA","gaus",-0.5,0.5);
1408       myGausA->SetParameters(hA[i-1]->GetMaximum(),0.0,0.01);
1409
1410       c40->cd(i);
1411       //hA[i-1]->Draw("HIST");
1412       hA[i-1]->Fit("myGausA","B","",-0.5,0.5);
1413       hA[i-1]->GetXaxis()->SetRangeUser(-1,1);
1414       gPad->SetLogy();
1415
1416       if(hA[i-1]->GetEntries() > 2)
1417     {
1418       auto fA = hA[i-1]->GetFunction("myGausA");
1419
1420       if(fA->GetParError(2)/fA->GetParameter(2) < 0.5)
1421         {
1422           jesVsENeg[i-1] = fA->GetParameter(1);
1423           jerVsENeg[i-1] = fA->GetParameter(2);
1424         }
1425       //cout << fA->GetParameter(0) << " " << fA->GetParameter(1) << " " << fA->GetParameter(2) << endl;
1426       //cout << fA->GetParError(0) << " " << fA->GetParError(1) << " " << fA->GetParError(2) << endl;
1427     }
1428     }
1429   if(PRINT) c40->Print((results_path+"/matchedJetResolutionVsEnergyNegEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0 fits
1430     delete c40;
1431
1432   TCanvas *c41 = new TCanvas("c41","Mid Rapidity Fit Results",800,600);
1433   c41->Clear();
1434   c41->Divide(5,4);
1435
1436   TH1D *hB[20];
1437   for(int i=1; i<21; i++)
1438     {
1439       hB[i-1] = (TH1D *)jetResVsEMidEtaNoDupHist->ProjectionY(Form("projYB_%d",i),i,i);
1440
1441       TF1 *myGausB = new TF1("myGausB","gaus",-0.5,0.5);
1442       myGausB->SetParameters(hB[i-1]->GetMaximum(),0.0,0.01);
1443
1444       c41->cd(i);
1445       //hA[i-1]->Draw("HIST");
1446       hB[i-1]->Fit("myGausB","B","",-0.5,0.5);
1447       hB[i-1]->GetXaxis()->SetRangeUser(-1,1);
1448       gPad->SetLogy();
1449
1450       if(hB[i-1]->GetEntries() > 2)
1451     {
1452       auto fB = hB[i-1]->GetFunction("myGausB");
1453
1454       if(fB->GetParError(2)/fB->GetParameter(2) < 0.5)
1455         {
1456           jesVsEMid[i-1] = fB->GetParameter(1);
1457           jerVsEMid[i-1] = fB->GetParameter(2);
1458         }
1459       //cout << fB->GetParameter(0) << " " << fB->GetParameter(1) << " " << fB->GetParameter(2) << endl;
1460       //cout << fB->GetParError(0) << " " << fB->GetParError(1) << " " << fB->GetParError(2) << endl;
1461     }
1462     }
1463   if(PRINT) c41->Print((results_path+"/matchedJetResolutionVsEnergyMidEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0 fits
1464     delete c41;
1465
1466   TCanvas *c42 = new TCanvas("c42","Positive Rapidity Fit Results",800,600);
1467   c42->Clear();
1468   c42->Divide(5,4);
1469
1470   TH1D *hC[20];
1471   for(int i=1; i<21; i++)
1472     {
1473       hC[i-1] = (TH1D *)jetResVsEPosEtaNoDupHist->ProjectionY(Form("projYC_%d",i),i,i);
1474
1475       TF1 *myGausC = new TF1("myGausC","gaus",-0.5,0.5);
1476       myGausC->SetParameters(hC[i-1]->GetMaximum(),0.0,0.01);
1477
1478       c42->cd(i);
1479       //hA[i-1]->Draw("HIST");
1480       hC[i-1]->Fit("myGausC","B","",-0.5,0.5);
1481       hC[i-1]->GetXaxis()->SetRangeUser(-1,1);
1482       gPad->SetLogy();
1483
1484       if(hC[i-1]->GetEntries() > 2)
1485     {
1486       auto fC = hC[i-1]->GetFunction("myGausC");
1487
1488       if(fC->GetParError(2)/fC->GetParameter(2) < 0.5)
1489         {
1490           jesVsEPos[i-1] = fC->GetParameter(1);
1491           jerVsEPos[i-1] = fC->GetParameter(2);
1492         }
1493       //cout << fC->GetParameter(0) << " " << fC->GetParameter(1) << " " << fC->GetParameter(2) << endl;
1494       //cout << fC->GetParError(0) << " " << fC->GetParError(1) << " " << fC->GetParError(2) << endl;
1495     }
1496     }
1497   if(PRINT) c42->Print((results_path+"/matchedJetResolutionVsEnergyPosEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5 fits
1498     delete c42;
1499   TCanvas *c43 = new TCanvas("c43","Positive JES/JER",800,600);
1500   c43->Clear();
1501   c43->Divide(1,1);
1502
1503   TGraph *gJESvsENeg = new TGraph(BINS,binCent,jesVsENeg);
1504   TGraph *gJERvsENeg = new TGraph(BINS,binCent,jerVsENeg);
1505
1506   TGraph *gJESvsEMid = new TGraph(BINS,binCent,jesVsEMid);
1507   TGraph *gJERvsEMid = new TGraph(BINS,binCent,jerVsEMid);
1508
1509   TGraph *gJESvsEPos = new TGraph(BINS,binCent,jesVsEPos);
1510   TGraph *gJERvsEPos = new TGraph(BINS,binCent,jerVsEPos);
1511
1512   TH2D *test43 = new TH2D("test43","Jet Energy Scale / Resolution Vs Eta;True Eta;JES/JER",1,0.,100.,1,-0.2,0.2);
1513   test43->Draw();
1514
1515   c43->cd(1);
1516   gJERvsENeg->Draw("*");
1517   gJERvsENeg->SetMarkerStyle(21);
1518   gJERvsENeg->SetMarkerSize(1);
1519   gJERvsENeg->SetMarkerColor(seabornBlue);
1520
1521   gJESvsENeg->Draw("*");
1522   gJESvsENeg->SetMarkerStyle(26);
1523   gJESvsENeg->SetMarkerSize(1);
1524   gJESvsENeg->SetMarkerColor(seabornBlue);
1525
1526   gJERvsEMid->Draw("*");
1527   gJERvsEMid->SetMarkerStyle(21);
1528   gJERvsEMid->SetMarkerSize(1);
1529   gJERvsEMid->SetMarkerColor(seabornRed);
1530
1531   gJESvsEMid->Draw("*");
1532   gJESvsEMid->SetMarkerStyle(26);
1533   gJESvsEMid->SetMarkerSize(1);
1534   gJESvsEMid->SetMarkerColor(seabornRed);
1535
1536   gJERvsEPos->Draw("*");
1537   gJERvsEPos->SetMarkerStyle(21);
1538   gJERvsEPos->SetMarkerSize(1);
1539   gJERvsEPos->SetMarkerColor(seabornGreen);
1540
1541   gJESvsEPos->Draw("*");
1542   gJESvsEPos->SetMarkerStyle(26);
1543   gJESvsEPos->SetMarkerSize(1);
1544   gJESvsEPos->SetMarkerColor(seabornGreen);
1545
1546   TLegend *legend = new TLegend(0.7,0.7,0.9,0.9);
1547 legend->AddEntry(gJERvsENeg, "JER, (-2.5 < #eta < -1)", "p");
1548 legend->AddEntry(gJESvsENeg, "JES, (-2.5 < #eta < -1)","p");
1549 legend->AddEntry(gJERvsEMid, "JER, (-1 < #eta < 1)", "p");
1550 legend->AddEntry(gJESvsEMid, "JES, (-1 < #eta < 1)", "p");
1551 legend->AddEntry(gJERvsEPos, "JER, (1 < #eta < 2.5) ", "p");
1552 legend->AddEntry(gJESvsEPos, "JES, (1 < #eta < 2.5)", "p");
1553 legend->Draw();
1554
1555   if(PRINT) c43->Print((results_path+"/matchedJetScaleResolutionSummary.png").c_str()); // Matched jet JER/JES summary
1556     delete c43;
1557
1558 delete mychain;
1559
1560   return 0;
1561 }