EcalPileUpDepMonitor.cc

Go to the documentation of this file.

/*
 * \file EcalPileUpDepMonitor.cc
 * \author Ben Carlson - CMU
 * Last Update:
 *
 */

#include "DQMOffline/Ecal/interface/EcalPileUpDepMonitor.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/VertexReco/interface/Vertex.h"

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

#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"

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

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

#include <cmath>

// Framework

EcalPileUpDepMonitor::EcalPileUpDepMonitor(const edm::ParameterSet &ps) {
  VertexCollection_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("VertexCollection"));

  if (ps.existsAs<edm::InputTag>("basicClusterCollection") &&
      !ps.getParameter<edm::InputTag>("basicClusterCollection").label().empty())
    basicClusterCollection_ =
        consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection"));
  else {
    basicClusterCollection_EE_ =
        consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection_EE"));
    basicClusterCollection_EB_ =
        consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection_EB"));
  }
  superClusterCollection_EB_ =
      consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("superClusterCollection_EB"));
  superClusterCollection_EE_ =
      consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("superClusterCollection_EE"));

  RecHitCollection_EB_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("RecHitCollection_EB"));
  RecHitCollection_EE_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("RecHitCollection_EE"));
  EleTag_ = consumes<reco::GsfElectronCollection>(ps.getParameter<edm::InputTag>("EleTag"));
}

EcalPileUpDepMonitor::~EcalPileUpDepMonitor() {}

