File indexing completed on 2024-04-06 12:25:27
0001
0002
0003
0004
0005
0006 #include "RecoJets/JetAnalyzers/interface/JetValidation.h"
0007 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
0008 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0009 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0010 #include "DataFormats/JetReco/interface/CaloJet.h"
0011 #include "DataFormats/JetReco/interface/GenJet.h"
0012 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include <TROOT.h>
0017 #include <TSystem.h>
0018 #include <TFile.h>
0019 #include <cmath>
0020 using namespace edm;
0021 using namespace reco;
0022 using namespace std;
0023
0024 JetValidation::JetValidation(edm::ParameterSet const& cfg) {
0025 dRmatch = cfg.getParameter<double>("dRmatch");
0026 PtMin = cfg.getParameter<double>("PtMin");
0027 Njets = cfg.getParameter<int>("Njets");
0028 MCarlo = cfg.getParameter<bool>("MCarlo");
0029 genAlgo = cfg.getParameter<string>("genAlgo");
0030 calAlgo = cfg.getParameter<string>("calAlgo");
0031 jetTracksAssociator = cfg.getParameter<string>("jetTracksAssociator");
0032 histoFileName = cfg.getParameter<string>("histoFileName");
0033 }
0034
0035 void JetValidation::beginJob() {
0036 m_file = new TFile(histoFileName.c_str(), "RECREATE");
0037
0038 m_HistNames1D["CaloJetMulti"] = new TH1F("CaloJetMulti", "Multiplicity of CaloJets", 100, 0, 100);
0039 m_HistNames1D["ptCalo"] = new TH1F("ptCalo", "p_{T} of CaloJets", 7000, 0, 7000);
0040 m_HistNames1D["etaCalo"] = new TH1F("etaCalo", "#eta of CaloJets", 100, -5.0, 5.0);
0041 m_HistNames1D["phiCalo"] = new TH1F("phiCalo", "#phi of CaloJets", 72, -M_PI, M_PI);
0042 m_HistNames1D["m2jCalo"] = new TH1F("m2jCalo", "Dijet Mass of leading CaloJets", 7000, 0, 14000);
0043 m_HistNames1D["nTracks"] = new TH1F("nTracks", "Number of tracks associated with a jet", 100, 0, 100);
0044 m_HistNames1D["chargeFraction"] = new TH1F("chargeFraction", "Fraction of charged tracks pt", 500, 0, 5);
0045 m_HistNames1D["emEnergyFraction"] = new TH1F("emEnergyFraction", "Jets EM Fraction", 110, 0, 1.1);
0046 m_HistNames1D["emEnergyInEB"] = new TH1F("emEnergyInEB", "Jets emEnergyInEB", 7000, 0, 14000);
0047 m_HistNames1D["emEnergyInEE"] = new TH1F("emEnergyInEE", "Jets emEnergyInEE", 7000, 0, 14000);
0048 m_HistNames1D["emEnergyInHF"] = new TH1F("emEnergyInHF", "Jets emEnergyInHF", 7000, 0, 14000);
0049 m_HistNames1D["hadEnergyInHB"] = new TH1F("hadEnergyInHB", "Jets hadEnergyInHB", 7000, 0, 14000);
0050 m_HistNames1D["hadEnergyInHE"] = new TH1F("hadEnergyInHE", "Jets hadEnergyInHE", 7000, 0, 14000);
0051 m_HistNames1D["hadEnergyInHF"] = new TH1F("hadEnergyInHF", "Jets hadEnergyInHF", 7000, 0, 14000);
0052 m_HistNames1D["hadEnergyInHO"] = new TH1F("hadEnergyInHO", "Jets hadEnergyInHO", 7000, 0, 14000);
0053 m_HistNamesProfile["EBfractionVsEta"] = new TProfile("EBfractionVsEta", "Jets EBfraction vs #eta", 100, -5.0, 5.0);
0054 m_HistNamesProfile["EEfractionVsEta"] = new TProfile("EEfractionVsEta", "Jets EEfraction vs #eta", 100, -5.0, 5.0);
0055 m_HistNamesProfile["HBfractionVsEta"] = new TProfile("HBfractionVsEta", "Jets HBfraction vs #eta", 100, -5.0, 5.0);
0056 m_HistNamesProfile["HOfractionVsEta"] = new TProfile("HOfractionVsEta", "Jets HOfraction vs #eta", 100, -5.0, 5.0);
0057 m_HistNamesProfile["HEfractionVsEta"] = new TProfile("HEfractionVsEta", "Jets HEfraction vs #eta", 100, -5.0, 5.0);
0058 m_HistNamesProfile["HFfractionVsEta"] = new TProfile("HFfractionVsEta", "Jets HFfraction vs #eta", 100, -5.0, 5.0);
0059 m_HistNamesProfile["CaloEnergyVsEta"] = new TProfile("CaloEnergyVsEta", "CaloJets Energy Vs. Eta", 100, -5.0, 5.0);
0060 m_HistNamesProfile["emEnergyVsEta"] = new TProfile("emEnergyVsEta", "Jets EM Energy Vs. Eta", 100, -5.0, 5.0);
0061 m_HistNamesProfile["hadEnergyVsEta"] = new TProfile("hadEnergyVsEta", "Jets HAD Energy Vs. Eta", 100, -5.0, 5.0);
0062 if (MCarlo) {
0063 m_HistNames1D["GenJetMulti"] = new TH1F("GenJetMulti", "Multiplicity of GenJets", 100, 0, 100);
0064 m_HistNames1D["ptHat"] = new TH1F("ptHat", "p_{T}hat", 7000, 0, 7000);
0065 m_HistNames1D["ptGen"] = new TH1F("ptGen", "p_{T} of GenJets", 7000, 0, 7000);
0066 m_HistNames1D["etaGen"] = new TH1F("etaGen", "#eta of GenJets", 100, -5.0, 5.0);
0067 m_HistNames1D["phiGen"] = new TH1F("phiGen", "#phi of GenJets", 72, -M_PI, M_PI);
0068 m_HistNames1D["m2jGen"] = new TH1F("m2jGen", "Dijet Mass of leading GenJets", 7000, 0, 14000);
0069 m_HistNames1D["dR"] = new TH1F("dR", "GenJets dR with matched CaloJet", 200, 0, 1);
0070 m_HistNamesProfile["GenEnergyVsEta"] = new TProfile("GenEnergyVsEta", "GenJets Energy Vs. Eta", 100, -5.0, 5.0);
0071 m_HistNamesProfile["respVsPtBarrel"] =
0072 new TProfile("respVsPtBarrel", "CaloJet Response of GenJets in Barrel", 7000, 0, 7000);
0073 m_HistNamesProfile["CaloErespVsEta"] =
0074 new TProfile("CaloErespVsEta", "Jets Energy Response Vs. Eta", 100, -5.0, 5.0);
0075 m_HistNamesProfile["emErespVsEta"] =
0076 new TProfile("emErespVsEta", "Jets EM Energy Response Vs. Eta", 100, -5.0, 5.0);
0077 m_HistNamesProfile["hadErespVsEta"] =
0078 new TProfile("hadErespVsEta", "Jets HAD Energy Response Vs. Eta", 100, -5.0, 5.0);
0079 }
0080 }
0081
0082 void JetValidation::analyze(edm::Event const& evt, edm::EventSetup const& iSetup) {
0083 math::XYZTLorentzVector p4jet[2];
0084 int jetInd, jetCounter, nTracks;
0085 double dRmin, dR, e, eta, emEB, emEE, emHF, hadHB, hadHE, hadHO, hadHF, pt, phi, pthat, chf;
0086 Handle<CaloJetCollection> caljets;
0087 Handle<GenJetCollection> genjets;
0088 Handle<double> genEventScale;
0089 Handle<JetTracksAssociation::Container> jetTracks;
0090 CaloJetCollection::const_iterator i_caljet;
0091 GenJetCollection::const_iterator i_genjet;
0092 evt.getByLabel(calAlgo, caljets);
0093 evt.getByLabel(jetTracksAssociator, jetTracks);
0094 jetInd = 0;
0095 jetCounter = 0;
0096 if (caljets->empty())
0097 cout << "WARNING: NO calo jets in event " << evt.id().event() << ", Run " << evt.id().run() << " !!!!" << endl;
0098 for (i_caljet = caljets->begin(); i_caljet != caljets->end() && jetInd < Njets; ++i_caljet) {
0099 e = i_caljet->energy();
0100 pt = i_caljet->pt();
0101 phi = i_caljet->phi();
0102 eta = i_caljet->eta();
0103 emEB = i_caljet->emEnergyInEB();
0104 emEE = i_caljet->emEnergyInEE();
0105 emHF = i_caljet->emEnergyInHF();
0106 hadHB = i_caljet->hadEnergyInHB();
0107 hadHE = i_caljet->hadEnergyInHE();
0108 hadHO = i_caljet->hadEnergyInHO();
0109 hadHF = i_caljet->hadEnergyInHF();
0110 nTracks = JetTracksAssociation::tracksNumber(*jetTracks, *i_caljet);
0111 chf = (JetTracksAssociation::tracksP4(*jetTracks, *i_caljet)).pt() / pt;
0112 if (jetInd < 2)
0113 p4jet[jetInd] = i_caljet->p4();
0114 if (pt > PtMin) {
0115 FillHist1D("ptCalo", pt);
0116 FillHist1D("etaCalo", eta);
0117 FillHist1D("phiCalo", phi);
0118 FillHist1D("emEnergyFraction", i_caljet->emEnergyFraction());
0119 FillHist1D("nTracks", nTracks);
0120 FillHist1D("chargeFraction", chf);
0121 FillHist1D("emEnergyInEB", emEB);
0122 FillHist1D("emEnergyInEE", emEE);
0123 FillHist1D("emEnergyInHF", emHF);
0124 FillHist1D("hadEnergyInHB", hadHB);
0125 FillHist1D("hadEnergyInHE", hadHE);
0126 FillHist1D("hadEnergyInHF", hadHF);
0127 FillHist1D("hadEnergyInHO", hadHO);
0128 FillHistProfile("EBfractionVsEta", eta, emEB / e);
0129 FillHistProfile("EEfractionVsEta", eta, emEE / e);
0130 FillHistProfile("HBfractionVsEta", eta, hadHB / e);
0131 FillHistProfile("HOfractionVsEta", eta, hadHO / e);
0132 FillHistProfile("HEfractionVsEta", eta, hadHE / e);
0133 FillHistProfile("HFfractionVsEta", eta, (hadHF + emHF) / e);
0134 FillHistProfile("CaloEnergyVsEta", eta, e);
0135 FillHistProfile("emEnergyVsEta", eta, emEB + emEE + emHF);
0136 FillHistProfile("hadEnergyVsEta", eta, hadHB + hadHO + hadHE + hadHF);
0137 jetCounter++;
0138 }
0139 jetInd++;
0140 }
0141 FillHist1D("CaloJetMulti", jetCounter);
0142 if (jetInd > 1)
0143 FillHist1D("m2jCalo", (p4jet[0] + p4jet[1]).mass());
0144 if (MCarlo) {
0145 evt.getByLabel(genAlgo, genjets);
0146 evt.getByLabel("genEventScale", genEventScale);
0147 pthat = *genEventScale;
0148 FillHist1D("ptHat", pthat);
0149 CaloJet MatchedJet;
0150 jetInd = 0;
0151 if (genjets->empty())
0152 cout << "WARNING: NO gen jets in event " << evt.id().event() << ", Run " << evt.id().run() << " !!!!" << endl;
0153 for (i_genjet = genjets->begin(); i_genjet != genjets->end() && jetInd < Njets; ++i_genjet) {
0154 if (jetInd < 2)
0155 p4jet[jetInd] = i_genjet->p4();
0156 FillHist1D("ptGen", i_genjet->pt());
0157 FillHist1D("etaGen", i_genjet->eta());
0158 FillHist1D("phiGen", i_genjet->phi());
0159 FillHistProfile("GenEnergyVsEta", i_genjet->eta(), i_genjet->energy());
0160 dRmin = 1000.0;
0161 for (i_caljet = caljets->begin(); i_caljet != caljets->end(); ++i_caljet) {
0162 dR = deltaR(i_caljet->eta(), i_caljet->phi(), i_genjet->eta(), i_genjet->phi());
0163 if (dR < dRmin) {
0164 dRmin = dR;
0165 MatchedJet = *i_caljet;
0166 }
0167 }
0168 FillHist1D("dR", dRmin);
0169 e = MatchedJet.energy();
0170 pt = MatchedJet.pt();
0171 eta = MatchedJet.eta();
0172 emEB = MatchedJet.emEnergyInEB();
0173 emEE = MatchedJet.emEnergyInEE();
0174 emHF = MatchedJet.emEnergyInHF();
0175 hadHB = MatchedJet.hadEnergyInHB();
0176 hadHE = MatchedJet.hadEnergyInHE();
0177 hadHO = MatchedJet.hadEnergyInHO();
0178 hadHF = MatchedJet.hadEnergyInHF();
0179 if (dRmin < dRmatch && pt > PtMin) {
0180 FillHistProfile("CaloErespVsEta", eta, e / i_genjet->energy());
0181 FillHistProfile("emErespVsEta", eta, (emEB + emEE + emHF) / i_genjet->energy());
0182 FillHistProfile("hadErespVsEta", eta, (hadHB + hadHO + hadHE + hadHF) / i_genjet->energy());
0183 if (fabs(i_genjet->eta()) < 1.)
0184 FillHistProfile("respVsPtBarrel", i_genjet->pt(), pt / i_genjet->pt());
0185 }
0186 jetInd++;
0187 }
0188 FillHist1D("GenJetMulti", jetInd);
0189 if (jetInd > 1)
0190 FillHist1D("m2jGen", (p4jet[0] + p4jet[1]).mass());
0191 }
0192 }
0193
0194 void JetValidation::endJob() {
0195
0196 if (m_file != nullptr) {
0197 m_file->cd();
0198 for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
0199 hid->second->Write();
0200 for (std::map<TString, TH2*>::iterator hid = m_HistNames2D.begin(); hid != m_HistNames2D.end(); hid++)
0201 hid->second->Write();
0202 for (std::map<TString, TProfile*>::iterator hid = m_HistNamesProfile.begin(); hid != m_HistNamesProfile.end();
0203 hid++)
0204 hid->second->Write();
0205 delete m_file;
0206 m_file = nullptr;
0207 }
0208 }
0209
0210 void JetValidation::FillHist1D(const TString& histName, const Double_t& value) {
0211 std::map<TString, TH1*>::iterator hid = m_HistNames1D.find(histName);
0212 if (hid == m_HistNames1D.end())
0213 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
0214 else
0215 hid->second->Fill(value);
0216 }
0217
0218 void JetValidation::FillHist2D(const TString& histName, const Double_t& valuex, const Double_t& valuey) {
0219 std::map<TString, TH2*>::iterator hid = m_HistNames2D.find(histName);
0220 if (hid == m_HistNames2D.end())
0221 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
0222 else
0223 hid->second->Fill(valuex, valuey);
0224 }
0225
0226 void JetValidation::FillHistProfile(const TString& histName, const Double_t& valuex, const Double_t& valuey) {
0227 std::map<TString, TProfile*>::iterator hid = m_HistNamesProfile.find(histName);
0228 if (hid == m_HistNamesProfile.end())
0229 std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
0230 else
0231 hid->second->Fill(valuex, valuey);
0232 }
0233 #include "FWCore/Framework/interface/MakerMacros.h"
0234 DEFINE_FWK_MODULE(JetValidation);
