JetMET/src/CaloTowerAnalyzer.cc

0001 #include "DQMOffline/JetMET/interface/CaloTowerAnalyzer.h"
0002 // author: Bobby Scurlock, University of Florida
0003 // first version 12/18/2006
0004 // modified: Mike Schmitt
0005 // date: 02.28.2007
0006 // note: code rewrite
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 //#include "Geometry/Vector/interface/GlobalPoint.h"
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   //Store number of events which pass each HLT bit
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   //--Store number of events used
0090   hCT_Nevents = ibooker.book1D("METTask_CT_Nevents", "", 1, 0, 1);
0091   //--Data integrated over all events and stored by CaloTower(ieta,iphi)
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   //--Data over eta-rings
0130
0131   // CaloTower values
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       // Integrated over phi
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         // Integrated over phi
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   }  // allhist
0167 }
0168
0169 void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0170   // Get HLT Results
0171   edm::Handle<edm::TriggerResults> TheHLTResults;
0172   iEvent.getByToken(HLTResultsLabel_, TheHLTResults);
0173
0174   // **** Get the TriggerResults container
0175   //triggerResultsToken_= consumes<edm::TriggerResults>(edm::InputTag(theTriggerResultsLabel));
0176   //edm::Handle<edm::TriggerResults> triggerResults;
0177   //iEvent.getByToken(triggerResultsToken_, triggerResults);
0178
0179   bool EventPasses = true;
0180   // Make sure handle is valid
0181   if (TheHLTResults.isValid() && hltselection_) {
0182     //Get HLT Names
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         //Change the default value since HLT requirement has been issued by the user
0188         if (index == 0)
0189           EventPasses = false;
0190         //Get the HLT bit and check to make sure it is valid
0191         unsigned int bit = TheTriggerNames.triggerIndex(HLTBitLabel_[index].label());
0192         if (bit < TheHLTResults->size()) {
0193           //If any of the HLT names given by the user accept, then the event passes
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   //----------GREG & CHRIS' idea---///
0211   float ETTowerMin = -1;  //GeV
0212   float METRingMin = -2;  // GeV
0213
0214   Nevents++;
0215   hCT_Nevents->Fill(0);
0216
0217   // ==========================================================
0218   // Retrieve!
0219   // ==========================================================
0220
0221   edm::Handle<edm::View<Candidate> > towers;
0222   iEvent.getByToken(caloTowersLabel_, towers);
0223
0224   if ((!towers.isValid())) {
0225     //DD:fix print label
0226     //edm::LogInfo("")<<"CaloTowers "<< caloTowersLabel_<<" not found!"<<std::endl;
0227     return;
0228   }
0229
0230   //HBHENoiseFilterResultToken_=consumes<HBHENoiseFilter>(edm::InputTag(HBHENoiseFilterResultLabel_));
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   // Fill Histograms!
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   //  for (CaloTowerCollection::const_iterator calotower = towers->begin(); calotower != towers->end(); calotower++)
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         //math::RhoEtaPhiVector Momentum = calotower->momentum();
0273         double Tower_ET = calotower->et();
0274         //double Tower_Energy  = calotower->energy();
0275         //    double Tower_Eta = calotower->eta();
0276         double Tower_Phi = calotower->phi();
0277         //double Tower_EMEnergy = calotower->emEnergy();
0278         //double Tower_HadEnergy = calotower->hadEnergy();
0279         double Tower_OuterEt = calotower->outerEt();
0280         double Tower_EMEt = calotower->emEt();
0281         double Tower_HadEt = calotower->hadEt();
0282         //int Tower_EMLV1 = calotower->emLvl1();
0283         //int Tower_HadLV1 = calotower->hadLv11();
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         // Fill Histograms
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       }  //end if (calotower) ..
0330     }  // end if(candidate) ...
0331
0332   }  // end loop over towers
0333
0334   // Fill eta-ring MET quantities
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     }  // ietaring
0351   }  // allhist_
0352
0353   edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
0354   edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
0355 }