File indexing completed on 2025-02-09 23:41:42
0001 #include "DQMOffline/JetMET/interface/CaloTowerAnalyzer.h"
0002
0003
0004
0005
0006
0007
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "FWCore/PluginManager/interface/ModuleDef.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "DataFormats/JetReco/interface/CaloJet.h"
0018 #include "DataFormats/METReco/interface/CaloMET.h"
0019 #include "DataFormats/METReco/interface/GenMET.h"
0020 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0021 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0022 #include "DataFormats/METReco/interface/GenMETCollection.h"
0023 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0024
0025 #include "DataFormats/Math/interface/LorentzVector.h"
0026 #include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"
0027
0028 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0029 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0030 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0031 #include "DataFormats/Common/interface/TriggerResults.h"
0032
0033 #include "FWCore/Common/interface/TriggerNames.h"
0034
0035 #include "DataFormats/DetId/interface/DetId.h"
0036 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0037
0038 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0039 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0040 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0041
0042
0043 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0044 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0045
0046 #include <vector>
0047 #include <utility>
0048 #include <ostream>
0049 #include <fstream>
0050 #include <string>
0051 #include <algorithm>
0052 #include <cmath>
0053 #include <memory>
0054 #include <TLorentzVector.h>
0055
0056 using namespace reco;
0057 using namespace edm;
0058 using namespace std;
0059
0060 CaloTowerAnalyzer::CaloTowerAnalyzer(const edm::ParameterSet& iConfig) {
0061 caloTowersLabel_ = consumes<edm::View<reco::Candidate> >(iConfig.getParameter<edm::InputTag>("CaloTowersLabel"));
0062 HLTResultsLabel_ = consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("HLTResultsLabel"));
0063 HBHENoiseFilterResultLabel_ = consumes<bool>(iConfig.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel"));
0064
0065 if (iConfig.exists("HLTBitLabels"))
0066 HLTBitLabel_ = iConfig.getParameter<std::vector<edm::InputTag> >("HLTBitLabels");
0067
0068 debug_ = iConfig.getParameter<bool>("Debug");
0069 finebinning_ = iConfig.getUntrackedParameter<bool>("FineBinning");
0070 allhist_ = iConfig.getUntrackedParameter<bool>("AllHist");
0071 FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");
0072
0073 hltselection_ = iConfig.getUntrackedParameter<bool>("HLTSelection");
0074 }
0075
0076 void CaloTowerAnalyzer::dqmbeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) { Nevents = 0; }
0077
0078 void CaloTowerAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0079 ibooker.setCurrentFolder(FolderName_);
0080
0081
0082 for (unsigned int i = 0; i < HLTBitLabel_.size(); i++) {
0083 if (!HLTBitLabel_[i].label().empty()) {
0084 hCT_NEvents_HLT.push_back(
0085 ibooker.book1D("METTask_CT_" + HLTBitLabel_[i].label(), HLTBitLabel_[i].label(), 2, -0.5, 1.5));
0086 }
0087 }
0088
0089
0090 hCT_Nevents = ibooker.book1D("METTask_CT_Nevents", "", 1, 0, 1);
0091
0092 hCT_et_ieta_iphi = ibooker.book2D("METTask_CT_et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0093 hCT_et_ieta_iphi->setOption("colz");
0094 hCT_et_ieta_iphi->setAxisTitle("ieta", 1);
0095 hCT_et_ieta_iphi->setAxisTitle("ephi", 2);
0096
0097 hCT_emEt_ieta_iphi = ibooker.book2D("METTask_CT_emEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0098 hCT_emEt_ieta_iphi->setOption("colz");
0099 hCT_emEt_ieta_iphi->setAxisTitle("ieta", 1);
0100 hCT_emEt_ieta_iphi->setAxisTitle("ephi", 2);
0101 hCT_hadEt_ieta_iphi = ibooker.book2D("METTask_CT_hadEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0102 hCT_hadEt_ieta_iphi->setOption("colz");
0103 hCT_hadEt_ieta_iphi->setAxisTitle("ieta", 1);
0104 hCT_hadEt_ieta_iphi->setAxisTitle("ephi", 2);
0105 hCT_outerEt_ieta_iphi = ibooker.book2D("METTask_CT_outerEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0106 hCT_outerEt_ieta_iphi->setOption("colz");
0107 hCT_outerEt_ieta_iphi->setAxisTitle("ieta", 1);
0108 hCT_outerEt_ieta_iphi->setAxisTitle("ephi", 2);
0109 hCT_Occ_ieta_iphi = ibooker.book2D("METTask_CT_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0110 hCT_Occ_ieta_iphi->setOption("colz");
0111 hCT_Occ_ieta_iphi->setAxisTitle("ieta", 1);
0112 hCT_Occ_ieta_iphi->setAxisTitle("ephi", 2);
0113
0114 hCT_Occ_EM_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_EM_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0115 hCT_Occ_EM_Et_ieta_iphi->setOption("colz");
0116 hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ieta", 1);
0117 hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ephi", 2);
0118
0119 hCT_Occ_HAD_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_HAD_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0120 hCT_Occ_HAD_Et_ieta_iphi->setOption("colz");
0121 hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ieta", 1);
0122 hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ephi", 2);
0123
0124 hCT_Occ_Outer_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_Outer_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0125 hCT_Occ_Outer_Et_ieta_iphi->setOption("colz");
0126 hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ieta", 1);
0127 hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ephi", 2);
0128
0129
0130
0131
0132 if (allhist_) {
0133 if (finebinning_) {
0134 hCT_etvsieta = ibooker.book2D("METTask_CT_etvsieta", "", 83, -41, 42, 10001, 0, 1001);
0135 hCT_Minetvsieta = ibooker.book2D("METTask_CT_Minetvsieta", "", 83, -41, 42, 10001, 0, 1001);
0136 hCT_Maxetvsieta = ibooker.book2D("METTask_CT_Maxetvsieta", "", 83, -41, 42, 10001, 0, 1001);
0137 hCT_emEtvsieta = ibooker.book2D("METTask_CT_emEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
0138 hCT_hadEtvsieta = ibooker.book2D("METTask_CT_hadEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
0139 hCT_outerEtvsieta = ibooker.book2D("METTask_CT_outerEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
0140
0141
0142 hCT_Occvsieta = ibooker.book2D("METTask_CT_Occvsieta", "", 83, -41, 42, 84, 0, 84);
0143 hCT_SETvsieta = ibooker.book2D("METTask_CT_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
0144 hCT_METvsieta = ibooker.book2D("METTask_CT_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
0145 hCT_METPhivsieta = ibooker.book2D("METTask_CT_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0146 hCT_MExvsieta = ibooker.book2D("METTask_CT_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
0147 hCT_MEyvsieta = ibooker.book2D("METTask_CT_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
0148 } else {
0149 if (allhist_) {
0150 hCT_etvsieta = ibooker.book2D("METTask_CT_etvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
0151 hCT_Minetvsieta = ibooker.book2D("METTask_CT_Minetvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
0152 hCT_Maxetvsieta = ibooker.book2D("METTask_CT_Maxetvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
0153 hCT_emEtvsieta = ibooker.book2D("METTask_CT_emEtvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
0154 hCT_hadEtvsieta = ibooker.book2D("METTask_CT_hadEtvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
0155 hCT_outerEtvsieta = ibooker.book2D("METTask_CT_outerEtvsieta", "", 83, -41, 42, 80, -0.5, 399.5);
0156
0157 }
0158
0159 hCT_Occvsieta = ibooker.book2D("METTask_CT_Occvsieta", "", 83, -41, 42, 73, -0.5, 72.5);
0160 hCT_SETvsieta = ibooker.book2D("METTask_CT_SETvsieta", "", 83, -41, 42, 200, -0.5, 1999.5);
0161 hCT_METvsieta = ibooker.book2D("METTask_CT_METvsieta", "", 83, -41, 42, 200, -0.5, 1999.5);
0162 hCT_METPhivsieta = ibooker.book2D("METTask_CT_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0163 hCT_MExvsieta = ibooker.book2D("METTask_CT_MExvsieta", "", 83, -41, 42, 100, -499.5, 499.5);
0164 hCT_MEyvsieta = ibooker.book2D("METTask_CT_MEyvsieta", "", 83, -41, 42, 100, -499.5, 499.5);
0165 }
0166 }
0167 }
0168
0169 void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0170
0171 edm::Handle<edm::TriggerResults> TheHLTResults;
0172 iEvent.getByToken(HLTResultsLabel_, TheHLTResults);
0173
0174
0175
0176
0177
0178
0179 bool EventPasses = true;
0180
0181 if (TheHLTResults.isValid() && hltselection_) {
0182
0183 const edm::TriggerNames& TheTriggerNames = iEvent.triggerNames(*TheHLTResults);
0184
0185 for (unsigned int index = 0; index < HLTBitLabel_.size(); index++) {
0186 if (!HLTBitLabel_[index].label().empty()) {
0187
0188 if (index == 0)
0189 EventPasses = false;
0190
0191 unsigned int bit = TheTriggerNames.triggerIndex(HLTBitLabel_[index].label());
0192 if (bit < TheHLTResults->size()) {
0193
0194 if (TheHLTResults->accept(bit) && !TheHLTResults->error(bit)) {
0195 EventPasses = true;
0196 hCT_NEvents_HLT[index]->Fill(1);
0197 } else
0198 hCT_NEvents_HLT[index]->Fill(0);
0199 } else {
0200 edm::LogInfo("OutputInfo") << "The HLT Trigger Name : " << HLTBitLabel_[index].label()
0201 << " is not valid for this trigger table " << std::endl;
0202 }
0203 }
0204 }
0205 }
0206
0207 if (!EventPasses && hltselection_)
0208 return;
0209
0210
0211 float ETTowerMin = -1;
0212 float METRingMin = -2;
0213
0214 Nevents++;
0215 hCT_Nevents->Fill(0);
0216
0217
0218
0219
0220
0221 edm::Handle<edm::View<Candidate> > towers;
0222 iEvent.getByToken(caloTowersLabel_, towers);
0223
0224 if ((!towers.isValid())) {
0225
0226
0227 return;
0228 }
0229
0230
0231 edm::Handle<bool> HBHENoiseFilterResultHandle;
0232 iEvent.getByToken(HBHENoiseFilterResultLabel_, HBHENoiseFilterResultHandle);
0233 bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
0234 if (!HBHENoiseFilterResultHandle.isValid()) {
0235 LogDebug("") << "CaloTowerAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
0236 }
0237
0238 bool bHcalNoiseFilter = HBHENoiseFilterResult;
0239
0240 if (!bHcalNoiseFilter)
0241 return;
0242
0243 edm::View<Candidate>::const_iterator towerCand = towers->begin();
0244
0245
0246
0247
0248
0249 int CTmin_iphi = 99, CTmax_iphi = -99;
0250 int CTmin_ieta = 99, CTmax_ieta = -99;
0251
0252 TLorentzVector vMET_EtaRing[83];
0253 int ActiveRing[83];
0254 int NActiveTowers[83];
0255 double SET_EtaRing[83];
0256 double MinEt_EtaRing[83];
0257 double MaxEt_EtaRing[83];
0258 for (int i = 0; i < 83; i++) {
0259 ActiveRing[i] = 0;
0260 NActiveTowers[i] = 0;
0261 SET_EtaRing[i] = 0;
0262 MinEt_EtaRing[i] = 0;
0263 MaxEt_EtaRing[i] = 0;
0264 }
0265
0266
0267 for (; towerCand != towers->end(); towerCand++) {
0268 const Candidate* candidate = &(*towerCand);
0269 if (candidate) {
0270 const CaloTower* calotower = dynamic_cast<const CaloTower*>(candidate);
0271 if (calotower) {
0272
0273 double Tower_ET = calotower->et();
0274
0275
0276 double Tower_Phi = calotower->phi();
0277
0278
0279 double Tower_OuterEt = calotower->outerEt();
0280 double Tower_EMEt = calotower->emEt();
0281 double Tower_HadEt = calotower->hadEt();
0282
0283
0284 int Tower_ieta = calotower->id().ieta();
0285 int Tower_iphi = calotower->id().iphi();
0286 int EtaRing = 41 + Tower_ieta;
0287 ActiveRing[EtaRing] = 1;
0288 NActiveTowers[EtaRing]++;
0289 SET_EtaRing[EtaRing] += Tower_ET;
0290 TLorentzVector v_;
0291 v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
0292 if (Tower_ET > ETTowerMin)
0293 vMET_EtaRing[EtaRing] -= v_;
0294
0295
0296 hCT_Occ_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
0297 if (calotower->emEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
0298 hCT_Occ_EM_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
0299 if (calotower->hadEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
0300 hCT_Occ_HAD_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
0301 if (calotower->outerEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
0302 hCT_Occ_Outer_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
0303
0304 hCT_et_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_ET);
0305 hCT_emEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_EMEt);
0306 hCT_hadEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_HadEt);
0307 hCT_outerEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_OuterEt);
0308
0309 if (allhist_) {
0310 hCT_etvsieta->Fill(Tower_ieta, Tower_ET);
0311 hCT_emEtvsieta->Fill(Tower_ieta, Tower_EMEt);
0312 hCT_hadEtvsieta->Fill(Tower_ieta, Tower_HadEt);
0313 hCT_outerEtvsieta->Fill(Tower_ieta, Tower_OuterEt);
0314 }
0315
0316 if (Tower_ET > MaxEt_EtaRing[EtaRing])
0317 MaxEt_EtaRing[EtaRing] = Tower_ET;
0318 if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET > 0)
0319 MinEt_EtaRing[EtaRing] = Tower_ET;
0320
0321 if (Tower_ieta < CTmin_ieta)
0322 CTmin_ieta = Tower_ieta;
0323 if (Tower_ieta > CTmax_ieta)
0324 CTmax_ieta = Tower_ieta;
0325 if (Tower_iphi < CTmin_iphi)
0326 CTmin_iphi = Tower_iphi;
0327 if (Tower_iphi > CTmax_iphi)
0328 CTmax_iphi = Tower_iphi;
0329 }
0330 }
0331
0332 }
0333
0334
0335 if (allhist_) {
0336 for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0337 hCT_Minetvsieta->Fill(iEtaRing - 41, MinEt_EtaRing[iEtaRing]);
0338 hCT_Maxetvsieta->Fill(iEtaRing - 41, MaxEt_EtaRing[iEtaRing]);
0339
0340 if (ActiveRing[iEtaRing]) {
0341 if (vMET_EtaRing[iEtaRing].Pt() > METRingMin) {
0342 hCT_METPhivsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Phi());
0343 hCT_MExvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Px());
0344 hCT_MEyvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Py());
0345 hCT_METvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Pt());
0346 }
0347 hCT_SETvsieta->Fill(iEtaRing - 41, SET_EtaRing[iEtaRing]);
0348 hCT_Occvsieta->Fill(iEtaRing - 41, NActiveTowers[iEtaRing]);
0349 }
0350 }
0351 }
0352
0353 edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
0354 edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
0355 }