Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 10:18:28 

0001 #include "common_bench/benchmark.h"
0002 #include "common_bench/mt.h"
0003 #include "common_bench/util.h"
0004 #include "common_bench/plot.h"
0005 
0006 #include <TCanvas.h>
0007 #include <TChain.h>
0008 #include <TFile.h>
0009 #include <TGraph.h>
0010 #include <TH1D.h>
0011 #include <TH2D.h>
0012 #include <TTreeReader.h>
0013 #include <TTreeReaderArray.h>
0014 #include <TLegend.h>
0015 #include <TVector3.h>
0016 
0017 #include "fmt/color.h"
0018 #include "fmt/core.h"
0019 
0020 #include "nlohmann/json.hpp"
0021 
0022 // Jet Benchmarks
0023 // Author: B. Page (bpage@bnl.gov)
0024 
0025 // To run: snakemake -c1 results/epic_craterlake/jets/10on100/minQ2=1
0026 
0027 int jets(const std::string& config_name)
0028 {
0029   // read our configuration
0030   std::ifstream  config_file{config_name};
0031   nlohmann::json config;
0032   config_file >> config;
0033 
0034   const std::string rec_file     = config["rec_file"];
0035   const std::string detector     = config["detector"];
0036   const std::string results_path = config["results_path"];
0037   const std::string plot_tag     = config["plot_tag"];
0038   const std::string test_tag     = config["test_tag"];
0039   const int ebeam                = config["ebeam"];
0040   const int pbeam                = config["pbeam"];
0041 
0042   fmt::print(fmt::emphasis::bold | fg(fmt::color::forest_green),
0043              "Running jet analysis...\n");
0044   fmt::print(" - Detector package: {}\n", detector);
0045   fmt::print(" - input file: {}\n", rec_file);
0046   fmt::print(" - test tag: {}\n", test_tag);
0047   fmt::print(" - ebeam: {}\n", ebeam);
0048   fmt::print(" - pbeam: {}\n", pbeam);
0049 
0050   const bool PRINT = true;
0051 
0052   // Input
0053   TChain *mychain = new TChain("events");
0054   mychain->Add(rec_file.c_str());
0055 
0056   const int seabornRed = TColor::GetColor(213, 94, 0);
0057 
0058 // Seaborn Green: #009E73 -> (0, 158, 115)
0059 const int seabornGreen = TColor::GetColor(0, 158, 115);
0060 
0061 // Seaborn Blue: #56B4E9 -> (86, 180, 233)
0062 const int seabornBlue = TColor::GetColor(100, 149, 237);
0063 
0064   // Output
0065   //TFile *ofile = TFile::Open("test_24-05-0.hist.root","RECREATE");
0066 
0067   // TTreeReader
0068   TTreeReader tree_reader(mychain);
0069 
0070   // Set Delta R Cut
0071   float DELTARCUT = 0.05;
0072 
0073   // Reco Jets
0074   TTreeReaderArray<int> recoType = {tree_reader, "ReconstructedChargedJets.type"};
0075   TTreeReaderArray<float> recoNRG = {tree_reader, "ReconstructedChargedJets.energy"};
0076   TTreeReaderArray<int> recoPDG = {tree_reader, "ReconstructedChargedJets.PDG"};
0077   TTreeReaderArray<float> recoMomX = {tree_reader, "ReconstructedChargedJets.momentum.x"};
0078   TTreeReaderArray<float> recoMomY = {tree_reader, "ReconstructedChargedJets.momentum.y"};
0079   TTreeReaderArray<float> recoMomZ = {tree_reader, "ReconstructedChargedJets.momentum.z"};
0080   TTreeReaderArray<float> recoM = {tree_reader, "ReconstructedChargedJets.mass"};
0081   TTreeReaderArray<unsigned int> partsBegin = {tree_reader, "ReconstructedChargedJets.particles_begin"};
0082   TTreeReaderArray<unsigned int> partsEnd = {tree_reader, "ReconstructedChargedJets.particles_end"};
0083 
0084   TTreeReaderArray<int> recoClustIndex = {tree_reader, "_ReconstructedChargedJets_clusters.index"};
0085   TTreeReaderArray<int> recoTrackIndex = {tree_reader, "_ReconstructedChargedJets_tracks.index"};
0086   TTreeReaderArray<int> recoPartIndex = {tree_reader, "_ReconstructedChargedJets_particles.index"};
0087 
0088   // Reconstructed Particles
0089   TTreeReaderArray<float> recoPartMomX = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
0090   TTreeReaderArray<float> recoPartMomY = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
0091   TTreeReaderArray<float> recoPartMomZ = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
0092   TTreeReaderArray<float> recoPartM = {tree_reader, "ReconstructedChargedParticles.mass"};
0093   TTreeReaderArray<int> recoPartPDG = {tree_reader, "ReconstructedChargedParticles.PDG"};
0094   TTreeReaderArray<float> recoPartNRG = {tree_reader, "ReconstructedChargedParticles.energy"};
0095 
0096   // Generated Jets
0097   TTreeReaderArray<int> genType = {tree_reader, "GeneratedChargedJets.type"};
0098   TTreeReaderArray<float> genNRG = {tree_reader, "GeneratedChargedJets.energy"};
0099   TTreeReaderArray<int> genPDG = {tree_reader, "GeneratedChargedJets.PDG"};
0100   TTreeReaderArray<float> genMomX = {tree_reader, "GeneratedChargedJets.momentum.x"};
0101   TTreeReaderArray<float> genMomY = {tree_reader, "GeneratedChargedJets.momentum.y"};
0102   TTreeReaderArray<float> genMomZ = {tree_reader, "GeneratedChargedJets.momentum.z"};
0103   TTreeReaderArray<float> genM = {tree_reader, "GeneratedChargedJets.mass"};
0104   TTreeReaderArray<unsigned int> genPartsBegin = {tree_reader, "GeneratedChargedJets.particles_begin"};
0105   TTreeReaderArray<unsigned int> genPartsEnd = {tree_reader, "GeneratedChargedJets.particles_end"};
0106   
0107   TTreeReaderArray<int> genPartIndex = {tree_reader, "_GeneratedChargedJets_particles.index"};
0108   //TTreeReaderArray<int> genChargedIndex = {tree_reader, "GeneratedChargedParticles_objIdx.index"};
0109   
0110   // MC
0111   //TTreeReaderArray<int> mcGenStat = {tree_reader, "MCParticles.generatorStatus"};
0112   TTreeReaderArray<float> mcMomX = {tree_reader, "GeneratedParticles.momentum.x"};
0113   TTreeReaderArray<float> mcMomY = {tree_reader, "GeneratedParticles.momentum.y"};
0114   TTreeReaderArray<float> mcMomZ = {tree_reader, "GeneratedParticles.momentum.z"};
0115   TTreeReaderArray<float> mcM = {tree_reader, "GeneratedParticles.mass"};
0116   TTreeReaderArray<int> pdg = {tree_reader, "GeneratedParticles.PDG"};
0117 
0118   // Define Histograms
0119   TH1D *counter = new TH1D("counter","",10,0.,10.);
0120 
0121   
0122   // Reco
0123   TH1D *numRecoChargedJetsECutHist = new TH1D("numRecoChargedJetsECut","",20,0.,20.);
0124   TH1D *recoChargedJetEHist = new TH1D("recoChargedJetE","",300,0.,300.);
0125   TH1D *recoChargedJetEtaECutHist = new TH1D("recoChargedJetEtaECut","",60,-3.,3.);
0126   TH2D *recoChargedJetEvsEtaHist = new TH2D("recoChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0127   TH2D *recoChargedJetPhiVsEtaECutHist = new TH2D("recoChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0128 
0129   TH1D *numRecoChargedJetsECutNoElecHist = new TH1D("numRecoChargedJetsECutNoElec","",20,0.,20.);
0130   TH1D *recoChargedJetENoElecHist = new TH1D("recoChargedJetENoElec","",300,0.,300.);
0131   TH1D *recoChargedJetEtaECutNoElecHist = new TH1D("recoChargedJetEtaECutNoElec","",60,-3.,3.);
0132   TH2D *recoChargedJetEvsEtaNoElecHist = new TH2D("recoChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0133   TH2D *recoChargedJetPhiVsEtaECutNoElecHist = new TH2D("recoChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0134 
0135   TH1D *numRecoChargedJetPartsHist = new TH1D("numRecoChargedJetParts","",20,0.,20.);
0136   TH1D *recoChargedJetPartEHist = new TH1D("recoChargedJetPartE","",500,0.,100.);
0137   TH1D *recoChargedJetPartEtaHist = new TH1D("recoChargedJetPartEta","",80,-4.,4.);
0138   TH2D *recoChargedJetPartEvsEtaHist = new TH2D("recoChargedJetPartEvsEta","",80,-4.,4.,500,0.,100.);
0139   TH2D *recoChargedJetPartPhiVsEtaHist = new TH2D("recoChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0140 
0141   TH1D *numRecoChargedJetPartsNoElecHist = new TH1D("numRecoChargedJetPartsNoElec","",20,0.,20.);
0142   TH1D *recoChargedJetPartENoElecHist = new TH1D("recoChargedJetPartENoElec","",500,0.,100.);
0143   TH1D *recoChargedJetPartEtaNoElecHist = new TH1D("recoChargedJetPartEtaNoElec","",80,-4.,4.);
0144   TH2D *recoChargedJetPartEvsEtaNoElecHist = new TH2D("recoChargedJetPartEvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0145   TH2D *recoChargedJetPartPhiVsEtaNoElecHist = new TH2D("recoChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0146 
0147   TH1D *recoChargedJetPartPairwiseDeltaRHist = new TH1D("recoChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0148 
0149   // Gen
0150   TH1D *numGenChargedJetsECutHist = new TH1D("numGenChargedJetsECut","",20,0.,20.);
0151   TH1D *genChargedJetEHist = new TH1D("genChargedJetE","",300,0.,300.);
0152   TH1D *genChargedJetEtaECutHist = new TH1D("genChargedJetEtaECut","",60,-3.,3.);
0153   TH2D *genChargedJetEvsEtaHist = new TH2D("genChargedJetEvsEta","",60,-3.,3.,300,0.,300.);
0154   TH2D *genChargedJetPhiVsEtaECutHist = new TH2D("genChargedJetPhiVsEtaECut","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0155 
0156   TH1D *numGenChargedJetsECutNoElecHist = new TH1D("numGenChargedJetsECutNoElec","",20,0.,20.);
0157   TH1D *genChargedJetENoElecHist = new TH1D("genChargedJetENoElec","",300,0.,300.);
0158   TH1D *genChargedJetEtaECutNoElecHist = new TH1D("genChargedJetEtaECutNoElec","",60,-3.,3.);
0159   TH2D *genChargedJetEvsEtaNoElecHist = new TH2D("genChargedJetEvsEtaNoElec","",60,-3.,3.,300,0.,300.);
0160   TH2D *genChargedJetPhiVsEtaECutNoElecHist = new TH2D("genChargedJetPhiVsEtaECutNoElec","",60,-3.,3.,100,-TMath::Pi(),TMath::Pi());
0161 
0162   TH1D *numGenChargedJetPartsHist = new TH1D("numGenChargedJetParts","",20,0.,20.);
0163   TH1D *genChargedJetPartEHist = new TH1D("genChargedJetPartE","",500,0.,100.);
0164   TH1D *genChargedJetPartEtaHist = new TH1D("genChargedJetPartEta","",80,-4.,4.);
0165   TH2D *genChargedJetPartEvsEtaHist = new TH2D("genChargedJetPartEvsEta","",80,-4.,4.,500,0.,100.);
0166   TH2D *genChargedJetPartPhiVsEtaHist = new TH2D("genChargedJetPartPhiVsEta","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0167 
0168   TH1D *numGenChargedJetPartsNoElecHist = new TH1D("numGenChargedJetPartsNoElec","",20,0.,20.);
0169   TH1D *genChargedJetPartENoElecHist = new TH1D("genChargedJetPartENoElec","",500,0.,100.);
0170   TH1D *genChargedJetPartEtaNoElecHist = new TH1D("genChargedJetPartEtaNoElec","",80,-4.,4.);
0171   TH2D *genChargedJetPartEvsEtaNoElecHist = new TH2D("genChargedJetPartEvsEtaNoElec","",80,-4.,4.,500,0.,100.);
0172   TH2D *genChargedJetPartPhiVsEtaNoElecHist = new TH2D("genChargedJetPartPhiVsEtaNoElec","",80,-4.,4.,100,-TMath::Pi(),TMath::Pi());
0173 
0174   TH1D *genChargedJetPartPairwiseDeltaRHist = new TH1D("genChargedJetPartPairwiseDeltaRHist","",5000,0.,5.);
0175 
0176   // Matched
0177   TH1D *matchJetDeltaRHist = new TH1D("matchJetDeltaR","",5000,0.,5.);
0178   TH1D *matchJetDeltaRBackHist = new TH1D("matchJetDeltaRBack","",5000,0.,5.);
0179   TH2D *recoVsGenChargedJetEtaHist = new TH2D("recoVsGenChargedJetEta","",80,-4.,4.,80,-4.,4.);
0180   TH2D *recoVsGenChargedJetPhiHist = new TH2D("recoVsGenChargedJetPhi","",100,-TMath::Pi(),TMath::Pi(),100,-TMath::Pi(),TMath::Pi());
0181   TH2D *recoVsGenChargedJetEHist = new TH2D("recoVsGenChargedJetE","",100,0.,100.,100,0.,100.);
0182   TH2D *recoVsGenChargedJetENoDRHist = new TH2D("recoVsGenChargedJetENoDRHist","",100,0.,100.,100,0.,100.);
0183   TH2D *recoVsGenChargedJetENoDupHist = new TH2D("recoVsGenChargedJetENoDup","",100,0.,100.,100,0.,100.);
0184 
0185   TH2D *jetResVsEtaHist = new TH2D("jetResVsEta","",80,-4.,4.,10000,-10.,10.);
0186   TH2D *jetResVsEHist = new TH2D("jetResVsE","",100,0.,100.,10000,-10.,10.);
0187   TH2D *jetResVsENegEtaHist = new TH2D("jetResVsENegEta","",20,0.,100.,10000,-10.,10.);
0188   TH2D *jetResVsEMidEtaHist = new TH2D("jetResVsEMidEta","",20,0.,100.,10000,-10.,10.);
0189   TH2D *jetResVsEPosEtaHist = new TH2D("jetResVsEPosEta","",20,0.,100.,10000,-10.,10.);
0190 
0191   TH2D *jetResVsENegEtaNoDupHist = new TH2D("jetResVsENegEtaNoDup","",20,0.,100.,10000,-10.,10.);
0192   TH2D *jetResVsEMidEtaNoDupHist = new TH2D("jetResVsEMidEtaNoDup","",20,0.,100.,10000,-10.,10.);
0193   TH2D *jetResVsEPosEtaNoDupHist = new TH2D("jetResVsEPosEtaNoDup","",20,0.,100.,10000,-10.,10.);
0194 
0195 
0196   // Loop Through Events
0197   int NEVENTS = 0;
0198   while(tree_reader.Next()) {
0199 
0200     if(NEVENTS%10000 == 0) cout << "Events Processed: " << NEVENTS << endl;
0201 
0202     counter->Fill(0);
0203 
0204     //////////////////////////////////////////////////////////////////////////
0205     //////////////////////  Analyze Reconstructed Jets  //////////////////////
0206     //////////////////////////////////////////////////////////////////////////
0207     int numRecoChargedJets = 0;
0208     int numRecoChargedJetsNoElec = 0;
0209     for(unsigned int i=0; i<recoType.GetSize(); i++)
0210       {
0211     TVector3 jetMom(recoMomX[i],recoMomY[i],recoMomZ[i]);
0212 
0213     counter->Fill(3);
0214 
0215     // Place eta cut to avoid edges of tracking acceptance
0216     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0217 
0218     // Place a minimum energy condition for several plots
0219     bool ECut = recoNRG[i] > 5.0;
0220 
0221     if(ECut) numRecoChargedJets++; 
0222 
0223     recoChargedJetEHist->Fill(recoNRG[i]);
0224     if(ECut) recoChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0225     recoChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0226     if(ECut) recoChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0227 
0228     bool noElectron = true;
0229     if(ECut)
0230       {
0231         // Find Jets with Electrons
0232         for(unsigned int m=partsBegin[i]; m<partsEnd[i]; m++)
0233           {
0234         if(recoPartM[recoPartIndex[m]] > 0.00050 && recoPartM[recoPartIndex[m]] < 0.00052)
0235           noElectron = false;
0236           }
0237 
0238         for(unsigned int j=partsBegin[i]; j<partsEnd[i]; j++)
0239           {
0240         // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
0241         // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
0242         double mX = recoPartMomX[recoPartIndex[j]];
0243         double mY = recoPartMomY[recoPartIndex[j]];
0244         double mZ = recoPartMomZ[recoPartIndex[j]];
0245         double mM = recoPartM[recoPartIndex[j]];
0246         double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0247         
0248         TVector3 partMom(mX,mY,mZ);
0249         
0250         recoChargedJetPartEHist->Fill(tmpE);
0251         recoChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0252         recoChargedJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
0253         recoChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0254 
0255         if(noElectron)
0256           {
0257             recoChargedJetPartENoElecHist->Fill(tmpE);
0258             recoChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0259             recoChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
0260             recoChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0261           }
0262 
0263         // Pairwise Distance Between Constituents
0264         if(j<(partsEnd[i]-1))
0265           {
0266             for(unsigned int k=j+1; k<partsEnd[i]; k++)
0267               {
0268             double mXB = recoPartMomX[recoPartIndex[k]];
0269             double mYB = recoPartMomY[recoPartIndex[k]];
0270             double mZB = recoPartMomZ[recoPartIndex[k]];
0271 
0272             TVector3 partMomB(mXB,mYB,mZB);
0273 
0274             double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0275             double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0276             double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0277 
0278             recoChargedJetPartPairwiseDeltaRHist->Fill(dR);
0279               }
0280           }
0281           }
0282         numRecoChargedJetPartsHist->Fill(partsEnd[i] - partsBegin[i]);
0283         if(noElectron) numRecoChargedJetPartsNoElecHist->Fill(partsEnd[i] - partsBegin[i]);
0284       }
0285 
0286     // No Electrons
0287     if(noElectron)
0288       {
0289         recoChargedJetENoElecHist->Fill(recoNRG[i]);
0290         if(ECut) recoChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0291         recoChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),recoNRG[i]);
0292         if(ECut) recoChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0293 
0294         if(ECut) numRecoChargedJetsNoElec++; 
0295       }
0296       }
0297     numRecoChargedJetsECutHist->Fill(numRecoChargedJets);
0298     numRecoChargedJetsECutNoElecHist->Fill(numRecoChargedJetsNoElec);
0299 
0300     //////////////////////////////////////////////////////////////////////////
0301     ////////////////////////  Analyze Generator Jets  ////////////////////////
0302     //////////////////////////////////////////////////////////////////////////
0303     int numGenChargedJets = 0;
0304     int numGenChargedJetsNoElec = 0;
0305     for(unsigned int i=0; i<genType.GetSize(); i++)
0306       {
0307     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0308 
0309     counter->Fill(4);
0310 
0311     // Place eta cut to avoid edges of tracking acceptance
0312     if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0313 
0314     // Place a minimum energy condition for several plots
0315     bool ECut = genNRG[i] > 5.0;
0316 
0317     if(ECut) numGenChargedJets++; 
0318 
0319     genChargedJetEHist->Fill(genNRG[i]);
0320     if(ECut) genChargedJetEtaECutHist->Fill(jetMom.PseudoRapidity());
0321     genChargedJetEvsEtaHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0322     if(ECut) genChargedJetPhiVsEtaECutHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0323 
0324     bool noElectron = true;
0325     if(ECut)
0326       {
0327         // Find Jets with Electrons
0328         for(unsigned int m=genPartsBegin[i]; m<genPartsEnd[i]; m++)
0329           {
0330         if(mcM[genPartIndex[m]] > 0.00050 && mcM[genPartIndex[m]] < 0.00052)
0331           noElectron = false;
0332           }
0333 
0334         for(unsigned int j=genPartsBegin[i]; j<genPartsEnd[i]; j++)
0335           {
0336         // partsbegin and partsEnd specify the entries from _ReconstructedChargedJets_particles.index that make up the jet
0337         // _ReconstructedChargedJets_particles.index stores the ReconstructedChargedParticles index of the jet constituent
0338         double mX = mcMomX[genPartIndex[j]];
0339         double mY = mcMomY[genPartIndex[j]];
0340         double mZ = mcMomZ[genPartIndex[j]];
0341         double mM = mcM[genPartIndex[j]];
0342         double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0343         
0344         TVector3 partMom(mX,mY,mZ);
0345         
0346         genChargedJetPartEHist->Fill(tmpE);
0347         genChargedJetPartEtaHist->Fill(partMom.PseudoRapidity());
0348         genChargedJetPartEvsEtaHist->Fill(partMom.PseudoRapidity(),tmpE);
0349         genChargedJetPartPhiVsEtaHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0350 
0351         if(noElectron)
0352           {
0353             genChargedJetPartENoElecHist->Fill(tmpE);
0354             genChargedJetPartEtaNoElecHist->Fill(partMom.PseudoRapidity());
0355             genChargedJetPartEvsEtaNoElecHist->Fill(partMom.PseudoRapidity(),tmpE);
0356             genChargedJetPartPhiVsEtaNoElecHist->Fill(partMom.PseudoRapidity(),partMom.Phi());
0357           }
0358 
0359         // Pairwise Distance Between Constituents
0360         if(j<(genPartsEnd[i]-1))
0361           {
0362             for(unsigned int k=j+1; k<genPartsEnd[i]; k++)
0363               {
0364             double mXB = mcMomX[genPartIndex[k]];
0365             double mYB = mcMomY[genPartIndex[k]];
0366             double mZB = mcMomZ[genPartIndex[k]];
0367 
0368             TVector3 partMomB(mXB,mYB,mZB);
0369 
0370             double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0371             double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0372             double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0373 
0374             genChargedJetPartPairwiseDeltaRHist->Fill(dR);
0375               }
0376           }
0377           }
0378         numGenChargedJetPartsHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
0379         if(noElectron) numGenChargedJetPartsNoElecHist->Fill(genPartsEnd[i] - genPartsBegin[i]);
0380       }
0381 
0382     // No Electrons
0383     if(noElectron)
0384       {
0385         genChargedJetENoElecHist->Fill(genNRG[i]);
0386         if(ECut) genChargedJetEtaECutNoElecHist->Fill(jetMom.PseudoRapidity());
0387         genChargedJetEvsEtaNoElecHist->Fill(jetMom.PseudoRapidity(),genNRG[i]);
0388         if(ECut) genChargedJetPhiVsEtaECutNoElecHist->Fill(jetMom.PseudoRapidity(),jetMom.Phi());
0389 
0390         if(ECut) numGenChargedJetsNoElec++; 
0391       }
0392       }
0393     numGenChargedJetsECutHist->Fill(numGenChargedJets);
0394     numGenChargedJetsECutNoElecHist->Fill(numGenChargedJetsNoElec);
0395 
0396     
0397     //////////////////////////////////////////////////////////////////////////
0398     /////////////////////////////  Matched Jets  /////////////////////////////
0399     //////////////////////////////////////////////////////////////////////////
0400     for(unsigned int i=0; i<genType.GetSize(); i++)
0401       {
0402     TVector3 jetMom(genMomX[i],genMomY[i],genMomZ[i]);
0403 
0404     // Place eta cut to avoid edges of tracking acceptance
0405     //if(TMath::Abs(jetMom.PseudoRapidity()) > 2.5) continue;
0406 
0407     // Place a minimum energy condition
0408     //if(genNRG[i] < 5.0) continue;
0409     
0410     // Don't Look at Electron Jets
0411     bool hasElectron = false;
0412     // Find Jets with Electrons
0413     for(unsigned int m=genPartsBegin[i]; m<genPartsEnd[i]; m++)
0414       {
0415         if(mcM[genPartIndex[m]] > 0.00050 && mcM[genPartIndex[m]] < 0.00052)
0416           hasElectron = true;
0417       }
0418     //if(hasElectron) continue;
0419 
0420     // Find Matching Reconstructed Jet
0421     double minDeltaR = 999.;
0422     int minIndex = -1;
0423     for(unsigned int j=0; j<recoType.GetSize(); j++)
0424       {
0425         TVector3 recoMom(recoMomX[j],recoMomY[j],recoMomZ[j]);
0426 
0427         double dEta = jetMom.PseudoRapidity() - recoMom.PseudoRapidity();
0428         double dPhi = TVector2::Phi_mpi_pi(jetMom.Phi() - recoMom.Phi());
0429         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0430 
0431         if(dR < minDeltaR)
0432           {
0433         minDeltaR = dR;
0434         minIndex = j;
0435           }
0436       }
0437 
0438     // Do Backwards Match
0439     double minDeltaRBack = 999.;
0440     double minIndexBack = -1;
0441     if(minIndex > -1)
0442       {
0443         TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0444         for(unsigned int j=0; j<genType.GetSize(); j++)
0445           {
0446         TVector3 genMom(genMomX[j],genMomY[j],genMomZ[j]);
0447         
0448         double dEta = recoMatchMom.PseudoRapidity() - genMom.PseudoRapidity();
0449         double dPhi = TVector2::Phi_mpi_pi(recoMatchMom.Phi() - genMom.Phi());
0450         double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0451         
0452         if(dR < minDeltaRBack)
0453           {
0454             minDeltaRBack = dR;
0455             minIndexBack = j;
0456           }
0457           }
0458       }
0459 
0460     // Look at Best Match
0461     if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1 && !hasElectron) matchJetDeltaRHist->Fill(minDeltaR);
0462     if(genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && minIndex > -1) matchJetDeltaRBackHist->Fill(minDeltaR);
0463     if(minIndex > -1 && genNRG[i] > 5.0 && TMath::Abs(jetMom.PseudoRapidity()) < 2.5 && !hasElectron)
0464       {
0465         TVector3 recoMatchMom(recoMomX[minIndex],recoMomY[minIndex],recoMomZ[minIndex]);
0466 
0467         recoVsGenChargedJetENoDRHist->Fill(genNRG[i],recoNRG[minIndex]);
0468 
0469         if(minDeltaR < DELTARCUT)
0470           {
0471         recoVsGenChargedJetEtaHist->Fill(jetMom.PseudoRapidity(),recoMatchMom.PseudoRapidity());
0472         recoVsGenChargedJetPhiHist->Fill(jetMom.Phi(),recoMatchMom.Phi());
0473         recoVsGenChargedJetEHist->Fill(genNRG[i],recoNRG[minIndex]);
0474         
0475         double jetERes = (recoNRG[minIndex] - genNRG[i])/genNRG[i];
0476         
0477         jetResVsEtaHist->Fill(jetMom.PseudoRapidity(),jetERes);
0478         jetResVsEHist->Fill(genNRG[i],jetERes);
0479         if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0480           jetResVsENegEtaHist->Fill(genNRG[i],jetERes);
0481         if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0482           jetResVsEMidEtaHist->Fill(genNRG[i],jetERes);
0483         if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0484           jetResVsEPosEtaHist->Fill(genNRG[i],jetERes);
0485         
0486         // Check for Duplicate Tracks
0487         bool noDuplicate = true;
0488         for(unsigned int j=partsBegin[minIndex]; j<partsEnd[minIndex]; j++)
0489           {
0490             double mX = recoPartMomX[recoPartIndex[j]];
0491             double mY = recoPartMomY[recoPartIndex[j]];
0492             double mZ = recoPartMomZ[recoPartIndex[j]];
0493             double mM = recoPartM[recoPartIndex[j]];
0494             double tmpE = TMath::Sqrt(mX*mX + mY*mY + mZ*mZ + mM*mM);
0495             
0496             TVector3 partMom(mX,mY,mZ);
0497             
0498             // Pairwise Distance Between Constituents
0499             if(j<(partsEnd[minIndex]-1))
0500               {
0501             for(unsigned int k=j+1; k<partsEnd[minIndex]; k++)
0502               {
0503                 double mXB = recoPartMomX[recoPartIndex[k]];
0504                 double mYB = recoPartMomY[recoPartIndex[k]];
0505                 double mZB = recoPartMomZ[recoPartIndex[k]];
0506                 
0507                 TVector3 partMomB(mXB,mYB,mZB);
0508                 
0509                 double dEta = partMom.PseudoRapidity() - partMomB.PseudoRapidity();
0510                 double dPhi = TVector2::Phi_mpi_pi(partMom.Phi() - partMomB.Phi());
0511                 double dR = TMath::Sqrt(dEta*dEta + dPhi*dPhi);
0512                 
0513                 if(dR < 0.02) noDuplicate = false;
0514               }
0515               }
0516           }
0517 
0518         if(noDuplicate)
0519           {
0520             recoVsGenChargedJetENoDupHist->Fill(genNRG[i],recoNRG[minIndex]);
0521 
0522             if(jetMom.PseudoRapidity() > -2.5 && jetMom.PseudoRapidity() < -1.0)
0523               jetResVsENegEtaNoDupHist->Fill(genNRG[i],jetERes);
0524             if(jetMom.PseudoRapidity() > -1.0 && jetMom.PseudoRapidity() < 1.0)
0525               jetResVsEMidEtaNoDupHist->Fill(genNRG[i],jetERes);
0526             if(jetMom.PseudoRapidity() > 1.0 && jetMom.PseudoRapidity() < 2.5)
0527               jetResVsEPosEtaNoDupHist->Fill(genNRG[i],jetERes);
0528           }
0529           }
0530       }
0531       }
0532 
0533     NEVENTS++;
0534   }
0535 
0536   
0537   gStyle->SetOptStat(0);
0538   ////////////////////////  Reconstructed Jets Plots  ////////////////////////
0539   // Reco Number
0540   TCanvas *c1 = new TCanvas("c1","Number Reco Jets",800,600);
0541   c1->Clear();
0542   c1->Divide(1,1);
0543 
0544   c1->cd(1);
0545   numRecoChargedJetsECutHist->Draw("HIST");
0546   numRecoChargedJetsECutNoElecHist->SetLineColor(seabornRed);
0547   numRecoChargedJetsECutNoElecHist->Draw("HISTSAME");
0548   numRecoChargedJetsECutHist->SetLineWidth(2); // Set line width to 2 (adjust as needed)
0549 numRecoChargedJetsECutNoElecHist->SetLineWidth(2);
0550   numRecoChargedJetsECutHist->SetTitle("Reconstructed Jets per Event (|eta| < 2.5 && E > 5);Number");
0551 
0552 TLegend *legend1 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0553 legend1->AddEntry(numRecoChargedJetsECutHist, "With Electrons", "l");
0554 legend1->AddEntry(numRecoChargedJetsECutNoElecHist, "No Electrons", "l");
0555 legend1->Draw();
0556 
0557 
0558 
0559   gPad->SetLogy();
0560   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
0561    delete c1;
0562 
0563   // Reco Energy
0564   TCanvas *c2 = new TCanvas("c2","Reco Jet Energy",800,600);
0565   c2->Clear();
0566   c2->Divide(1,1);
0567 
0568   c2->cd(1);
0569   recoChargedJetEHist->Draw("HIST");
0570   recoChargedJetENoElecHist->SetLineColor(seabornRed);
0571   recoChargedJetENoElecHist->Draw("HISTSAME");
0572 
0573   recoChargedJetEHist->SetLineWidth(2);
0574   recoChargedJetENoElecHist->SetLineWidth(2);
0575   recoChargedJetEHist->SetTitle("Reconstructed Jet Energy (|eta| < 2.5);Energy [GeV]");
0576 
0577 TLegend *legend2 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0578 legend2->AddEntry(recoChargedJetEHist, "With Electrons", "l");
0579 legend2->AddEntry(recoChargedJetENoElecHist, "No Electrons", "l");
0580 legend2->Draw();
0581 
0582   gPad->SetLogy();
0583   if(PRINT) c2->Print((results_path+"/recoJetEnergy.png").c_str()); // Energy spectrum of reconstructed jets with Abs(eta) < 2.5
0584 
0585     delete c2;
0586   // Reco Eta
0587   TCanvas *c3 = new TCanvas("c3","Reco Jet Eta",800,600);
0588   c3->Clear();
0589   c3->Divide(1,1);
0590 
0591   c3->cd(1);
0592   recoChargedJetEtaECutHist->Draw("HIST");
0593   recoChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
0594   recoChargedJetEtaECutNoElecHist->Draw("HISTSAME");
0595 
0596   recoChargedJetEtaECutHist->SetLineWidth(2);
0597   recoChargedJetEtaECutNoElecHist->SetLineWidth(2);
0598   recoChargedJetEtaECutHist->SetTitle("Reconstructed Jet Eta (E > 5);Eta");
0599 
0600 //add legend 
0601 TLegend *legend3 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0602 legend3->AddEntry(recoChargedJetEtaECutHist, "With Electrons", "l");
0603 legend3->AddEntry(recoChargedJetEtaECutNoElecHist, "No Electrons", "l");
0604 legend3->Draw();
0605 
0606   gPad->SetLogy();
0607   if(PRINT) c3->Print((results_path+"/recoJetEta.png").c_str()); // Eta spectrum of reconstructed jets with energy > 5 GeV
0608     delete c3;
0609   // Reco E Vs Eta
0610   TCanvas *c4 = new TCanvas("c4","Reco Jet E Vs Eta",800,600);
0611   c4->Clear();
0612   c4->Divide(1,1);
0613 
0614   c4->cd(1);
0615   recoChargedJetEvsEtaHist->Draw("COLZ");
0616   recoChargedJetEvsEtaHist->SetTitle("Reconstructed Jet Energy Vs Eta;Eta;Energy [GeV]");
0617   gPad->SetLogz();
0618   if(PRINT) c4->Print((results_path+"/recoJetEnergyvsEta.png").c_str()); // Energy vs eta of reconstructed jets
0619 
0620   // Reco Phi Vs Eta
0621   TCanvas *c5 = new TCanvas("c5","Reco Jet Phi Vs Eta",800,600);
0622   c5->Clear();
0623   c5->Divide(1,1);
0624 
0625   c5->cd(1);
0626   recoChargedJetPhiVsEtaECutHist->Draw("COLZ");
0627   recoChargedJetPhiVsEtaECutHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5);Eta;Phi");
0628   gPad->SetLogz();
0629   if(PRINT) c5->Print((results_path+"/recoJetPhivsEta.png").c_str()); // Phi vs eta of reconstructed jets
0630 
0631   // Num Particles Per Reco Jet
0632   TCanvas *c6 = new TCanvas("c6","Number Constituents Per Reco Jet",800,600);
0633   c6->Clear();
0634   c6->Divide(1,1);
0635 
0636   c6->cd(1);
0637   numRecoChargedJetPartsHist->Draw("HIST");
0638   numRecoChargedJetPartsNoElecHist->SetLineColor(seabornRed);
0639   numRecoChargedJetPartsNoElecHist->Draw("HISTSAME");
0640 
0641   numRecoChargedJetPartsHist->SetLineWidth(2);
0642   numRecoChargedJetPartsNoElecHist->SetLineWidth(2);
0643   numRecoChargedJetPartsHist->SetTitle("Number of Constituents Per Reco Jet;Number of Constituents");
0644 
0645 TLegend *legend6 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0646 legend6->AddEntry(numRecoChargedJetPartsHist, "With Electrons", "l");
0647 legend6->AddEntry(numRecoChargedJetPartsNoElecHist, "No Electrons", "l");
0648 legend6->Draw();
0649 
0650   gPad->SetLogy();
0651   if(PRINT) c6->Print((results_path+"/numConstituentsPerRecoJet.png").c_str()); // Number of constituents in reconstructed jets
0652 
0653   // Reco Part Energy
0654   TCanvas *c7 = new TCanvas("c7","Reco Jet Constituent Energy",800,600);
0655   c7->Clear();
0656   c7->Divide(1,1);
0657 
0658   c7->cd(1);
0659   recoChargedJetPartEHist->Draw("HIST");
0660   recoChargedJetPartENoElecHist->SetLineColor(seabornRed);
0661   recoChargedJetPartENoElecHist->Draw("HISTSAME");
0662 
0663   recoChargedJetPartEHist->SetLineWidth(2);
0664   recoChargedJetPartENoElecHist->SetLineWidth(2);
0665   recoChargedJetPartEHist->SetTitle("Reconstructed Jet Constituent Energy;Energy [GeV]");
0666 
0667   TLegend *legend7 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0668   legend7->AddEntry(recoChargedJetPartEHist, "With Electrons", "l");
0669   legend7->AddEntry(recoChargedJetPartENoElecHist, "No Electrons", "l");
0670   legend7->Draw();
0671 
0672   gPad->SetLogy();
0673   if(PRINT) c7->Print((results_path+"/recoJetConstituentEnergy.png").c_str()); // Energy of reconstructed jet constituents
0674 
0675   // Reco Part Eta
0676   TCanvas *c8 = new TCanvas("c8","Reco Jet Constituent Eta",800,600);
0677   c8->Clear();
0678   c8->Divide(1,1);
0679 
0680   c8->cd(1);
0681   recoChargedJetPartEtaHist->Draw("HIST");
0682   recoChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
0683   recoChargedJetPartEtaNoElecHist->Draw("HISTSAME");
0684 
0685   recoChargedJetPartEtaHist->SetLineWidth(2);
0686   recoChargedJetPartEtaNoElecHist->SetLineWidth(2);
0687 
0688   recoChargedJetPartEtaHist->SetTitle("Reconstructed Jet Constituent Eta;Eta");
0689 
0690   TLegend *legend8 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0691   legend8->AddEntry(recoChargedJetPartEtaHist, "With Electrons", "l");
0692   legend8->AddEntry(recoChargedJetPartEtaNoElecHist, "No Electrons", "l");
0693   legend8->Draw();
0694 
0695   gPad->SetLogy();
0696   if(PRINT) c8->Print((results_path+"/recoJetConstituentEta.png").c_str()); // Eta of reconstructed jet constituents
0697 
0698   // Reco Part E Vs Eta
0699   TCanvas *c9 = new TCanvas("c9","Reco Jet Constituent E Vs Eta",800,600);
0700   c9->Clear();
0701   c9->Divide(1,1);
0702 
0703   c9->cd(1);
0704   recoChargedJetPartEvsEtaHist->Draw("COLZ");
0705   recoChargedJetPartEvsEtaHist->SetTitle("Reconstructed Jet Constituent Energy Vs Eta;Eta;Energy [GeV]");
0706   gPad->SetLogz();
0707   if(PRINT) c9->Print((results_path+"/recoJetConstituentEnergyVsEta.png").c_str()); // Energy vs eta of reconstructed jet constituents
0708 
0709   // Reco Part Phi Vs Eta
0710   TCanvas *c10 = new TCanvas("c10","Reco Jet Constituent Phi Vs Eta",800,600);
0711   c10->Clear();
0712   c10->Divide(1,1);
0713 
0714   c10->cd(1);
0715   recoChargedJetPartPhiVsEtaHist->Draw("COLZ");
0716   recoChargedJetPartPhiVsEtaHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta;Eta;Phi");
0717   gPad->SetLogz();
0718   if(PRINT) c10->Print((results_path+"/recoJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of reconstructed jet constituents
0719 
0720   // Reco Constituent Pairwise delta R
0721   TCanvas *c11 = new TCanvas("c11","Reco Jet Constituent Pairwise Delta R",800,600);
0722   c11->Clear();
0723   c11->Divide(1,1);
0724 
0725   c11->cd(1);
0726   recoChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
0727   recoChargedJetPartPairwiseDeltaRHist->SetTitle("Pairwise Constituent Delta R;Delta R");
0728   recoChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
0729   gPad->SetLogy();
0730   if(PRINT) c11->Print((results_path+"/recoJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in reconstructed jets
0731 
0732   // Reco E Vs Eta No Electron Jets
0733   TCanvas *c12 = new TCanvas("c12","Reco Jet E Vs Eta (No Electrons)",800,600);
0734   c12->Clear();
0735   c12->Divide(1,1);
0736 
0737   c12->cd(1);
0738   recoChargedJetEvsEtaNoElecHist->Draw("COLZ");
0739   recoChargedJetEvsEtaNoElecHist->SetTitle("Reconstructed Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
0740   gPad->SetLogz();
0741   if(PRINT) c12->Print((results_path+"/recoJetEnergyVsEtaNoElectron.png").c_str()); // Reconstructed jet energy - no jets containing electrons included
0742 
0743   // Reco Phi Vs Eta No Electron Jets
0744   TCanvas *c13 = new TCanvas("c13","Reco Jet Phi Vs Eta (No Electrons)",800,600);
0745   c13->Clear();
0746   c13->Divide(1,1);
0747 
0748   c13->cd(1);
0749   recoChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
0750   recoChargedJetPhiVsEtaECutNoElecHist->SetTitle("Reconstructed Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
0751   gPad->SetLogz();
0752   if(PRINT) c13->Print((results_path+"/recoJetPhiVsEtaNoElectron.png").c_str()); // Reconstructed Jet phi vs eta - no jets containing electrons included
0753 
0754   // Reco Part E Vs Eta No Electron Jets
0755   TCanvas *c14 = new TCanvas("c14","Reco Jet Constituent E Vs Eta (No Electrons)",800,600);
0756   c14->Clear();
0757   c14->Divide(1,1);
0758 
0759   c14->cd(1);
0760   recoChargedJetPartEvsEtaNoElecHist->Draw("COLZ");
0761   recoChargedJetPartEvsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
0762   gPad->SetLogz();
0763   if(PRINT) c14->Print((results_path+"/recoJetConstituentEnergyVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent energy vs eta - no jets containing electrons included
0764 
0765   // Reco Part Phi Vs Eta No Electron Jets
0766   TCanvas *c15 = new TCanvas("c15","Reco Jet Constituent Phi Vs Eta (No Electrons)",800,600);
0767   c15->Clear();
0768   c15->Divide(1,1);
0769 
0770   c15->cd(1);
0771   recoChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
0772   recoChargedJetPartPhiVsEtaNoElecHist->SetTitle("Reconstructed Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
0773   gPad->SetLogz();
0774   if(PRINT) c15->Print((results_path+"/recoJetConstituentPhiVsEtaNoElectron.png").c_str()); // Reconstructed jet constituent phi vs eta - no jets containing electrons included
0775 
0776   
0777   ////////////////////////  Generated Jets Plots  ////////////////////////
0778   // Gen Number
0779   TCanvas *c16 = new TCanvas("c16","Number Gen Jets",800,600);
0780   c16->Clear();
0781   c16->Divide(1,1);
0782 
0783   c16->cd(1);
0784   numGenChargedJetsECutHist->Draw("HIST");
0785   numGenChargedJetsECutNoElecHist->SetLineColor(seabornRed);
0786   numGenChargedJetsECutNoElecHist->Draw("HISTSAME");
0787 
0788   numGenChargedJetsECutHist->SetLineWidth(2);
0789   numGenChargedJetsECutNoElecHist->SetLineWidth(2);
0790 
0791   numGenChargedJetsECutHist->SetTitle("Generator Jets per Event (|eta| < 2.5 && E > 5);Number");
0792 
0793   TLegend *legend16 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0794   legend16->AddEntry(numGenChargedJetsECutHist, "With Electrons", "l");
0795   legend16->AddEntry(numGenChargedJetsECutNoElecHist, "No Electrons", "l");
0796   legend16->Draw();
0797 
0798   gPad->SetLogy();
0799   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
0800 
0801   // Gen Energy
0802   TCanvas *c17 = new TCanvas("c17","Gen Jet Energy",800,600);
0803   c17->Clear();
0804   c17->Divide(1,1);
0805 
0806   c17->cd(1);
0807   genChargedJetEHist->Draw("HIST");
0808   genChargedJetENoElecHist->SetLineColor(seabornRed);
0809   genChargedJetENoElecHist->Draw("HISTSAME");
0810 
0811   genChargedJetEHist->SetLineWidth(2);
0812   genChargedJetENoElecHist->SetLineWidth(2);
0813 
0814   genChargedJetEHist->SetTitle("Generator Jet Energy (|eta| < 2.5);Energy [GeV]");
0815   TLegend *legend17 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0816   legend17->AddEntry(genChargedJetEHist, "With Electrons", "l");
0817   legend17->AddEntry(genChargedJetENoElecHist, "No Electrons", "l");
0818   legend17->Draw();
0819 
0820   gPad->SetLogy();
0821   if(PRINT) c17->Print((results_path+"/genJetEnergy.png").c_str()); // Energy spectrum of generated jets with Abs(eta) < 2.5
0822 
0823   // Gen Eta
0824   TCanvas *c18 = new TCanvas("c18","Gen Jet Eta",800,600);
0825   c18->Clear();
0826   c18->Divide(1,1);
0827 
0828   c18->cd(1);
0829   genChargedJetEtaECutHist->Draw("HIST");
0830   genChargedJetEtaECutNoElecHist->SetLineColor(seabornRed);
0831   genChargedJetEtaECutNoElecHist->Draw("HISTSAME");
0832 
0833   genChargedJetEtaECutHist->SetLineWidth(2);
0834   genChargedJetEtaECutNoElecHist->SetLineWidth(2);
0835 
0836   genChargedJetEtaECutHist->SetTitle("Generator Jet Eta (E > 5);Eta");
0837 
0838   TLegend *legend18 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0839   legend18->AddEntry(genChargedJetEtaECutHist, "With Electrons", "l");
0840   legend18->AddEntry(genChargedJetEtaECutNoElecHist, "No Electrons", "l");
0841   legend18->Draw();
0842 
0843   gPad->SetLogy();
0844   if(PRINT) c18->Print((results_path+"/genJetEta.png").c_str()); // Eta spectrum of generator jets with energy > 5 GeV
0845 
0846   // Gen E Vs Eta
0847   TCanvas *c19 = new TCanvas("c19","Gen Jet E Vs Eta",800,600);
0848   c19->Clear();
0849   c19->Divide(1,1);
0850 
0851   c19->cd(1);
0852   genChargedJetEvsEtaHist->Draw("COLZ");
0853   genChargedJetEvsEtaHist->SetTitle("Generator Jet Energy Vs Eta;Eta;Energy [GeV]");
0854   gPad->SetLogz();
0855   if(PRINT) c19->Print((results_path+"/genJetEnergyvsEta.png").c_str()); // Energy vs eta of generator jets
0856 
0857   // Gen Phi Vs Eta
0858   TCanvas *c20 = new TCanvas("c20","Gen Jet Phi Vs Eta",800,600);
0859   c20->Clear();
0860   c20->Divide(1,1);
0861 
0862   c20->cd(1);
0863   genChargedJetPhiVsEtaECutHist->Draw("COLZ");
0864   genChargedJetPhiVsEtaECutHist->SetTitle("Generator Jet Phi Vs Eta (E > 5);Eta;Phi");
0865   gPad->SetLogz();
0866   if(PRINT) c20->Print((results_path+"/genJetPhiVsEta.png").c_str()); // Phi vs eta of generator jets
0867 
0868   // Num Particles Per Gen Jet
0869   TCanvas *c21 = new TCanvas("c21","Number Constituents Per Gen Jet",800,600);
0870   c21->Clear();
0871   c21->Divide(1,1);
0872 
0873   c21->cd(1);
0874   numGenChargedJetPartsHist->Draw("HIST");
0875   numGenChargedJetPartsNoElecHist->SetLineColor(seabornRed);
0876   numGenChargedJetPartsNoElecHist->Draw("HISTSAME");
0877 
0878   numGenChargedJetPartsHist->SetLineWidth(2);
0879   numGenChargedJetPartsNoElecHist->SetLineWidth(2);
0880 
0881   numGenChargedJetPartsHist->SetTitle("Number of Constituents Per Gen Jet;Number of Constituents");
0882 
0883   TLegend *legend21 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0884   legend21->AddEntry(numGenChargedJetPartsHist, "With Electrons", "l");
0885   legend21->AddEntry(numGenChargedJetPartsNoElecHist, "No Electrons", "l");
0886   legend21->Draw();
0887   gPad->SetLogy();
0888   if(PRINT) c21->Print((results_path+"/numConstituentsPerGenJet.png").c_str()); // Number of constituents in generator jets
0889 
0890   // Gen Part Energy
0891   TCanvas *c22 = new TCanvas("c22","Gen Jet Constituent Energy",800,600);
0892   c22->Clear();
0893   c22->Divide(1,1);
0894 
0895   c22->cd(1);
0896   genChargedJetPartEHist->Draw("HIST");
0897   genChargedJetPartENoElecHist->SetLineColor(seabornRed);
0898   genChargedJetPartENoElecHist->Draw("HISTSAME");
0899 
0900   genChargedJetPartEHist->SetLineWidth(2);
0901   genChargedJetPartENoElecHist->SetLineWidth(2);
0902 
0903   genChargedJetPartEHist->SetTitle("Generator Jet Constituent Energy;Energy [GeV]");
0904 
0905   TLegend *legend22 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0906   legend22->AddEntry(genChargedJetPartEHist, "With Electrons", "l");
0907   legend22->AddEntry(genChargedJetPartENoElecHist, "No Electrons", "l");
0908   legend22->Draw();
0909 
0910   gPad->SetLogy();
0911   if(PRINT) c22->Print((results_path+"/genJetConstituentEnergy.png").c_str()); // Energy of generator jet constituents
0912 
0913   // Gen Part Eta
0914   TCanvas *c23 = new TCanvas("c23","Gen Jet Constituent Eta",800,600);
0915   c23->Clear();
0916   c23->Divide(1,1);
0917 
0918   c23->cd(1);
0919   genChargedJetPartEtaHist->Draw("HIST");
0920   genChargedJetPartEtaNoElecHist->SetLineColor(seabornRed);
0921   genChargedJetPartEtaNoElecHist->Draw("HISTSAME");
0922 
0923   genChargedJetPartEtaHist->SetLineWidth(2);
0924   genChargedJetPartEtaNoElecHist->SetLineWidth(2);
0925 
0926   genChargedJetPartEtaHist->SetTitle("Generator Jet Constituent Eta;Eta");
0927 
0928   TLegend *legend23 = new TLegend(0.7, 0.7, 0.9, 0.9); // Adjust the coordinates as needed
0929   legend23->AddEntry(genChargedJetPartEtaHist, "With Electrons", "l");
0930   legend23->AddEntry(genChargedJetPartEtaNoElecHist, "No Electrons", "l");
0931   legend23->Draw();
0932 
0933   gPad->SetLogy();
0934   if(PRINT) c23->Print((results_path+"/genJetConstituentEta.png").c_str()); // Eta of generator jet constituents
0935 
0936   // Gen Part E Vs Eta
0937   TCanvas *c24 = new TCanvas("c24","Gen Jet Constituent E Vs Eta",800,600);
0938   c24->Clear();
0939   c24->Divide(1,1);
0940 
0941   c24->cd(1);
0942   genChargedJetPartEvsEtaHist->Draw("COLZ");
0943   genChargedJetPartEvsEtaHist->SetTitle("Generator Jet Constituent Energy Vs Eta;Eta;Energy [GeV]");
0944   gPad->SetLogz();
0945   if(PRINT) c24->Print((results_path+"/genJetConstituentEnergyVsEta.png").c_str()); // Energy vs eta of generator jet constituents
0946 
0947   // Gen Part Phi Vs Eta
0948   TCanvas *c25 = new TCanvas("c25","Gen Jet Constituent Phi Vs Eta",800,600);
0949   c25->Clear();
0950   c25->Divide(1,1);
0951 
0952   c25->cd(1);
0953   genChargedJetPartPhiVsEtaHist->Draw("COLZ");
0954   genChargedJetPartPhiVsEtaHist->SetTitle("Generator Jet Constituent Phi Vs Eta;Eta;Phi");
0955   gPad->SetLogz();
0956   if(PRINT) c25->Print((results_path+"/genJetConstituentPhiVsEta.png").c_str()); // Phi vs eta of generator jet constituents
0957 
0958   // Gen Constituent Pairwise delta R
0959   TCanvas *c26 = new TCanvas("c26","Gen Jet Constituent Pairwise Delta R",800,600);
0960   c26->Clear();
0961   c26->Divide(1,1);
0962 
0963   c26->cd(1);
0964   genChargedJetPartPairwiseDeltaRHist->Draw("COLZ");
0965   genChargedJetPartPairwiseDeltaRHist->SetTitle("Generator Jet Pairwise Constituent Delta R;Delta R");
0966   genChargedJetPartPairwiseDeltaRHist->GetXaxis()->SetRangeUser(0,0.5);
0967   gPad->SetLogy();
0968   if(PRINT) c26->Print((results_path+"/genJetConstituentPairwiseDR.png").c_str()); // Distance between each pair of constituents in generator jets
0969 
0970   // Gen E Vs Eta No Electron Jets
0971   TCanvas *c27 = new TCanvas("c27","Gen Jet E Vs Eta (No Electrons)",800,600);
0972   c27->Clear();
0973   c27->Divide(1,1);
0974 
0975   c27->cd(1);
0976   genChargedJetEvsEtaNoElecHist->Draw("COLZ");
0977   genChargedJetEvsEtaNoElecHist->SetTitle("Generator Jet Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
0978   gPad->SetLogz();
0979   if(PRINT) c27->Print((results_path+"/genJetEnergyVsEtaNoElectron.png").c_str()); // Generator jet energy vs eta - no jets containing electrons included
0980 
0981   // Gen Phi Vs Eta No Electron Jets
0982   TCanvas *c28 = new TCanvas("c28","Gen Jet Phi Vs Eta (No Electrons)",800,600);
0983   c28->Clear();
0984   c28->Divide(1,1);
0985 
0986   c28->cd(1);
0987   genChargedJetPhiVsEtaECutNoElecHist->Draw("COLZ");
0988   genChargedJetPhiVsEtaECutNoElecHist->SetTitle("Generator Jet Phi Vs Eta (E > 5) (No Electrons);Eta;Phi");
0989   gPad->SetLogz();
0990   if(PRINT) c28->Print((results_path+"/genJetPhiVsEtaNoElectron.png").c_str()); // Generator Jet phi vs eta - no jets containing electrons included
0991 
0992   // Gen Part E Vs Eta No Electron Jets
0993   TCanvas *c29 = new TCanvas("c29","Gen Jet Constituent E Vs Eta (No Electrons)",800,600);
0994   c29->Clear();
0995   c29->Divide(1,1);
0996 
0997   c29->cd(1);
0998   genChargedJetPartEvsEtaNoElecHist->Draw("COLZ");
0999   genChargedJetPartEvsEtaNoElecHist->SetTitle("Generator Jet Constituent Energy Vs Eta (No Electrons);Eta;Energy [GeV]");
1000   gPad->SetLogz();
1001   if(PRINT) c29->Print((results_path+"/genJetConstituentEnergyVsEtaNoElectron.png").c_str()); // Generator jet constituent energy vs eta - no jets containing electrons included
1002 
1003   // Gen Part Phi Vs Eta No Electron Jets
1004   TCanvas *c30 = new TCanvas("c30","Gen Jet Constituent Phi Vs Eta (No Electrons)",800,600);
1005   c30->Clear();
1006   c30->Divide(1,1);
1007 
1008   c30->cd(1);
1009   genChargedJetPartPhiVsEtaNoElecHist->Draw("COLZ");
1010   genChargedJetPartPhiVsEtaNoElecHist->SetTitle("Generator Jet Constituent Phi Vs Eta (No Electrons);Eta;Phi");
1011   gPad->SetLogz();
1012   //c30->Print((results_path+"/recoJetEvsEta.png").c_str());
1013   if(PRINT) c30->Print((results_path+"/genJetConstituentPhiVsEtaNoElectron.png").c_str()); // Generator jet constituent phi vs eta - no jets containing electrons included
1014 
1015   
1016   ////////////////////////  Matched Jets Plots  ////////////////////////
1017   // Matched Delta R
1018   TCanvas *c31 = new TCanvas("c31","Gen - Reco Delta R",800,600);
1019   c31->Clear();
1020   c31->Divide(1,1);
1021 
1022   c31->cd(1);
1023   matchJetDeltaRHist->Draw("HIST");
1024   matchJetDeltaRBackHist->SetLineColor(seabornRed);
1025   //matchJetDeltaRBackHist->Draw("HISTSAME");
1026   matchJetDeltaRHist->SetTitle("Matched Gen - Reco Jet Delta R;Delta R");
1027   gPad->SetLogy();
1028   if(PRINT) c31->Print((results_path+"/genRecoJetDeltaR.png").c_str()); // Distance between closest generated and reconstructed jet pair
1029 
1030   // Matched Reco Vs Gen Eta
1031   TCanvas *c32 = new TCanvas("c32","Reco Vs Gen Eta",800,600);
1032   c32->Clear();
1033   c32->Divide(1,1);
1034 
1035   c32->cd(1);
1036   recoVsGenChargedJetEtaHist->Draw("COLZ");
1037   recoVsGenChargedJetEtaHist->SetTitle("Reconstructed Vs Generator Jet Eta;Gen Eta;Reco Eta");
1038   gPad->SetLogz();
1039   if(PRINT) c32->Print((results_path+"/matchedRecoVsGenJetEta.png").c_str()); // Matched Reconstructed Vs Generator Jet eta
1040 
1041   // Matched Reco Vs Gen Phi
1042   TCanvas *c33 = new TCanvas("c33","Reco Vs Gen Phi",800,600);
1043   c33->Clear();
1044   c33->Divide(1,1);
1045 
1046   c33->cd(1);
1047   recoVsGenChargedJetPhiHist->Draw("COLZ");
1048   recoVsGenChargedJetPhiHist->SetTitle("Reconstructed Vs Generator Jet Phi;Gen Phi;Reco Phi");
1049   gPad->SetLogz();
1050   if(PRINT) c33->Print((results_path+"/matchedRecoVsGenJetPhi.png").c_str()); // Matched reconstructed vs generator jet phi
1051 
1052   // Matched Reco Vs Gen Energy
1053   TCanvas *c34 = new TCanvas("c34","Reco Vs Gen Energy",800,600);
1054   c34->Clear();
1055   c34->Divide(1,1);
1056 
1057   TF1 *f1_34 = new TF1("f1_34","1.0*x + 0.0",1,100);
1058   TF1 *f2_34 = new TF1("f2_34","2.0*x + 0.0",1,100);
1059   TF1 *f3_34 = new TF1("f3_34","3.0*x + 0.0",1,100);
1060 
1061   c34->cd(1);
1062   recoVsGenChargedJetEHist->Draw("COLZ");
1063   recoVsGenChargedJetEHist->SetTitle("Reconstructed Vs Generator Jet Energy;Gen E;Reco E");
1064   f1_34->Draw("SAME");
1065   f2_34->Draw("SAME");
1066   f3_34->Draw("SAME");
1067   gPad->SetLogz();
1068   if(PRINT) c34->Print((results_path+"/matchedRecoVsGenJetEnergy.png").c_str()); // Matched reconstructed vs generator jet energy
1069 
1070   // Jet Res Vs Gen Eta
1071   TCanvas *c35 = new TCanvas("c35","Jet Res Vs Gen Eta",800,600);
1072   c35->Clear();
1073   c35->Divide(1,1);
1074 
1075   c35->cd(1);
1076   jetResVsEtaHist->Draw("COLZ");
1077   jetResVsEtaHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Eta;Gen Eta;Res");
1078   gPad->SetLogz();
1079   if(PRINT) c35->Print((results_path+"/matchedJetResolutionVsEta.png").c_str()); // Matched jet resolution vs generator jet eta
1080 
1081   // Jet Res Vs Gen E
1082   TCanvas *c36 = new TCanvas("c36","Jet Res Vs Gen E",800,600);
1083   c36->Clear();
1084   c36->Divide(1,1);
1085 
1086   c36->cd(1);
1087   jetResVsEHist->Draw("COLZ");
1088   jetResVsEHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy;Gen E;Res");
1089   gPad->SetLogz();
1090   if(PRINT) c36->Print((results_path+"/matchedJetResolutionVsEnergy.png").c_str()); // Matched jet resolution vs generator jet energy
1091 
1092   // Jet Res Vs Gen E Neg Eta
1093   TCanvas *c37 = new TCanvas("c37","Jet Res Vs Gen E (-2.5 < eta < -1.0)",800,600);
1094   c37->Clear();
1095   c37->Divide(1,1);
1096 
1097   c37->cd(1);
1098   jetResVsENegEtaNoDupHist->Draw("COLZ");
1099   jetResVsENegEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-2.5 < eta < -1.0) No Duplicate;Gen E;Res");
1100   gPad->SetLogz();
1101   if(PRINT) c37->Print((results_path+"/matchedJetResolutionVsEnergyNegEta.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0
1102 
1103   // Jet Res Vs Gen E Mid Eta
1104   TCanvas *c38 = new TCanvas("c38","Jet Res Vs Gen E (-1.0 < eta < 1.0)",800,600);
1105   c38->Clear();
1106   c38->Divide(1,1);
1107 
1108   c38->cd(1);
1109   jetResVsEMidEtaNoDupHist->Draw("COLZ");
1110   jetResVsEMidEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (-1.0 < eta < 1.0) No Duplicate;Gen E;Res");
1111   gPad->SetLogz();
1112   if(PRINT) c38->Print((results_path+"/matchedJetResolutionVsEnergyMidEta.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0
1113     delete c38;
1114   // Jet Res Vs Gen E Pos Eta
1115   TCanvas *c39 = new TCanvas("c39","Jet Res Vs Gen E (1.0 < eta < 2.5)",800,600);
1116   c39->Clear();
1117   c39->Divide(1,1);
1118 
1119   c39->cd(1);
1120   jetResVsEPosEtaNoDupHist->Draw("COLZ");
1121   jetResVsEPosEtaNoDupHist->SetTitle("(Reco - Gen)/Gen Jet Energy Vs Gen Energy (1.0 < eta < 2.5) No Duplicate;Gen E;Res");
1122   gPad->SetLogz();
1123   if(PRINT) c39->Print((results_path+"/matchedJetResolutionVsEnergyPosEta.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5
1124   delete c39;
1125 
1126   
1127   // Generate Resolution Plots
1128   const int BINS = 20;
1129   double binCent[BINS];
1130   double jesVsENeg[BINS];
1131   double jesVsEMid[BINS];
1132   double jesVsEPos[BINS];
1133   double jerVsENeg[BINS];
1134   double jerVsEMid[BINS];
1135   double jerVsEPos[BINS];
1136 
1137   std::fill(std::begin(binCent), std::end(binCent), -999.);
1138   std::fill(std::begin(jesVsENeg), std::end(jesVsENeg), -999.);
1139   std::fill(std::begin(jesVsEMid), std::end(jesVsEMid), -999.);
1140   std::fill(std::begin(jesVsEPos), std::end(jesVsEPos), -999.);
1141   std::fill(std::begin(jerVsENeg), std::end(jerVsENeg), -999.);
1142   std::fill(std::begin(jerVsEMid), std::end(jerVsEMid), -999.);
1143   std::fill(std::begin(jerVsEPos), std::end(jerVsEPos), -999.);
1144 
1145   TH1D *pxA = jetResVsENegEtaNoDupHist->ProjectionX("pxA",1,10000);
1146   for(int i=0; i<BINS; i++)
1147     {
1148       binCent[i] = pxA->GetBinCenter(i+1);
1149     }
1150 
1151   TCanvas *c40 = new TCanvas("c40","Negative Rapidity Fit Results",800,600);
1152   c40->Clear();
1153   c40->Divide(5,4);
1154 
1155   TH1D *hA[20];
1156   for(int i=1; i<21; i++)
1157     {
1158       hA[i-1] = (TH1D *)jetResVsENegEtaNoDupHist->ProjectionY(Form("projYA_%d",i),i,i);
1159 
1160       TF1 *myGausA = new TF1("myGausA","gaus",-0.5,0.5);
1161       myGausA->SetParameters(hA[i-1]->GetMaximum(),0.0,0.01);
1162 
1163       c40->cd(i);
1164       //hA[i-1]->Draw("HIST");
1165       hA[i-1]->Fit("myGausA","B","",-0.5,0.5);
1166       hA[i-1]->GetXaxis()->SetRangeUser(-1,1);
1167       gPad->SetLogy();
1168 
1169       if(hA[i-1]->GetEntries() > 2)
1170     {
1171       auto fA = hA[i-1]->GetFunction("myGausA");
1172 
1173       if(fA->GetParError(2)/fA->GetParameter(2) < 0.5)
1174         {
1175           jesVsENeg[i-1] = fA->GetParameter(1);
1176           jerVsENeg[i-1] = fA->GetParameter(2);
1177         }
1178       //cout << fA->GetParameter(0) << " " << fA->GetParameter(1) << " " << fA->GetParameter(2) << endl;
1179       //cout << fA->GetParError(0) << " " << fA->GetParError(1) << " " << fA->GetParError(2) << endl;
1180     }
1181     }
1182   if(PRINT) c40->Print((results_path+"/matchedJetResolutionVsEnergyNegEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -2.5 < eta < -1.0 fits
1183     delete c40;
1184 
1185   TCanvas *c41 = new TCanvas("c41","Mid Rapidity Fit Results",800,600);
1186   c41->Clear();
1187   c41->Divide(5,4);
1188 
1189   TH1D *hB[20];
1190   for(int i=1; i<21; i++)
1191     {
1192       hB[i-1] = (TH1D *)jetResVsEMidEtaNoDupHist->ProjectionY(Form("projYB_%d",i),i,i);
1193 
1194       TF1 *myGausB = new TF1("myGausB","gaus",-0.5,0.5);
1195       myGausB->SetParameters(hB[i-1]->GetMaximum(),0.0,0.01);
1196 
1197       c41->cd(i);
1198       //hA[i-1]->Draw("HIST");
1199       hB[i-1]->Fit("myGausB","B","",-0.5,0.5);
1200       hB[i-1]->GetXaxis()->SetRangeUser(-1,1);
1201       gPad->SetLogy();
1202 
1203       if(hB[i-1]->GetEntries() > 2)
1204     {
1205       auto fB = hB[i-1]->GetFunction("myGausB");
1206 
1207       if(fB->GetParError(2)/fB->GetParameter(2) < 0.5)
1208         {
1209           jesVsEMid[i-1] = fB->GetParameter(1);
1210           jerVsEMid[i-1] = fB->GetParameter(2);
1211         }
1212       //cout << fB->GetParameter(0) << " " << fB->GetParameter(1) << " " << fB->GetParameter(2) << endl;
1213       //cout << fB->GetParError(0) << " " << fB->GetParError(1) << " " << fB->GetParError(2) << endl;
1214     }
1215     }
1216   if(PRINT) c41->Print((results_path+"/matchedJetResolutionVsEnergyMidEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy -1.0 < eta < 1.0 fits
1217     delete c41;
1218 
1219   TCanvas *c42 = new TCanvas("c42","Positive Rapidity Fit Results",800,600);
1220   c42->Clear();
1221   c42->Divide(5,4);
1222 
1223   TH1D *hC[20];
1224   for(int i=1; i<21; i++)
1225     {
1226       hC[i-1] = (TH1D *)jetResVsEPosEtaNoDupHist->ProjectionY(Form("projYC_%d",i),i,i);
1227 
1228       TF1 *myGausC = new TF1("myGausC","gaus",-0.5,0.5);
1229       myGausC->SetParameters(hC[i-1]->GetMaximum(),0.0,0.01);
1230 
1231       c42->cd(i);
1232       //hA[i-1]->Draw("HIST");
1233       hC[i-1]->Fit("myGausC","B","",-0.5,0.5);
1234       hC[i-1]->GetXaxis()->SetRangeUser(-1,1);
1235       gPad->SetLogy();
1236 
1237       if(hC[i-1]->GetEntries() > 2)
1238     {
1239       auto fC = hC[i-1]->GetFunction("myGausC");
1240 
1241       if(fC->GetParError(2)/fC->GetParameter(2) < 0.5)
1242         {
1243           jesVsEPos[i-1] = fC->GetParameter(1);
1244           jerVsEPos[i-1] = fC->GetParameter(2);
1245         }
1246       //cout << fC->GetParameter(0) << " " << fC->GetParameter(1) << " " << fC->GetParameter(2) << endl;
1247       //cout << fC->GetParError(0) << " " << fC->GetParError(1) << " " << fC->GetParError(2) << endl;
1248     }
1249     }
1250   if(PRINT) c42->Print((results_path+"/matchedJetResolutionVsEnergyPosEtaFitSummary.png").c_str()); // Matched jet resolution vs generator jet energy 1.0 < eta < 2.5 fits
1251     delete c42;
1252   TCanvas *c43 = new TCanvas("c43","Positive JES/JER",800,600);
1253   c43->Clear();
1254   c43->Divide(1,1);
1255 
1256   TGraph *gJESvsENeg = new TGraph(BINS,binCent,jesVsENeg);
1257   TGraph *gJERvsENeg = new TGraph(BINS,binCent,jerVsENeg);
1258 
1259   TGraph *gJESvsEMid = new TGraph(BINS,binCent,jesVsEMid);
1260   TGraph *gJERvsEMid = new TGraph(BINS,binCent,jerVsEMid);
1261 
1262   TGraph *gJESvsEPos = new TGraph(BINS,binCent,jesVsEPos);
1263   TGraph *gJERvsEPos = new TGraph(BINS,binCent,jerVsEPos);
1264 
1265   TH2D *test43 = new TH2D("test43","Jet Energy Scale / Resolution Vs Eta;True Eta;JES/JER",1,0.,100.,1,-0.2,0.2);
1266   test43->Draw();
1267 
1268   c43->cd(1);
1269   gJERvsENeg->Draw("*");
1270   gJERvsENeg->SetMarkerStyle(21);
1271   gJERvsENeg->SetMarkerSize(1);
1272   gJERvsENeg->SetMarkerColor(seabornBlue);
1273 
1274   gJESvsENeg->Draw("*");
1275   gJESvsENeg->SetMarkerStyle(26);
1276   gJESvsENeg->SetMarkerSize(1);
1277   gJESvsENeg->SetMarkerColor(seabornBlue);
1278 
1279   gJERvsEMid->Draw("*");
1280   gJERvsEMid->SetMarkerStyle(21);
1281   gJERvsEMid->SetMarkerSize(1);
1282   gJERvsEMid->SetMarkerColor(seabornRed);
1283 
1284   gJESvsEMid->Draw("*");
1285   gJESvsEMid->SetMarkerStyle(26);
1286   gJESvsEMid->SetMarkerSize(1);
1287   gJESvsEMid->SetMarkerColor(seabornRed);
1288 
1289   gJERvsEPos->Draw("*");
1290   gJERvsEPos->SetMarkerStyle(21);
1291   gJERvsEPos->SetMarkerSize(1);
1292   gJERvsEPos->SetMarkerColor(seabornGreen);
1293 
1294   gJESvsEPos->Draw("*");
1295   gJESvsEPos->SetMarkerStyle(26);
1296   gJESvsEPos->SetMarkerSize(1);
1297   gJESvsEPos->SetMarkerColor(seabornGreen);
1298 
1299   TLegend *legend = new TLegend(0.7,0.7,0.9,0.9); 
1300 legend->AddEntry(gJERvsENeg, "JER, (-2.5 < #eta < -1)", "p");
1301 legend->AddEntry(gJESvsENeg, "JES, (-2.5 < #eta < -1)","p");
1302 legend->AddEntry(gJERvsEMid, "JER, (-1 < #eta < 1)", "p");
1303 legend->AddEntry(gJESvsEMid, "JES, (-1 < #eta < 1)", "p");
1304 legend->AddEntry(gJERvsEPos, "JER, (1 < #eta < 2.5) ", "p");
1305 legend->AddEntry(gJESvsEPos, "JES, (1 < #eta < 2.5)", "p");
1306 legend->Draw();
1307 
1308   if(PRINT) c43->Print((results_path+"/matchedJetScaleResolutionSummary.png").c_str()); // Matched jet JER/JES summary
1309     delete c43;
1310 
1311 delete mychain;
1312 
1313   return 0;
1314 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!