void EcalPileUpDepMonitor::bookHistograms(DQMStore::IBooker &ibooker,
                                          edm::Run const &,
                                          edm::EventSetup const &eventSetup) {
  eventSetup.get<CaloGeometryRecord>().get(geomH);
  eventSetup.get<CaloTopologyRecord>().get(caloTop);

  ibooker.cd();
  ibooker.setCurrentFolder("Ecal/EcalPileUpDepMonitor");

  std::string prof_name = "bcEB_PV";
  std::string title = "Basic clusters EB vs. PV";
  bcEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  bcEB_PV->setAxisTitle("N_{pv}", 1);
  bcEB_PV->setAxisTitle("Basic Clusters", 2);

  prof_name = "bcEE_PV";
  title = "Basic Clusters EE vs. PV";
  bcEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  bcEE_PV->setAxisTitle("N_{pv}", 1);
  bcEE_PV->setAxisTitle("Basic Clusters", 2);

  prof_name = "scEB_PV";
  title = "Super Clusters EB vs. PV";
  scEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  scEB_PV->setAxisTitle("N_{pv}", 1);
  scEB_PV->setAxisTitle("Super Clusters", 2);

  prof_name = "scEE_PV";
  title = "Super Clusters EE vs. PV";
  scEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  scEE_PV->setAxisTitle("N_{pv}", 1);
  scEE_PV->setAxisTitle("Super Clusters", 2);

  prof_name = "scEtEB_PV";
  title = "Super Clusters Et EB vs. PV";
  scEtEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  scEtEB_PV->setAxisTitle("N_{pv}", 1);
  scEtEB_PV->setAxisTitle("Super Cluster E_{T} [GeV]", 2);

  prof_name = "scEtEE_PV";
  title = "Super Clusters Et EE vs. PV";
  scEtEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  scEtEE_PV->setAxisTitle("N_{pv}", 1);
  scEtEE_PV->setAxisTitle("Super Cluster E_{T} [GeV]", 2);

  prof_name = "recHitEtEB_PV";
  title = "Reconstructed Hit Et EB vs. PV";
  recHitEtEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  recHitEtEB_PV->setAxisTitle("N_{pv}", 1);
  recHitEtEB_PV->setAxisTitle("Reconstructed hit E_{T} [GeV]", 2);

  prof_name = "recHitEtEE_PV";
  title = "Reconstructed Hit Et EE vs. PV";
  recHitEtEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
  recHitEtEE_PV->setAxisTitle("N_{pv}", 1);
  recHitEtEE_PV->setAxisTitle("Reconstructed hit E_{T} [GeV]", 2);

  prof_name = "emIso_PV";
  title = "EM Isolation vs. PV";
  emIso_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350);
  emIso_PV->setAxisTitle("N_{pv}", 1);
  emIso_PV->setAxisTitle("EM_{Isolation} [GeV]", 2);

  prof_name = "emIso";
  title = "EM Isolation";
  emIso = ibooker.book1D(prof_name, title, 50, 0, 50);
  emIso->setAxisTitle("EM_{Isolation} [GeV]", 1);
  emIso->setAxisTitle("Events", 2);

  prof_name = "scHitEtEB";
  title = "Super Cluster Hit Et EB";
  scHitEtEB = ibooker.book1D(prof_name, title, 50, 0, 100);
  scHitEtEB->setAxisTitle("super cluster hit E_{T} [GeV]", 1);
  scHitEtEB->setAxisTitle("Events", 2);

  prof_name = "scHitEtEE";
  title = "Super Cluster Hit Et EE";
  scHitEtEE = ibooker.book1D(prof_name, title, 50, 0, 100);
  scHitEtEE->setAxisTitle("super cluster hit E_{T} [GeV]", 1);
  scHitEtEE->setAxisTitle("Events", 2);

  //   prof_name="scHitE_EB";
  //   title="Super Cluster Hit E EB";
  //   scHitE_EB=ibooker.book1D(prof_name,title,50,0,100);
  //   scHitE_EB->setAxisTitle("super cluster hit E [GeV]",1);
  //   scHitE_EB->setAxisTitle("Events",2);

  //   prof_name="scHitE_EE";
  //   title="Super Cluster Hit E EE";
  //   scHitE_EE=ibooker.book1D(prof_name,title,50,0,100);
  //   scHitE_EE->setAxisTitle("super cluster hit E [GeV]",1);
  //   scHitE_EE->setAxisTitle("Events",2);

  // SC eta
  //   prof_name="scEta_EB";
  //   title="Super Cluster #eta EB";
  //   scEta_EB=ibooker.book1D(prof_name,title,50,-6,6);
  //   scEta_EB->setAxisTitle("#eta",1);
  //   scEta_EB->setAxisTitle("Events",2);

  //   prof_name="scEta_EE";
  //   title="Super Cluster #eta EE";
  //   scEta_EE=ibooker.book1D(prof_name,title,50,-6,6);
  //   scEta_EE->setAxisTitle("#eta",1);
  //   scEta_EE->setAxisTitle("Events",2);

  // SC phi
  //   prof_name="scPhi_EB";
  //   title="Super Cluster #phi EB";
  //   scPhi_EB=ibooker.book1D(prof_name,title,50,-3.14,3.14);
  //   scPhi_EB->setAxisTitle("super cluster #phi",1);
  //   scPhi_EB->setAxisTitle("Events",2);

  //   prof_name="scPhi_EE";
  //   title="Super Cluster #phi EE";
  //   scPhi_EE=ibooker.book1D(prof_name,title,50,-3.14,3.14);
  //   scPhi_EE->setAxisTitle("super cluster #phi",1);
  //   scPhi_EE->setAxisTitle("Events",2);

  // sc sigma eta eta / eta phi

  prof_name = "scSigmaIetaIeta_EB";
  title = "Super Cluster sigmaIetaIeta EB";
  scSigmaIetaIeta_EB = ibooker.book1D(prof_name, title, 50, 0, 0.03);
  scSigmaIetaIeta_EB->setAxisTitle("#sigma_{i#etai#eta}", 1);
  scSigmaIetaIeta_EB->setAxisTitle("Events", 2);

  prof_name = "scSigmaIetaIeta_EE";
  title = "Super Cluster sigmaIetaIeta EE";
  scSigmaIetaIeta_EE = ibooker.book1D(prof_name, title, 50, 0, 0.1);
  scSigmaIetaIeta_EE->setAxisTitle("#sigma_{i#etai#eta}", 1);
  scSigmaIetaIeta_EE->setAxisTitle("Events", 2);

  // phi
  prof_name = "scSigmaIetaIphi_EB";
  title = "Super Cluster sigmaIetaIphi EB";
  scSigmaIetaIphi_EB = ibooker.book1D(prof_name, title, 50, -5.e-4, 5.e-4);
  scSigmaIetaIphi_EB->setAxisTitle("#sigma_{i#etai#phi}", 1);
  scSigmaIetaIphi_EB->setAxisTitle("Events", 2);

  prof_name = "scSigmaIetaIphi_EE";
  title = "Super Cluster sigmaIetaIphi EE";
  scSigmaIetaIphi_EE = ibooker.book1D(prof_name, title, 50, -2.5e-3, 2.5e-3);
  scSigmaIetaIphi_EE->setAxisTitle("#sigma_{i#etai#phi}", 1);
  scSigmaIetaIphi_EE->setAxisTitle("Events", 2);

  // R9
  prof_name = "r9_EB";
  title = "r9 EB";
  r9_EB = ibooker.book1D(prof_name, title, 50, 0, 1.5);
  r9_EB->setAxisTitle("R_{9}", 1);
  r9_EB->setAxisTitle("Events", 2);

  prof_name = "r9_EE";
  title = "r9 EE";
  r9_EE = ibooker.book1D(prof_name, title, 50, 0, 1.5);
  r9_EE->setAxisTitle("R_{9}", 1);
  r9_EE->setAxisTitle("Events", 2);

  // Rec Hit

  prof_name = "recHitEtEB";
  title = "RecHit Et EB";
  recHitEtEB = ibooker.book1D(prof_name, title, 50, 0, 400);
  recHitEtEB->setAxisTitle("Reconstructed Hit E_{T} [GeV]", 1);
  recHitEtEB->setAxisTitle("Events", 2);

  prof_name = "recHitEtEE";
  title = "RecHit Et EE";
  recHitEtEE = ibooker.book1D(prof_name, title, 50, 0, 400);
  recHitEtEE->setAxisTitle("Reconstructed Hit E_{T} [GeV]", 1);
  recHitEtEE->setAxisTitle("Events", 2);
}

