Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:56

0001 // JetValidation.cc

0002 // Description:  Some Basic validation plots for jets.

0003 // Author: K. Kousouris

0004 // Date:  27 - August - 2008

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   /////////// Write Histograms in output ROOT file ////////

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);