ECALMultifitAnalyzer_HI.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    DQMOffline/ECALMultifitAnalyzer_HI
// Class:      ECALMultifitAnalyzer_HI
//
//
// Original Author:  R. Alex Barbieri
//         Created:  Tue, 17 Nov 2015 17:18:50 GMT
// Modified for DQM: Raghav Kunnawalkam Elayavalli
//                   Wednesday 18 Nov 2015 9:45:AM GVA time
//
//

// system include files
#include <memory>
#include <vector>

// user include files
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"

#include "Geometry/Records/interface/CaloTopologyRecord.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "TH2F.h"

//
// class declaration
//

class ECALMultifitAnalyzer_HI : public DQMEDAnalyzer {
public:
  explicit ECALMultifitAnalyzer_HI(const edm::ParameterSet &);
  ~ECALMultifitAnalyzer_HI() override {}

private:
  void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
  void analyze(const edm::Event &, const edm::EventSetup &) override;

  edm::EDGetTokenT<std::vector<reco::Photon>> recoPhotonsCollection_;
  edm::EDGetTokenT<reco::CaloJetCollection> caloJetToken_;
  edm::EDGetTokenT<EcalRecHitCollection> RecHitCollection_EB_;
  edm::EDGetTokenT<EcalRecHitCollection> RecHitCollection_EE_;
  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomH;

  double mRechitEnergyThreshold;
  double mRecoPhotonPtThreshold;
  double mRecoJetPtThreshold;
  double mDeltaRPhotonThreshold;
  double mDeltaRJetThreshold;

  MonitorElement *eb_chi2;
  MonitorElement *eb_chi2_eta;
  MonitorElement *eb_chi2_e5;
  MonitorElement *eb_chi2_e5_eta;
  MonitorElement *eb_errors;
  MonitorElement *eb_errors_eta;
  MonitorElement *eb_errors_e5;
  MonitorElement *eb_errors_e5_eta;
  MonitorElement *eb_chi2_photon15;
  MonitorElement *eb_errors_photon15;
  MonitorElement *eb_chi2_jet30;
  MonitorElement *eb_errors_jet30;

  MonitorElement *ee_chi2;
  MonitorElement *ee_chi2_eta;
  MonitorElement *ee_chi2_e5;
  MonitorElement *ee_chi2_e5_eta;
  MonitorElement *ee_errors;
  MonitorElement *ee_errors_eta;
  MonitorElement *ee_errors_e5;
  MonitorElement *ee_errors_e5_eta;
  MonitorElement *ee_chi2_photon15;
  MonitorElement *ee_errors_photon15;
  MonitorElement *ee_chi2_jet30;
  MonitorElement *ee_errors_jet30;
};

//
// constructors and destructor
//
ECALMultifitAnalyzer_HI::ECALMultifitAnalyzer_HI(const edm::ParameterSet &iConfig)
    : recoPhotonsCollection_(consumes<std::vector<reco::Photon>>(iConfig.getParameter<edm::InputTag>("recoPhotonSrc"))),
      caloJetToken_(consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("recoJetSrc"))),
      RecHitCollection_EB_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("RecHitCollection_EB"))),
      RecHitCollection_EE_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("RecHitCollection_EE"))),
      geomH(esConsumes()),
      mRechitEnergyThreshold(iConfig.getParameter<double>("rechitEnergyThreshold")),
      mRecoPhotonPtThreshold(iConfig.getParameter<double>("recoPhotonPtThreshold")),
      mRecoJetPtThreshold(iConfig.getParameter<double>("recoJetPtThreshold")),
      mDeltaRPhotonThreshold(iConfig.getParameter<double>("deltaRPhotonThreshold")),
      mDeltaRJetThreshold(iConfig.getParameter<double>("deltaRJetThreshold")) {}

//
// member functions
//

