CMS 3D CMS Logo
JetValidation.cc
Go to the documentation of this file.
1 // JetValidation.cc
2 // Description: Some Basic validation plots for jets.
3 // Author: K. Kousouris
4 // Date: 27 - August - 2008
5 //
6 #include "RecoJets/JetAnalyzers/interface/JetValidation.h"
7 #include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
8 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
9 #include "DataFormats/JetReco/interface/GenJetCollection.h"
10 #include "DataFormats/JetReco/interface/CaloJet.h"
11 #include "DataFormats/JetReco/interface/GenJet.h"
12 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
13 #include "DataFormats/Math/interface/deltaR.h"
14 #include "FWCore/Framework/interface/Event.h"
15 #include "FWCore/ParameterSet/interface/ParameterSet.h"
16 #include <TROOT.h>
17 #include <TSystem.h>
18 #include <TFile.h>
19 #include <cmath>
20 using namespace edm;
21 using namespace reco;
22 using namespace std;
24 JetValidation::JetValidation(edm::ParameterSet const& cfg) {
25  dRmatch = cfg.getParameter<double>("dRmatch");
26  PtMin = cfg.getParameter<double>("PtMin");
27  Njets = cfg.getParameter<int>("Njets");
28  MCarlo = cfg.getParameter<bool>("MCarlo");
29  genAlgo = cfg.getParameter<string>("genAlgo");
30  calAlgo = cfg.getParameter<string>("calAlgo");
31  jetTracksAssociator = cfg.getParameter<string>("jetTracksAssociator");
32  histoFileName = cfg.getParameter<string>("histoFileName");
33 }
35 void JetValidation::beginJob() {
36  m_file = new TFile(histoFileName.c_str(), "RECREATE");
37 
38  m_HistNames1D["CaloJetMulti"] = new TH1F("CaloJetMulti", "Multiplicity of CaloJets", 100, 0, 100);
39  m_HistNames1D["ptCalo"] = new TH1F("ptCalo", "p_{T} of CaloJets", 7000, 0, 7000);
40  m_HistNames1D["etaCalo"] = new TH1F("etaCalo", "#eta of CaloJets", 100, -5.0, 5.0);
41  m_HistNames1D["phiCalo"] = new TH1F("phiCalo", "#phi of CaloJets", 72, -M_PI, M_PI);
42  m_HistNames1D["m2jCalo"] = new TH1F("m2jCalo", "Dijet Mass of leading CaloJets", 7000, 0, 14000);
43  m_HistNames1D["nTracks"] = new TH1F("nTracks", "Number of tracks associated with a jet", 100, 0, 100);
44  m_HistNames1D["chargeFraction"] = new TH1F("chargeFraction", "Fraction of charged tracks pt", 500, 0, 5);
45  m_HistNames1D["emEnergyFraction"] = new TH1F("emEnergyFraction", "Jets EM Fraction", 110, 0, 1.1);
46  m_HistNames1D["emEnergyInEB"] = new TH1F("emEnergyInEB", "Jets emEnergyInEB", 7000, 0, 14000);
47  m_HistNames1D["emEnergyInEE"] = new TH1F("emEnergyInEE", "Jets emEnergyInEE", 7000, 0, 14000);
48  m_HistNames1D["emEnergyInHF"] = new TH1F("emEnergyInHF", "Jets emEnergyInHF", 7000, 0, 14000);
49  m_HistNames1D["hadEnergyInHB"] = new TH1F("hadEnergyInHB", "Jets hadEnergyInHB", 7000, 0, 14000);
50  m_HistNames1D["hadEnergyInHE"] = new TH1F("hadEnergyInHE", "Jets hadEnergyInHE", 7000, 0, 14000);
51  m_HistNames1D["hadEnergyInHF"] = new TH1F("hadEnergyInHF", "Jets hadEnergyInHF", 7000, 0, 14000);
52  m_HistNames1D["hadEnergyInHO"] = new TH1F("hadEnergyInHO", "Jets hadEnergyInHO", 7000, 0, 14000);
53  m_HistNamesProfile["EBfractionVsEta"] = new TProfile("EBfractionVsEta", "Jets EBfraction vs #eta", 100, -5.0, 5.0);
54  m_HistNamesProfile["EEfractionVsEta"] = new TProfile("EEfractionVsEta", "Jets EEfraction vs #eta", 100, -5.0, 5.0);
55  m_HistNamesProfile["HBfractionVsEta"] = new TProfile("HBfractionVsEta", "Jets HBfraction vs #eta", 100, -5.0, 5.0);
56  m_HistNamesProfile["HOfractionVsEta"] = new TProfile("HOfractionVsEta", "Jets HOfraction vs #eta", 100, -5.0, 5.0);
57  m_HistNamesProfile["HEfractionVsEta"] = new TProfile("HEfractionVsEta", "Jets HEfraction vs #eta", 100, -5.0, 5.0);
58  m_HistNamesProfile["HFfractionVsEta"] = new TProfile("HFfractionVsEta", "Jets HFfraction vs #eta", 100, -5.0, 5.0);
59  m_HistNamesProfile["CaloEnergyVsEta"] = new TProfile("CaloEnergyVsEta", "CaloJets Energy Vs. Eta", 100, -5.0, 5.0);
60  m_HistNamesProfile["emEnergyVsEta"] = new TProfile("emEnergyVsEta", "Jets EM Energy Vs. Eta", 100, -5.0, 5.0);
61  m_HistNamesProfile["hadEnergyVsEta"] = new TProfile("hadEnergyVsEta", "Jets HAD Energy Vs. Eta", 100, -5.0, 5.0);
62  if (MCarlo) {
63  m_HistNames1D["GenJetMulti"] = new TH1F("GenJetMulti", "Multiplicity of GenJets", 100, 0, 100);
64  m_HistNames1D["ptHat"] = new TH1F("ptHat", "p_{T}hat", 7000, 0, 7000);
65  m_HistNames1D["ptGen"] = new TH1F("ptGen", "p_{T} of GenJets", 7000, 0, 7000);
66  m_HistNames1D["etaGen"] = new TH1F("etaGen", "#eta of GenJets", 100, -5.0, 5.0);
67  m_HistNames1D["phiGen"] = new TH1F("phiGen", "#phi of GenJets", 72, -M_PI, M_PI);
68  m_HistNames1D["m2jGen"] = new TH1F("m2jGen", "Dijet Mass of leading GenJets", 7000, 0, 14000);
69  m_HistNames1D["dR"] = new TH1F("dR", "GenJets dR with matched CaloJet", 200, 0, 1);
70  m_HistNamesProfile["GenEnergyVsEta"] = new TProfile("GenEnergyVsEta", "GenJets Energy Vs. Eta", 100, -5.0, 5.0);
71  m_HistNamesProfile["respVsPtBarrel"] =
72  new TProfile("respVsPtBarrel", "CaloJet Response of GenJets in Barrel", 7000, 0, 7000);
73  m_HistNamesProfile["CaloErespVsEta"] =
74  new TProfile("CaloErespVsEta", "Jets Energy Response Vs. Eta", 100, -5.0, 5.0);
75  m_HistNamesProfile["emErespVsEta"] =
76  new TProfile("emErespVsEta", "Jets EM Energy Response Vs. Eta", 100, -5.0, 5.0);
77  m_HistNamesProfile["hadErespVsEta"] =
78  new TProfile("hadErespVsEta", "Jets HAD Energy Response Vs. Eta", 100, -5.0, 5.0);
79  }
80 }
82 void JetValidation::analyze(edm::Event const& evt, edm::EventSetup const& iSetup) {
83  math::XYZTLorentzVector p4jet[2];
84  int jetInd, jetCounter, nTracks;
85  double dRmin, dR, e, eta, emEB, emEE, emHF, hadHB, hadHE, hadHO, hadHF, pt, phi, pthat, chf;
86  Handle<CaloJetCollection> caljets;
87  Handle<GenJetCollection> genjets;
88  Handle<double> genEventScale;
89  Handle<JetTracksAssociation::Container> jetTracks;
90  CaloJetCollection::const_iterator i_caljet;
91  GenJetCollection::const_iterator i_genjet;
92  evt.getByLabel(calAlgo, caljets);
93  evt.getByLabel(jetTracksAssociator, jetTracks);
94  jetInd = 0;
95  jetCounter = 0;
96  if (caljets->empty())
97  cout << "WARNING: NO calo jets in event " << evt.id().event() << ", Run " << evt.id().run() << " !!!!" << endl;
98  for (i_caljet = caljets->begin(); i_caljet != caljets->end() && jetInd < Njets; ++i_caljet) {
99  e = i_caljet->energy();
100  pt = i_caljet->pt();
101  phi = i_caljet->phi();
102  eta = i_caljet->eta();
103  emEB = i_caljet->emEnergyInEB();
104  emEE = i_caljet->emEnergyInEE();
105  emHF = i_caljet->emEnergyInHF();
106  hadHB = i_caljet->hadEnergyInHB();
107  hadHE = i_caljet->hadEnergyInHE();
108  hadHO = i_caljet->hadEnergyInHO();
109  hadHF = i_caljet->hadEnergyInHF();
110  nTracks = JetTracksAssociation::tracksNumber(*jetTracks, *i_caljet);
111  chf = (JetTracksAssociation::tracksP4(*jetTracks, *i_caljet)).pt() / pt;
112  if (jetInd < 2)
113  p4jet[jetInd] = i_caljet->p4();
114  if (pt > PtMin) {
115  FillHist1D("ptCalo", pt);
116  FillHist1D("etaCalo", eta);
117  FillHist1D("phiCalo", phi);
118  FillHist1D("emEnergyFraction", i_caljet->emEnergyFraction());
119  FillHist1D("nTracks", nTracks);
120  FillHist1D("chargeFraction", chf);
121  FillHist1D("emEnergyInEB", emEB);
122  FillHist1D("emEnergyInEE", emEE);
123  FillHist1D("emEnergyInHF", emHF);
124  FillHist1D("hadEnergyInHB", hadHB);
125  FillHist1D("hadEnergyInHE", hadHE);
126  FillHist1D("hadEnergyInHF", hadHF);
127  FillHist1D("hadEnergyInHO", hadHO);
128  FillHistProfile("EBfractionVsEta", eta, emEB / e);
129  FillHistProfile("EEfractionVsEta", eta, emEE / e);
130  FillHistProfile("HBfractionVsEta", eta, hadHB / e);
131  FillHistProfile("HOfractionVsEta", eta, hadHO / e);
132  FillHistProfile("HEfractionVsEta", eta, hadHE / e);
133  FillHistProfile("HFfractionVsEta", eta, (hadHF + emHF) / e);
134  FillHistProfile("CaloEnergyVsEta", eta, e);
135  FillHistProfile("emEnergyVsEta", eta, emEB + emEE + emHF);
136  FillHistProfile("hadEnergyVsEta", eta, hadHB + hadHO + hadHE + hadHF);
137  jetCounter++;
138  }
139  jetInd++;
140  }
141  FillHist1D("CaloJetMulti", jetCounter);
142  if (jetInd > 1)
143  FillHist1D("m2jCalo", (p4jet[0] + p4jet[1]).mass());
144  if (MCarlo) {
145  evt.getByLabel(genAlgo, genjets);
146  evt.getByLabel("genEventScale", genEventScale);
147  pthat = *genEventScale;
148  FillHist1D("ptHat", pthat);
149  CaloJet MatchedJet;
150  jetInd = 0;
151  if (genjets->empty())
152  cout << "WARNING: NO gen jets in event " << evt.id().event() << ", Run " << evt.id().run() << " !!!!" << endl;
153  for (i_genjet = genjets->begin(); i_genjet != genjets->end() && jetInd < Njets; ++i_genjet) {
154  if (jetInd < 2)
155  p4jet[jetInd] = i_genjet->p4();
156  FillHist1D("ptGen", i_genjet->pt());
157  FillHist1D("etaGen", i_genjet->eta());
158  FillHist1D("phiGen", i_genjet->phi());
159  FillHistProfile("GenEnergyVsEta", i_genjet->eta(), i_genjet->energy());
160  dRmin = 1000.0;
161  for (i_caljet = caljets->begin(); i_caljet != caljets->end(); ++i_caljet) {
162  dR = deltaR(i_caljet->eta(), i_caljet->phi(), i_genjet->eta(), i_genjet->phi());
163  if (dR < dRmin) {
164  dRmin = dR;
165  MatchedJet = *i_caljet;
166  }
167  }
168  FillHist1D("dR", dRmin);
169  e = MatchedJet.energy();
170  pt = MatchedJet.pt();
171  eta = MatchedJet.eta();
172  emEB = MatchedJet.emEnergyInEB();
173  emEE = MatchedJet.emEnergyInEE();
174  emHF = MatchedJet.emEnergyInHF();
175  hadHB = MatchedJet.hadEnergyInHB();
176  hadHE = MatchedJet.hadEnergyInHE();
177  hadHO = MatchedJet.hadEnergyInHO();
178  hadHF = MatchedJet.hadEnergyInHF();
179  if (dRmin < dRmatch && pt > PtMin) {
180  FillHistProfile("CaloErespVsEta", eta, e / i_genjet->energy());
181  FillHistProfile("emErespVsEta", eta, (emEB + emEE + emHF) / i_genjet->energy());
182  FillHistProfile("hadErespVsEta", eta, (hadHB + hadHO + hadHE + hadHF) / i_genjet->energy());
183  if (fabs(i_genjet->eta()) < 1.)
184  FillHistProfile("respVsPtBarrel", i_genjet->pt(), pt / i_genjet->pt());
185  }
186  jetInd++;
187  }
188  FillHist1D("GenJetMulti", jetInd);
189  if (jetInd > 1)
190  FillHist1D("m2jGen", (p4jet[0] + p4jet[1]).mass());
191  }
192 }
194 void JetValidation::endJob() {
196  if (m_file != nullptr) {
197  m_file->cd();
198  for (std::map<TString, TH1*>::iterator hid = m_HistNames1D.begin(); hid != m_HistNames1D.end(); hid++)
199  hid->second->Write();
200  for (std::map<TString, TH2*>::iterator hid = m_HistNames2D.begin(); hid != m_HistNames2D.end(); hid++)
201  hid->second->Write();
202  for (std::map<TString, TProfile*>::iterator hid = m_HistNamesProfile.begin(); hid != m_HistNamesProfile.end();
203  hid++)
204  hid->second->Write();
205  delete m_file;
206  m_file = nullptr;
207  }
208 }
210 void JetValidation::FillHist1D(const TString& histName, const Double_t& value) {
211  std::map<TString, TH1*>::iterator hid = m_HistNames1D.find(histName);
212  if (hid == m_HistNames1D.end())
213  std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
214  else
215  hid->second->Fill(value);
216 }
218 void JetValidation::FillHist2D(const TString& histName, const Double_t& valuex, const Double_t& valuey) {
219  std::map<TString, TH2*>::iterator hid = m_HistNames2D.find(histName);
220  if (hid == m_HistNames2D.end())
221  std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
222  else
223  hid->second->Fill(valuex, valuey);
224 }
226 void JetValidation::FillHistProfile(const TString& histName, const Double_t& valuex, const Double_t& valuey) {
227  std::map<TString, TProfile*>::iterator hid = m_HistNamesProfile.find(histName);
228  if (hid == m_HistNamesProfile.end())
229  std::cout << "%fillHist -- Could not find histogram with name: " << histName << std::endl;
230  else
231  hid->second->Fill(valuex, valuey);
232 }
233 #include "FWCore/Framework/interface/MakerMacros.h"
234 DEFINE_FWK_MODULE(JetValidation);
JetValidation::FillHist1D
void FillHist1D(const TString &histName, const Double_t &x)
Definition: JetValidation.cc:210
JetTracksAssociation.h
CaloJetCollection.h
GenJetCollection.h
reco::CaloJet
Jets made from CaloTowers.
Definition: CaloJet.h:27
CaloJet.h
DiDispStaMuonMonitor_cfi.pt
pt
Definition: DiDispStaMuonMonitor_cfi.py:39
edm
HLT enums.
Definition: AlignableModifier.h:19
gather_cfg.cout
cout
Definition: gather_cfg.py:144
reco::JetTracksAssociation::tracksP4
LorentzVector tracksP4(const Container &, const reco::JetBaseRef)
Get LorentzVector as sum of all tracks associated with jet.
Definition: JetTracksAssociation.cc:14
JetValidation.h
reco::CaloJet::hadEnergyInHF
float hadEnergyInHF() const
Definition: CaloJet.h:105
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
reco::CaloJet::hadEnergyInHB
float hadEnergyInHB() const
Definition: CaloJet.h:99
BeamSpotPI::nTracks
Definition: BeamSpotPayloadInspectorHelper.h:42
edm::Handle
Definition: AssociativeIterator.h:50
reco::CaloJet::emEnergyInEE
float emEnergyInEE() const
Definition: CaloJet.h:109
MakerMacros.h
mix_2012_lumiLevel_15_20_50ns_PoissonOOTPU_cfi.histoFileName
histoFileName
Definition: mix_2012_lumiLevel_15_20_50ns_PoissonOOTPU_cfi.py:52
reco::CaloJet::hadEnergyInHE
float hadEnergyInHE() const
Definition: CaloJet.h:103
JetAlgoHelper.h
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
HLT_FULL_cff.PtMin
PtMin
Definition: HLT_FULL_cff.py:9904
PVValHelper::eta
Definition: PVValidationHelpers.h:70
ECF_cff.Njets
Njets
Definition: ECF_cff.py:21
reco::CaloJet::emEnergyInEB
float emEnergyInEB() const
Definition: CaloJet.h:107
JetValidation::JetValidation
JetValidation(edm::ParameterSet const &cfg)
Definition: JetValidation.cc:24
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:496
reco::CaloJet::hadEnergyInHO
float hadEnergyInHO() const
Definition: CaloJet.h:101
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:38
PbPb_ZMuSkimMuonDPG_cff.deltaR
deltaR
Definition: PbPb_ZMuSkimMuonDPG_cff.py:63
edm::ParameterSet
Definition: ParameterSet.h:47
JetValidation::beginJob
void beginJob() override
Definition: JetValidation.cc:35
Event.h
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
deltaR.h
JetValidation::FillHistProfile
void FillHistProfile(const TString &histName, const Double_t &x, const Double_t &y)
Definition: JetValidation.cc:226
HLT_FULL_cff.jetTracks
jetTracks
Definition: HLT_FULL_cff.py:50870
JetValidation::endJob
void endJob() override
Definition: JetValidation.cc:194
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:40
value
Definition: value.py:1
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
edm::EventSetup
Definition: EventSetup.h:58
JetValidation
Definition: JetValidation.h:17
looper.cfg
cfg
Definition: looper.py:297
std
Definition: JetResolutionObject.h:76
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
EgHLTOffHistBins_cfi.mass
mass
Definition: EgHLTOffHistBins_cfi.py:34
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:125
JetValidation::analyze
void analyze(edm::Event const &e, edm::EventSetup const &iSetup) override
Definition: JetValidation.cc:82
GenJet.h
HGC3DClusterGenMatchSelector_cfi.dR
dR
Definition: HGC3DClusterGenMatchSelector_cfi.py:7
ParameterSet.h
reco::JetTracksAssociation::tracksNumber
int tracksNumber(const Container &, const reco::JetBaseRef)
Get number of tracks associated with jet.
Definition: JetTracksAssociation.cc:7
edm::Event
Definition: Event.h:73
reco::CaloJet::emEnergyInHF
float emEnergyInHF() const
Definition: CaloJet.h:111
JetValidation::FillHist2D
void FillHist2D(const TString &histName, const Double_t &x, const Double_t &y)
Definition: JetValidation.cc:218
MillePedeFileConverter_cfg.e
e
Definition: MillePedeFileConverter_cfg.py:37
Generated for CMSSW Reference Manual by doxygen 1.8.16