Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-12-18 10:45:04 

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