void EcalPileUpDepMonitor::analyze(const edm::Event &e, const edm::EventSetup &) {
  const CaloGeometry *geom = geomH.product();
  // Vertex collection:
  //-----------------------------------------
  edm::Handle<reco::VertexCollection> PVCollection_h;
  e.getByToken(VertexCollection_, PVCollection_h);
  if (!PVCollection_h.isValid()) {
    edm::LogWarning("VertexCollection") << "VertexCollection not found";
  }
  //-----------------gsfElectrons -------------------------
  edm::Handle<reco::GsfElectronCollection> electronCollection_h;
  e.getByToken(EleTag_, electronCollection_h);
  if (!electronCollection_h.isValid()) {
    edm::LogWarning("EBRecoSummary") << "Electrons not found";
  }

  if (basicClusterCollection_.isUninitialized()) {
    //----------------- Basic Cluster Collection Ecal Barrel  ---------
    edm::Handle<edm::View<reco::CaloCluster>> basicClusters_EB_h;
    e.getByToken(basicClusterCollection_EB_, basicClusters_EB_h);
    if (!basicClusters_EB_h.isValid()) {
      edm::LogWarning("EBRecoSummary") << "basicClusters_EB_h not found";
    }

    bcEB_PV->Fill(PVCollection_h->size(), basicClusters_EB_h->size());

    //----------------- Basic Cluster Collection Ecal Endcal  ---------

    edm::Handle<edm::View<reco::CaloCluster>> basicClusters_EE_h;
    e.getByToken(basicClusterCollection_EE_, basicClusters_EE_h);
    if (!basicClusters_EE_h.isValid()) {
      edm::LogWarning("EERecoSummary") << "basicClusters_EE_h not found";
    }

    bcEE_PV->Fill(PVCollection_h->size(), basicClusters_EE_h->size());
  } else {
    //----------------- Basic Cluster Collection Ecal Barrel  ---------
    edm::Handle<edm::View<reco::CaloCluster>> basicClusters_h;
    e.getByToken(basicClusterCollection_, basicClusters_h);
    if (!basicClusters_h.isValid()) {
      edm::LogWarning("EBRecoSummary") << "basicClusters_h not found";
    }

    int nBarrel(0);
    int nEndcap(0);
    for (edm::View<reco::CaloCluster>::const_iterator bcItr(basicClusters_h->begin()); bcItr != basicClusters_h->end();
         ++bcItr) {
      if (bcItr->caloID().detector(reco::CaloID::DET_ECAL_BARREL))
        ++nBarrel;
      if (bcItr->caloID().detector(reco::CaloID::DET_ECAL_ENDCAP))
        ++nEndcap;
    }

    bcEB_PV->Fill(PVCollection_h->size(), nBarrel);
    bcEE_PV->Fill(PVCollection_h->size(), nEndcap);
  }

  //----------------- Reconstructed Hit Ecal barrel

  edm::Handle<EcalRecHitCollection> RecHitsEB;
  e.getByToken(RecHitCollection_EB_, RecHitsEB);
  if (!RecHitsEB.isValid()) {
    edm::LogWarning("EBRecoSummary") << "RecHitsEB not found";
  }

  //----------------- Reconstructed Hit Ecal Endcap

  edm::Handle<EcalRecHitCollection> RecHitsEE;
  e.getByToken(RecHitCollection_EE_, RecHitsEE);
  if (!RecHitsEE.isValid()) {
    edm::LogWarning("EBRecoSummary") << "RecHitsEB not found";
  }
  //----------------- Super Cluster Collection Ecal Endcap  ---------

  edm::Handle<reco::SuperClusterCollection> superClusters_EE_h;
  e.getByToken(superClusterCollection_EE_, superClusters_EE_h);
  if (!superClusters_EE_h.isValid()) {
    edm::LogWarning("EERecoSummary") << "superClusters_EE_h not found";
  }

  //--------- Fill Isolation -----------------

  if (electronCollection_h.isValid()) {
    for (reco::GsfElectronCollection::const_iterator recoElectron = electronCollection_h->begin();
         recoElectron != electronCollection_h->end();
         recoElectron++) {
      double IsoEcal = recoElectron->dr03EcalRecHitSumEt();  
      emIso_PV->Fill(PVCollection_h->size(), IsoEcal);
      emIso->Fill(IsoEcal);
    }
  }

  // fill super clusters EE
  scEE_PV->Fill(PVCollection_h->size(), superClusters_EE_h->size());

  for (reco::SuperClusterCollection::const_iterator itSC = superClusters_EE_h->begin();
       itSC != superClusters_EE_h->end();
       ++itSC) {
    double scEE_Et = itSC->energy() * sin(2. * atan(exp(-itSC->position().eta())));
    //    double scEE_E=itSC->energy();

    // fill super cluster endcap eta/phi
    //     scEta_EE->Fill(itSC->position().eta());
    //     scPhi_EE->Fill(itSC->position().phi());

    // get sigma eta_eta etc

    CaloTopology const *p_topology = caloTop.product();  // get calo topology
    const EcalRecHitCollection *eeRecHits = RecHitsEE.product();

    reco::BasicCluster const &seedCluster(*itSC->seed());
    std::vector<float> cov = EcalClusterTools::localCovariances(seedCluster, eeRecHits, p_topology);
    float sigmaIetaIeta = std::sqrt(cov[0]);
    float sigmaIetaIphi = cov[1];

    float e3x3 = EcalClusterTools::e3x3(seedCluster, eeRecHits, p_topology);
    float r9 = e3x3 / itSC->rawEnergy();

    r9_EE->Fill(r9);
    scSigmaIetaIeta_EE->Fill(sigmaIetaIeta);
    scSigmaIetaIphi_EE->Fill(sigmaIetaIphi);

    // std::cout  << " sigmaIetaIeta: " << sigmaIetaIeta << std::endl;
    scEtEE_PV->Fill(PVCollection_h->size(), scEE_Et);
    scHitEtEE->Fill(scEE_Et);  // super cluster Et historam
    //    scHitE_EE->Fill(scEE_E); //super cluster energy histogram

  }  // sc-EE loop

  //----------------- Super Cluster Collection Ecal Barrel  ---------

  edm::Handle<reco::SuperClusterCollection> superClusters_EB_h;
  e.getByToken(superClusterCollection_EB_, superClusters_EB_h);
  if (!superClusters_EB_h.isValid()) {
    edm::LogWarning("EBRecoSummary") << "superClusters_EB_h not found";
  }
  scEB_PV->Fill(PVCollection_h->size(), superClusters_EB_h->size());

  for (reco::SuperClusterCollection::const_iterator itSC = superClusters_EB_h->begin();
       itSC != superClusters_EB_h->end();
       ++itSC) {
    double scEB_Et = itSC->energy() * sin(2. * atan(exp(-itSC->position().eta())));  // super cluster transverse energy
    //    double scEB_E= itSC->energy(); // super cluster energy

    // fill super cluster Barrel eta/phi
    //     scEta_EB->Fill(itSC->position().eta()); //super cluster eta
    //     scPhi_EB->Fill(itSC->position().phi()); // super cluster phi

    // sigma ietaieta etc

    CaloTopology const *p_topology = caloTop.product();  // get calo topology
    const EcalRecHitCollection *ebRecHits = RecHitsEB.product();

    reco::BasicCluster const &seedCluster(*itSC->seed());
    std::vector<float> cov = EcalClusterTools::localCovariances(seedCluster, ebRecHits, p_topology);
    float sigmaIetaIeta = std::sqrt(cov[0]);
    float sigmaIetaIphi = cov[1];

    float e3x3 = EcalClusterTools::e3x3(seedCluster, ebRecHits, p_topology);
    float r9 = e3x3 / itSC->rawEnergy();

    r9_EB->Fill(r9);
    scSigmaIetaIeta_EB->Fill(sigmaIetaIeta);
    scSigmaIetaIphi_EB->Fill(sigmaIetaIphi);

    scEtEB_PV->Fill(PVCollection_h->size(), scEB_Et);
    scHitEtEB->Fill(scEB_Et);
    //    scHitE_EB->Fill(scEB_E);
  }  // sc-EB loop

  //-------------------Compute scalar sum of reconstructed hit Et
  double RecHitEt_EB = 0;

  for (EcalRecHitCollection::const_iterator itr = RecHitsEB->begin(); itr != RecHitsEB->end(); ++itr) {
    // RecHitEt_EB +=itr->energy();

    GlobalPoint const &position = geom->getGeometry(itr->detid())->getPosition();
    RecHitEt_EB += itr->energy() * sin(position.theta());
  }  // EB Rec Hit

  recHitEtEB->Fill(RecHitEt_EB);
  recHitEtEB_PV->Fill(PVCollection_h->size(), RecHitEt_EB);

  //-------------------Compute scalar sum of reconstructed hit Et
  double RecHitEt_EE = 0;

  for (EcalRecHitCollection::const_iterator itr = RecHitsEE->begin(); itr != RecHitsEE->end(); ++itr) {
    GlobalPoint const &position = geom->getGeometry(itr->detid())->getPosition();
    RecHitEt_EE += itr->energy() * sin(position.theta());
  }  // EB Rec Hit

  recHitEtEE->Fill(RecHitEt_EE);
  recHitEtEE_PV->Fill(PVCollection_h->size(), RecHitEt_EE);
}

DEFINE_FWK_MODULE(EcalPileUpDepMonitor);