// ------------ method called for each event  ------------
void ECALMultifitAnalyzer_HI::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
  using namespace edm;

  const CaloGeometry *geom = &iSetup.getData(geomH);

  Handle<std::vector<reco::Photon>> recoPhotonsHandle;
  iEvent.getByToken(recoPhotonsCollection_, recoPhotonsHandle);

  Handle<reco::CaloJetCollection> recoJetHandle;
  iEvent.getByToken(caloJetToken_, recoJetHandle);

  Handle<EcalRecHitCollection> ebHandle;
  iEvent.getByToken(RecHitCollection_EB_, ebHandle);

  Handle<EcalRecHitCollection> eeHandle;
  iEvent.getByToken(RecHitCollection_EE_, eeHandle);

  for (EcalRecHitCollection::const_iterator hit = ebHandle->begin(); hit != ebHandle->end(); ++hit) {
    eb_chi2->Fill(hit->chi2());
    eb_errors->Fill(hit->energyError());
    double eta = geom->getGeometry(hit->detid())->getPosition().eta();
    double phi = geom->getGeometry(hit->detid())->getPosition().phi();
    eb_chi2_eta->Fill(eta, hit->chi2());
    eb_errors_eta->Fill(eta, hit->energyError());
    if (hit->energy() > mRechitEnergyThreshold) {
      eb_chi2_e5->Fill(hit->chi2());
      eb_errors_e5->Fill(hit->energyError());
      eb_chi2_e5_eta->Fill(eta, hit->chi2());
      eb_errors_e5_eta->Fill(eta, hit->energyError());
    }

    for (std::vector<reco::Photon>::const_iterator pho = recoPhotonsHandle->begin(); pho != recoPhotonsHandle->end();
         ++pho) {
      if (pho->et() < mRecoPhotonPtThreshold)
        continue;
      double dr = reco::deltaR(eta, phi, pho->eta(), pho->phi());
      if (dr < mDeltaRPhotonThreshold) {
        eb_chi2_photon15->Fill(hit->chi2());
        eb_errors_photon15->Fill(hit->energyError());
      }
    }
    for (std::vector<reco::CaloJet>::const_iterator jet = recoJetHandle->begin(); jet != recoJetHandle->end(); ++jet) {
      if (jet->pt() < mRecoJetPtThreshold)
        continue;
      double dr = reco::deltaR(eta, phi, jet->eta(), jet->phi());
      if (dr < mDeltaRJetThreshold) {
        eb_chi2_jet30->Fill(hit->chi2());
        eb_errors_jet30->Fill(hit->energyError());
      }
    }
  }

  for (EcalRecHitCollection::const_iterator hit = eeHandle->begin(); hit != eeHandle->end(); ++hit) {
    ee_chi2->Fill(hit->chi2());
    ee_errors->Fill(hit->energyError());
    double eta = geom->getGeometry(hit->detid())->getPosition().eta();
    double phi = geom->getGeometry(hit->detid())->getPosition().phi();
    ee_chi2_eta->Fill(eta, hit->chi2());
    ee_errors_eta->Fill(eta, hit->energyError());
    if (hit->energy() > mRechitEnergyThreshold) {
      ee_chi2_e5->Fill(hit->chi2());
      ee_errors_e5->Fill(hit->energyError());
      ee_chi2_e5_eta->Fill(eta, hit->chi2());
      ee_errors_e5_eta->Fill(eta, hit->energyError());
    }

    for (std::vector<reco::Photon>::const_iterator pho = recoPhotonsHandle->begin(); pho != recoPhotonsHandle->end();
         ++pho) {
      if (pho->et() < mRecoPhotonPtThreshold)
        continue;
      double dr = reco::deltaR(eta, phi, pho->eta(), pho->phi());
      if (dr < mDeltaRPhotonThreshold) {
        ee_chi2_photon15->Fill(hit->chi2());
        ee_errors_photon15->Fill(hit->energyError());
      }
    }
    for (std::vector<reco::CaloJet>::const_iterator jet = recoJetHandle->begin(); jet != recoJetHandle->end(); ++jet) {
      if (jet->pt() < mRecoJetPtThreshold)
        continue;
      double dr = reco::deltaR(eta, phi, jet->eta(), jet->phi());
      if (dr < mDeltaRJetThreshold) {
        ee_chi2_jet30->Fill(hit->chi2());
        ee_errors_jet30->Fill(hit->energyError());
      }
    }
  }
}

