CaloTowerAnalyzer.cc

Go to the documentation of this file.

#include "DQMOffline/JetMET/interface/CaloTowerAnalyzer.h"
// author: Bobby Scurlock, University of Florida
// first version 12/18/2006
// modified: Mike Schmitt
// date: 02.28.2007
// note: code rewrite

#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/METReco/interface/CaloMET.h"
#include "DataFormats/METReco/interface/GenMET.h"
#include "DataFormats/METReco/interface/CaloMETCollection.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/METReco/interface/GenMETCollection.h"
#include "DataFormats/METReco/interface/CaloMETCollection.h"

#include "DataFormats/Math/interface/LorentzVector.h"
#include "RecoMET/METAlgorithms/interface/CaloSpecificAlgo.h"

#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/Common/interface/TriggerResults.h"

#include "FWCore/Common/interface/TriggerNames.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"

#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"

//#include "Geometry/Vector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include <vector>
#include <utility>
#include <ostream>
#include <fstream>
#include <string>
#include <algorithm>
#include <cmath>
#include <memory>
#include <TLorentzVector.h>

using namespace reco;
using namespace edm;
using namespace std;

CaloTowerAnalyzer::CaloTowerAnalyzer(const edm::ParameterSet& iConfig) {
  caloTowersLabel_ = consumes<edm::View<reco::Candidate> >(iConfig.getParameter<edm::InputTag>("CaloTowersLabel"));
  HLTResultsLabel_ = consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("HLTResultsLabel"));
  HBHENoiseFilterResultLabel_ = consumes<bool>(iConfig.getParameter<edm::InputTag>("HBHENoiseFilterResultLabel"));

  if (iConfig.exists("HLTBitLabels"))
    HLTBitLabel_ = iConfig.getParameter<std::vector<edm::InputTag> >("HLTBitLabels");

  debug_ = iConfig.getParameter<bool>("Debug");
  finebinning_ = iConfig.getUntrackedParameter<bool>("FineBinning");
  allhist_ = iConfig.getUntrackedParameter<bool>("AllHist");
  FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");

  hltselection_ = iConfig.getUntrackedParameter<bool>("HLTSelection");
}

void CaloTowerAnalyzer::dqmbeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) { Nevents = 0; }

void CaloTowerAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
  ibooker.setCurrentFolder(FolderName_);

  //Store number of events which pass each HLT bit
  for (unsigned int i = 0; i < HLTBitLabel_.size(); i++) {
    if (!HLTBitLabel_[i].label().empty()) {
      hCT_NEvents_HLT.push_back(
          ibooker.book1D("METTask_CT_" + HLTBitLabel_[i].label(), HLTBitLabel_[i].label(), 2, -0.5, 1.5));
    }
  }

  //--Store number of events used
  hCT_Nevents = ibooker.book1D("METTask_CT_Nevents", "", 1, 0, 1);
  //--Data integrated over all events and stored by CaloTower(ieta,iphi)
  hCT_et_ieta_iphi = ibooker.book2D("METTask_CT_et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_et_ieta_iphi->setOption("colz");
  hCT_et_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_et_ieta_iphi->setAxisTitle("ephi", 2);

  hCT_emEt_ieta_iphi = ibooker.book2D("METTask_CT_emEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_emEt_ieta_iphi->setOption("colz");
  hCT_emEt_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_emEt_ieta_iphi->setAxisTitle("ephi", 2);
  hCT_hadEt_ieta_iphi = ibooker.book2D("METTask_CT_hadEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_hadEt_ieta_iphi->setOption("colz");
  hCT_hadEt_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_hadEt_ieta_iphi->setAxisTitle("ephi", 2);
  hCT_outerEt_ieta_iphi = ibooker.book2D("METTask_CT_outerEt_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_outerEt_ieta_iphi->setOption("colz");
  hCT_outerEt_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_outerEt_ieta_iphi->setAxisTitle("ephi", 2);
  hCT_Occ_ieta_iphi = ibooker.book2D("METTask_CT_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_Occ_ieta_iphi->setOption("colz");
  hCT_Occ_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_Occ_ieta_iphi->setAxisTitle("ephi", 2);

  hCT_Occ_EM_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_EM_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_Occ_EM_Et_ieta_iphi->setOption("colz");
  hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_Occ_EM_Et_ieta_iphi->setAxisTitle("ephi", 2);

  hCT_Occ_HAD_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_HAD_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_Occ_HAD_Et_ieta_iphi->setOption("colz");
  hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_Occ_HAD_Et_ieta_iphi->setAxisTitle("ephi", 2);

  hCT_Occ_Outer_Et_ieta_iphi = ibooker.book2D("METTask_CT_Occ_Outer_Et_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
  hCT_Occ_Outer_Et_ieta_iphi->setOption("colz");
  hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ieta", 1);
  hCT_Occ_Outer_Et_ieta_iphi->setAxisTitle("ephi", 2);

  //--Data over eta-rings

  // CaloTower values
  if (allhist_) {
    if (finebinning_) {
      hCT_etvsieta = ibooker.book2D("METTask_CT_etvsieta", "", 83, -41, 42, 10001, 0, 1001);
      hCT_Minetvsieta = ibooker.book2D("METTask_CT_Minetvsieta", "", 83, -41, 42, 10001, 0, 1001);
      hCT_Maxetvsieta = ibooker.book2D("METTask_CT_Maxetvsieta", "", 83, -41, 42, 10001, 0, 1001);
      hCT_emEtvsieta = ibooker.book2D("METTask_CT_emEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
      hCT_hadEtvsieta = ibooker.book2D("METTask_CT_hadEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
      hCT_outerEtvsieta = ibooker.book2D("METTask_CT_outerEtvsieta", "", 83, -41, 42, 10001, 0, 1001);
      // Integrated over phi

      hCT_Occvsieta = ibooker.book2D("METTask_CT_Occvsieta", "", 83, -41, 42, 84, 0, 84);
      hCT_SETvsieta = ibooker.book2D("METTask_CT_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
      hCT_METvsieta = ibooker.book2D("METTask_CT_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
      hCT_METPhivsieta = ibooker.book2D("METTask_CT_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
      hCT_MExvsieta = ibooker.book2D("METTask_CT_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
      hCT_MEyvsieta = ibooker.book2D("METTask_CT_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
    } else {
      if (allhist_) {
        hCT_etvsieta = ibooker.book2D("METTask_CT_etvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
        hCT_Minetvsieta = ibooker.book2D("METTask_CT_Minetvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
        hCT_Maxetvsieta = ibooker.book2D("METTask_CT_Maxetvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
        hCT_emEtvsieta = ibooker.book2D("METTask_CT_emEtvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
        hCT_hadEtvsieta = ibooker.book2D("METTask_CT_hadEtvsieta", "", 83, -41, 42, 200, -0.5, 999.5);
        hCT_outerEtvsieta = ibooker.book2D("METTask_CT_outerEtvsieta", "", 83, -41, 42, 80, -0.5, 399.5);
        // Integrated over phi
      }

      hCT_Occvsieta = ibooker.book2D("METTask_CT_Occvsieta", "", 83, -41, 42, 73, -0.5, 72.5);
      hCT_SETvsieta = ibooker.book2D("METTask_CT_SETvsieta", "", 83, -41, 42, 200, -0.5, 1999.5);
      hCT_METvsieta = ibooker.book2D("METTask_CT_METvsieta", "", 83, -41, 42, 200, -0.5, 1999.5);
      hCT_METPhivsieta = ibooker.book2D("METTask_CT_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
      hCT_MExvsieta = ibooker.book2D("METTask_CT_MExvsieta", "", 83, -41, 42, 100, -499.5, 499.5);
      hCT_MEyvsieta = ibooker.book2D("METTask_CT_MEyvsieta", "", 83, -41, 42, 100, -499.5, 499.5);
    }
  }  // allhist
}

void CaloTowerAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // Get HLT Results
  edm::Handle<edm::TriggerResults> TheHLTResults;
  iEvent.getByToken(HLTResultsLabel_, TheHLTResults);

  // **** Get the TriggerResults container
  //triggerResultsToken_= consumes<edm::TriggerResults>(edm::InputTag(theTriggerResultsLabel));
  //edm::Handle<edm::TriggerResults> triggerResults;
  //iEvent.getByToken(triggerResultsToken_, triggerResults);

  bool EventPasses = true;
  // Make sure handle is valid
  if (TheHLTResults.isValid() && hltselection_) {
    //Get HLT Names
    const edm::TriggerNames& TheTriggerNames = iEvent.triggerNames(*TheHLTResults);

    for (unsigned int index = 0; index < HLTBitLabel_.size(); index++) {
      if (!HLTBitLabel_[index].label().empty()) {
        //Change the default value since HLT requirement has been issued by the user
        if (index == 0)
          EventPasses = false;
        //Get the HLT bit and check to make sure it is valid
        unsigned int bit = TheTriggerNames.triggerIndex(HLTBitLabel_[index].label());
        if (bit < TheHLTResults->size()) {
          //If any of the HLT names given by the user accept, then the event passes
          if (TheHLTResults->accept(bit) && !TheHLTResults->error(bit)) {
            EventPasses = true;
            hCT_NEvents_HLT[index]->Fill(1);
          } else
            hCT_NEvents_HLT[index]->Fill(0);
        } else {
          edm::LogInfo("OutputInfo") << "The HLT Trigger Name : " << HLTBitLabel_[index].label()
                                     << " is not valid for this trigger table " << std::endl;
        }
      }
    }
  }

  if (!EventPasses && hltselection_)
    return;

  //----------GREG & CHRIS' idea---///
  float ETTowerMin = -1;  //GeV
  float METRingMin = -2;  // GeV

  Nevents++;
  hCT_Nevents->Fill(0);

  // ==========================================================
  // Retrieve!
  // ==========================================================

  edm::Handle<edm::View<Candidate> > towers;
  iEvent.getByToken(caloTowersLabel_, towers);

  if ((!towers.isValid())) {
    //DD:fix print label
    //edm::LogInfo("")<<"CaloTowers "<< caloTowersLabel_<<" not found!"<<std::endl;
    return;
  }

  //HBHENoiseFilterResultToken_=consumes<HBHENoiseFilter>(edm::InputTag(HBHENoiseFilterResultLabel_));
  edm::Handle<bool> HBHENoiseFilterResultHandle;
  iEvent.getByToken(HBHENoiseFilterResultLabel_, HBHENoiseFilterResultHandle);
  bool HBHENoiseFilterResult = *HBHENoiseFilterResultHandle;
  if (!HBHENoiseFilterResultHandle.isValid()) {
    LogDebug("") << "CaloTowerAnalyzer: Could not find HBHENoiseFilterResult" << std::endl;
  }

  bool bHcalNoiseFilter = HBHENoiseFilterResult;

  if (!bHcalNoiseFilter)
    return;

  edm::View<Candidate>::const_iterator towerCand = towers->begin();

  // ==========================================================
  // Fill Histograms!
  // ==========================================================

  int CTmin_iphi = 99, CTmax_iphi = -99;
  int CTmin_ieta = 99, CTmax_ieta = -99;

  TLorentzVector vMET_EtaRing[83];
  int ActiveRing[83];
  int NActiveTowers[83];
  double SET_EtaRing[83];
  double MinEt_EtaRing[83];
  double MaxEt_EtaRing[83];
  for (int i = 0; i < 83; i++) {
    ActiveRing[i] = 0;
    NActiveTowers[i] = 0;
    SET_EtaRing[i] = 0;
    MinEt_EtaRing[i] = 0;
    MaxEt_EtaRing[i] = 0;
  }

  //  for (CaloTowerCollection::const_iterator calotower = towers->begin(); calotower != towers->end(); calotower++)
  for (; towerCand != towers->end(); towerCand++) {
    const Candidate* candidate = &(*towerCand);
    if (candidate) {
      const CaloTower* calotower = dynamic_cast<const CaloTower*>(candidate);
      if (calotower) {
        //math::RhoEtaPhiVector Momentum = calotower->momentum();
        double Tower_ET = calotower->et();
        //double Tower_Energy  = calotower->energy();
        //    double Tower_Eta = calotower->eta();
        double Tower_Phi = calotower->phi();
        //double Tower_EMEnergy = calotower->emEnergy();
        //double Tower_HadEnergy = calotower->hadEnergy();
        double Tower_OuterEt = calotower->outerEt();
        double Tower_EMEt = calotower->emEt();
        double Tower_HadEt = calotower->hadEt();
        //int Tower_EMLV1 = calotower->emLvl1();
        //int Tower_HadLV1 = calotower->hadLv11();
        int Tower_ieta = calotower->id().ieta();
        int Tower_iphi = calotower->id().iphi();
        int EtaRing = 41 + Tower_ieta;
        ActiveRing[EtaRing] = 1;
        NActiveTowers[EtaRing]++;
        SET_EtaRing[EtaRing] += Tower_ET;
        TLorentzVector v_;
        v_.SetPtEtaPhiE(Tower_ET, 0, Tower_Phi, Tower_ET);
        if (Tower_ET > ETTowerMin)
          vMET_EtaRing[EtaRing] -= v_;

        // Fill Histograms
        hCT_Occ_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
        if (calotower->emEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
          hCT_Occ_EM_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
        if (calotower->hadEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
          hCT_Occ_HAD_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);
        if (calotower->outerEt() > 0 && calotower->emEt() + calotower->hadEt() > 0.3)
          hCT_Occ_Outer_Et_ieta_iphi->Fill(Tower_ieta, Tower_iphi);

        hCT_et_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_ET);
        hCT_emEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_EMEt);
        hCT_hadEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_HadEt);
        hCT_outerEt_ieta_iphi->Fill(Tower_ieta, Tower_iphi, Tower_OuterEt);

        if (allhist_) {
          hCT_etvsieta->Fill(Tower_ieta, Tower_ET);
          hCT_emEtvsieta->Fill(Tower_ieta, Tower_EMEt);
          hCT_hadEtvsieta->Fill(Tower_ieta, Tower_HadEt);
          hCT_outerEtvsieta->Fill(Tower_ieta, Tower_OuterEt);
        }

        if (Tower_ET > MaxEt_EtaRing[EtaRing])
          MaxEt_EtaRing[EtaRing] = Tower_ET;
        if (Tower_ET < MinEt_EtaRing[EtaRing] && Tower_ET > 0)
          MinEt_EtaRing[EtaRing] = Tower_ET;

        if (Tower_ieta < CTmin_ieta)
          CTmin_ieta = Tower_ieta;
        if (Tower_ieta > CTmax_ieta)
          CTmax_ieta = Tower_ieta;
        if (Tower_iphi < CTmin_iphi)
          CTmin_iphi = Tower_iphi;
        if (Tower_iphi > CTmax_iphi)
          CTmax_iphi = Tower_iphi;
      }  //end if (calotower) ..
    }    // end if(candidate) ...

  }  // end loop over towers

  // Fill eta-ring MET quantities
  if (allhist_) {
    for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
      hCT_Minetvsieta->Fill(iEtaRing - 41, MinEt_EtaRing[iEtaRing]);
      hCT_Maxetvsieta->Fill(iEtaRing - 41, MaxEt_EtaRing[iEtaRing]);

      if (ActiveRing[iEtaRing]) {
        if (vMET_EtaRing[iEtaRing].Pt() > METRingMin) {
          hCT_METPhivsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Phi());
          hCT_MExvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Px());
          hCT_MEyvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Py());
          hCT_METvsieta->Fill(iEtaRing - 41, vMET_EtaRing[iEtaRing].Pt());
        }
        hCT_SETvsieta->Fill(iEtaRing - 41, SET_EtaRing[iEtaRing]);
        hCT_Occvsieta->Fill(iEtaRing - 41, NActiveTowers[iEtaRing]);
      }
    }  // ietaring
  }    // allhist_

  edm::LogInfo("OutputInfo") << "CT ieta range: " << CTmin_ieta << " " << CTmax_ieta;
  edm::LogInfo("OutputInfo") << "CT iphi range: " << CTmin_iphi << " " << CTmax_iphi;
}