// ------------ method called once each job just before starting event loop
// ------------
void ECALMultifitAnalyzer_HI::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
  iBooker.setCurrentFolder("EcalCalibration/HIN_MultiFitAnalyzer");

  eb_chi2 = nullptr;
  eb_chi2_eta = nullptr;
  eb_chi2_e5 = nullptr;
  eb_chi2_e5_eta = nullptr;
  eb_errors = nullptr;
  eb_errors_eta = nullptr;
  eb_errors_e5 = nullptr;
  eb_errors_e5_eta = nullptr;
  eb_chi2_photon15 = nullptr;
  eb_errors_photon15 = nullptr;
  eb_chi2_jet30 = nullptr;
  eb_errors_jet30 = nullptr;

  ee_chi2 = nullptr;
  ee_chi2_eta = nullptr;
  ee_chi2_e5 = nullptr;
  ee_chi2_e5_eta = nullptr;
  ee_errors = nullptr;
  ee_errors_eta = nullptr;
  ee_errors_e5 = nullptr;
  ee_errors_e5_eta = nullptr;
  ee_chi2_photon15 = nullptr;
  ee_errors_photon15 = nullptr;
  ee_chi2_jet30 = nullptr;
  ee_errors_jet30 = nullptr;

  const int nBins = 500;
  const float maxChi2 = 70;
  const float maxError = 0.5;

  TH2F *hProfile_Chi2 = new TH2F("hProfile_Chi2", "", nBins, -5, 5, nBins, 0, maxChi2);
  TH2F *hProfile_Err = new TH2F("hProfile_Err", "", nBins, -5, 5, nBins, 0, maxError);

  eb_chi2 = iBooker.book1D("rechit_eb_chi2", "Rechit eb_chi2;chi2 fit value;", nBins, 0, maxChi2);
  eb_chi2_eta = iBooker.book2D("rechit_eb_eta_Vs_mean_Chi2", hProfile_Chi2);
  eb_chi2_e5 = iBooker.book1D(Form("rechit_eb_chi2_e%2.0f", mRechitEnergyThreshold),
                              Form("Rechit eb_chi2, e>%2.0fGeV;chi2 fit value;", mRechitEnergyThreshold),
                              nBins,
                              0,
                              maxChi2);
  eb_chi2_e5_eta = iBooker.book2D(Form("rechit_eb_chi2_e%2.0f_eta", mRechitEnergyThreshold), hProfile_Chi2);
  eb_errors = iBooker.book1D("rechit_eb_errors", "Rechit eb_errors;error on the energy;", nBins, 0, maxError);
  eb_errors_eta = iBooker.book2D("rechit_eb_errors_eta", hProfile_Err);
  eb_errors_e5 = iBooker.book1D(Form("rechit_eb_errors_e%2.0f", mRechitEnergyThreshold),
                                "Rechit eb_errors, e>5GeV;error on the energy;",
                                nBins,
                                0,
                                maxError);
  eb_errors_e5_eta = iBooker.book2D(Form("rechit_eb_errors_e%2.0f_eta", mRechitEnergyThreshold), hProfile_Err);
  eb_chi2_photon15 = iBooker.book1D(Form("rechit_eb_chi2_photon%2.0f", mRecoPhotonPtThreshold),
                                    "Rechit eb_chi2 near photons;chi2 fit value;",
                                    nBins,
                                    0,
                                    maxChi2);
  eb_errors_photon15 = iBooker.book1D(Form("rechit_eb_errors_photon%2.0f", mRecoPhotonPtThreshold),
                                      "Rechit eb_errors near photons;error on the energy;",
                                      nBins,
                                      0,
                                      maxError);
  eb_chi2_jet30 = iBooker.book1D(Form("rechit_eb_chi2_jet%2.0f", mRecoJetPtThreshold),
                                 "Rechit eb_chi2 near jets;chi2 fit value;",
                                 nBins,
                                 0,
                                 maxChi2);
  eb_errors_jet30 = iBooker.book1D(Form("rechit_eb_errors_jet%2.0f", mRecoJetPtThreshold),
                                   "Rechit eb_errors near jets;error on the energy;",
                                   nBins,
                                   0,
                                   maxError);

  ee_chi2 = iBooker.book1D("rechit_ee_chi2", "Rechit ee_chi2;chi2 fit value;", nBins, 0, maxChi2);
  ee_chi2_eta = iBooker.book2D("rechit_ee_chi2_eta", hProfile_Chi2);
  ee_chi2_e5 = iBooker.book1D(Form("rechit_ee_chi2_e%2.0f", mRechitEnergyThreshold),
                              "Rechit ee_chi2, e>5GeV;chi2 fit value;",
                              nBins,
                              0,
                              maxChi2);
  ee_chi2_e5_eta = iBooker.book2D(Form("rechit_ee_chi2_e%2.0f_eta", mRechitEnergyThreshold), hProfile_Chi2);
  ee_errors = iBooker.book1D("rechit_ee_errors", "Rechit ee_errors;error on the energy;", nBins, 0, maxError);
  ee_errors_eta = iBooker.book2D("rechit_ee_errors_eta", hProfile_Err);
  ee_errors_e5 = iBooker.book1D(Form("rechit_ee_errors_e%2.0f", mRechitEnergyThreshold),
                                "Rechit ee_errors, e>5GeV;error on the energy;",
                                nBins,
                                0,
                                maxError);
  ee_errors_e5_eta = iBooker.book2D(Form("rechit_ee_errors_e%2.0f_eta", mRechitEnergyThreshold), hProfile_Err);
  ee_chi2_photon15 = iBooker.book1D(Form("rechit_ee_chi2_photon%2.0f", mRecoPhotonPtThreshold),
                                    "Rechit ee_chi2 near photons;chi2 fit value;",
                                    nBins,
                                    0,
                                    maxChi2);
  ee_errors_photon15 = iBooker.book1D(Form("rechit_ee_errors_photon%2.0f", mRecoPhotonPtThreshold),
                                      "Rechit ee_errors near photons;error on the energy;",
                                      nBins,
                                      0,
                                      maxError);
  ee_chi2_jet30 = iBooker.book1D(Form("rechit_ee_chi2_jet%2.0f", mRecoJetPtThreshold),
                                 "Rechit ee_chi2 near jets;chi2 fit value;",
                                 nBins,
                                 0,
                                 maxChi2);
  ee_errors_jet30 = iBooker.book1D(Form("rechit_ee_errors_jet%2.0f", mRecoJetPtThreshold),
                                   "Rechit ee_errors near jets;error on the energy;",
                                   nBins,
                                   0,
                                   maxError);

  delete hProfile_Chi2;
  delete hProfile_Err;
}

// define this as a plug-in
DEFINE_FWK_MODULE(ECALMultifitAnalyzer_HI);