Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:09:20

0001 /*
0002  * \file EcalPileUpDepMonitor.cc
0003  * \author Ben Carlson - CMU
0004  * Last Update:
0005  *
0006  */
0007 
0008 #include "DQMOffline/Ecal/interface/EcalPileUpDepMonitor.h"
0009 
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0012 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0013 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0014 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0015 #include "DataFormats/VertexReco/interface/Vertex.h"
0016 
0017 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0018 
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 
0025 #include <cmath>
0026 
0027 // Framework
0028 
0029 EcalPileUpDepMonitor::EcalPileUpDepMonitor(const edm::ParameterSet &ps) : geomH(esConsumes()), caloTop(esConsumes()) {
0030   VertexCollection_ = consumes<reco::VertexCollection>(ps.getParameter<edm::InputTag>("VertexCollection"));
0031 
0032   if (ps.existsAs<edm::InputTag>("basicClusterCollection") &&
0033       !ps.getParameter<edm::InputTag>("basicClusterCollection").label().empty())
0034     basicClusterCollection_ =
0035         consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection"));
0036   else {
0037     basicClusterCollection_EE_ =
0038         consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection_EE"));
0039     basicClusterCollection_EB_ =
0040         consumes<edm::View<reco::CaloCluster>>(ps.getParameter<edm::InputTag>("basicClusterCollection_EB"));
0041   }
0042   superClusterCollection_EB_ =
0043       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("superClusterCollection_EB"));
0044   superClusterCollection_EE_ =
0045       consumes<reco::SuperClusterCollection>(ps.getParameter<edm::InputTag>("superClusterCollection_EE"));
0046 
0047   RecHitCollection_EB_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("RecHitCollection_EB"));
0048   RecHitCollection_EE_ = consumes<EcalRecHitCollection>(ps.getParameter<edm::InputTag>("RecHitCollection_EE"));
0049   EleTag_ = consumes<reco::GsfElectronCollection>(ps.getParameter<edm::InputTag>("EleTag"));
0050 }
0051 
0052 EcalPileUpDepMonitor::~EcalPileUpDepMonitor() {}
0053 
0054 void EcalPileUpDepMonitor::bookHistograms(DQMStore::IBooker &ibooker,
0055                                           edm::Run const &,
0056                                           edm::EventSetup const &eventSetup) {
0057   ibooker.cd();
0058   ibooker.setCurrentFolder("Ecal/EcalPileUpDepMonitor");
0059 
0060   std::string prof_name = "bcEB_PV";
0061   std::string title = "Basic clusters EB vs. PV";
0062   bcEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0063   bcEB_PV->setAxisTitle("N_{pv}", 1);
0064   bcEB_PV->setAxisTitle("Basic Clusters", 2);
0065 
0066   prof_name = "bcEE_PV";
0067   title = "Basic Clusters EE vs. PV";
0068   bcEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0069   bcEE_PV->setAxisTitle("N_{pv}", 1);
0070   bcEE_PV->setAxisTitle("Basic Clusters", 2);
0071 
0072   prof_name = "scEB_PV";
0073   title = "Super Clusters EB vs. PV";
0074   scEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0075   scEB_PV->setAxisTitle("N_{pv}", 1);
0076   scEB_PV->setAxisTitle("Super Clusters", 2);
0077 
0078   prof_name = "scEE_PV";
0079   title = "Super Clusters EE vs. PV";
0080   scEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0081   scEE_PV->setAxisTitle("N_{pv}", 1);
0082   scEE_PV->setAxisTitle("Super Clusters", 2);
0083 
0084   prof_name = "scEtEB_PV";
0085   title = "Super Clusters Et EB vs. PV";
0086   scEtEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0087   scEtEB_PV->setAxisTitle("N_{pv}", 1);
0088   scEtEB_PV->setAxisTitle("Super Cluster E_{T} [GeV]", 2);
0089 
0090   prof_name = "scEtEE_PV";
0091   title = "Super Clusters Et EE vs. PV";
0092   scEtEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0093   scEtEE_PV->setAxisTitle("N_{pv}", 1);
0094   scEtEE_PV->setAxisTitle("Super Cluster E_{T} [GeV]", 2);
0095 
0096   prof_name = "recHitEtEB_PV";
0097   title = "Reconstructed Hit Et EB vs. PV";
0098   recHitEtEB_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0099   recHitEtEB_PV->setAxisTitle("N_{pv}", 1);
0100   recHitEtEB_PV->setAxisTitle("Reconstructed hit E_{T} [GeV]", 2);
0101 
0102   prof_name = "recHitEtEE_PV";
0103   title = "Reconstructed Hit Et EE vs. PV";
0104   recHitEtEE_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350.);
0105   recHitEtEE_PV->setAxisTitle("N_{pv}", 1);
0106   recHitEtEE_PV->setAxisTitle("Reconstructed hit E_{T} [GeV]", 2);
0107 
0108   prof_name = "emIso_PV";
0109   title = "EM Isolation vs. PV";
0110   emIso_PV = ibooker.bookProfile(prof_name, title, 50, 0., 50., 50, 0., 350);
0111   emIso_PV->setAxisTitle("N_{pv}", 1);
0112   emIso_PV->setAxisTitle("EM_{Isolation} [GeV]", 2);
0113 
0114   prof_name = "emIso";
0115   title = "EM Isolation";
0116   emIso = ibooker.book1D(prof_name, title, 50, 0, 50);
0117   emIso->setAxisTitle("EM_{Isolation} [GeV]", 1);
0118   emIso->setAxisTitle("Events", 2);
0119 
0120   prof_name = "scHitEtEB";
0121   title = "Super Cluster Hit Et EB";
0122   scHitEtEB = ibooker.book1D(prof_name, title, 50, 0, 100);
0123   scHitEtEB->setAxisTitle("super cluster hit E_{T} [GeV]", 1);
0124   scHitEtEB->setAxisTitle("Events", 2);
0125 
0126   prof_name = "scHitEtEE";
0127   title = "Super Cluster Hit Et EE";
0128   scHitEtEE = ibooker.book1D(prof_name, title, 50, 0, 100);
0129   scHitEtEE->setAxisTitle("super cluster hit E_{T} [GeV]", 1);
0130   scHitEtEE->setAxisTitle("Events", 2);
0131 
0132   //   prof_name="scHitE_EB";
0133   //   title="Super Cluster Hit E EB";
0134   //   scHitE_EB=ibooker.book1D(prof_name,title,50,0,100);
0135   //   scHitE_EB->setAxisTitle("super cluster hit E [GeV]",1);
0136   //   scHitE_EB->setAxisTitle("Events",2);
0137 
0138   //   prof_name="scHitE_EE";
0139   //   title="Super Cluster Hit E EE";
0140   //   scHitE_EE=ibooker.book1D(prof_name,title,50,0,100);
0141   //   scHitE_EE->setAxisTitle("super cluster hit E [GeV]",1);
0142   //   scHitE_EE->setAxisTitle("Events",2);
0143 
0144   // SC eta
0145   //   prof_name="scEta_EB";
0146   //   title="Super Cluster #eta EB";
0147   //   scEta_EB=ibooker.book1D(prof_name,title,50,-6,6);
0148   //   scEta_EB->setAxisTitle("#eta",1);
0149   //   scEta_EB->setAxisTitle("Events",2);
0150 
0151   //   prof_name="scEta_EE";
0152   //   title="Super Cluster #eta EE";
0153   //   scEta_EE=ibooker.book1D(prof_name,title,50,-6,6);
0154   //   scEta_EE->setAxisTitle("#eta",1);
0155   //   scEta_EE->setAxisTitle("Events",2);
0156 
0157   // SC phi
0158   //   prof_name="scPhi_EB";
0159   //   title="Super Cluster #phi EB";
0160   //   scPhi_EB=ibooker.book1D(prof_name,title,50,-3.14,3.14);
0161   //   scPhi_EB->setAxisTitle("super cluster #phi",1);
0162   //   scPhi_EB->setAxisTitle("Events",2);
0163 
0164   //   prof_name="scPhi_EE";
0165   //   title="Super Cluster #phi EE";
0166   //   scPhi_EE=ibooker.book1D(prof_name,title,50,-3.14,3.14);
0167   //   scPhi_EE->setAxisTitle("super cluster #phi",1);
0168   //   scPhi_EE->setAxisTitle("Events",2);
0169 
0170   // sc sigma eta eta / eta phi
0171 
0172   prof_name = "scSigmaIetaIeta_EB";
0173   title = "Super Cluster sigmaIetaIeta EB";
0174   scSigmaIetaIeta_EB = ibooker.book1D(prof_name, title, 50, 0, 0.03);
0175   scSigmaIetaIeta_EB->setAxisTitle("#sigma_{i#etai#eta}", 1);
0176   scSigmaIetaIeta_EB->setAxisTitle("Events", 2);
0177 
0178   prof_name = "scSigmaIetaIeta_EE";
0179   title = "Super Cluster sigmaIetaIeta EE";
0180   scSigmaIetaIeta_EE = ibooker.book1D(prof_name, title, 50, 0, 0.1);
0181   scSigmaIetaIeta_EE->setAxisTitle("#sigma_{i#etai#eta}", 1);
0182   scSigmaIetaIeta_EE->setAxisTitle("Events", 2);
0183 
0184   // phi
0185   prof_name = "scSigmaIetaIphi_EB";
0186   title = "Super Cluster sigmaIetaIphi EB";
0187   scSigmaIetaIphi_EB = ibooker.book1D(prof_name, title, 50, -5.e-4, 5.e-4);
0188   scSigmaIetaIphi_EB->setAxisTitle("#sigma_{i#etai#phi}", 1);
0189   scSigmaIetaIphi_EB->setAxisTitle("Events", 2);
0190 
0191   prof_name = "scSigmaIetaIphi_EE";
0192   title = "Super Cluster sigmaIetaIphi EE";
0193   scSigmaIetaIphi_EE = ibooker.book1D(prof_name, title, 50, -2.5e-3, 2.5e-3);
0194   scSigmaIetaIphi_EE->setAxisTitle("#sigma_{i#etai#phi}", 1);
0195   scSigmaIetaIphi_EE->setAxisTitle("Events", 2);
0196 
0197   // R9
0198   prof_name = "r9_EB";
0199   title = "r9 EB";
0200   r9_EB = ibooker.book1D(prof_name, title, 50, 0, 1.5);
0201   r9_EB->setAxisTitle("R_{9}", 1);
0202   r9_EB->setAxisTitle("Events", 2);
0203 
0204   prof_name = "r9_EE";
0205   title = "r9 EE";
0206   r9_EE = ibooker.book1D(prof_name, title, 50, 0, 1.5);
0207   r9_EE->setAxisTitle("R_{9}", 1);
0208   r9_EE->setAxisTitle("Events", 2);
0209 
0210   // Rec Hit
0211 
0212   prof_name = "recHitEtEB";
0213   title = "RecHit Et EB";
0214   recHitEtEB = ibooker.book1D(prof_name, title, 50, 0, 400);
0215   recHitEtEB->setAxisTitle("Reconstructed Hit E_{T} [GeV]", 1);
0216   recHitEtEB->setAxisTitle("Events", 2);
0217 
0218   prof_name = "recHitEtEE";
0219   title = "RecHit Et EE";
0220   recHitEtEE = ibooker.book1D(prof_name, title, 50, 0, 400);
0221   recHitEtEE->setAxisTitle("Reconstructed Hit E_{T} [GeV]", 1);
0222   recHitEtEE->setAxisTitle("Events", 2);
0223 }
0224 
0225 void EcalPileUpDepMonitor::analyze(const edm::Event &e, const edm::EventSetup &es) {
0226   const CaloGeometry *geom = &es.getData(geomH);
0227   // Vertex collection:
0228   //-----------------------------------------
0229   edm::Handle<reco::VertexCollection> PVCollection_h;
0230   e.getByToken(VertexCollection_, PVCollection_h);
0231   if (!PVCollection_h.isValid()) {
0232     edm::LogWarning("VertexCollection") << "VertexCollection not found";
0233   }
0234   //-----------------gsfElectrons -------------------------
0235   edm::Handle<reco::GsfElectronCollection> electronCollection_h;
0236   e.getByToken(EleTag_, electronCollection_h);
0237   if (!electronCollection_h.isValid()) {
0238     edm::LogWarning("EBRecoSummary") << "Electrons not found";
0239   }
0240 
0241   if (basicClusterCollection_.isUninitialized()) {
0242     //----------------- Basic Cluster Collection Ecal Barrel  ---------
0243     edm::Handle<edm::View<reco::CaloCluster>> basicClusters_EB_h;
0244     e.getByToken(basicClusterCollection_EB_, basicClusters_EB_h);
0245     if (!basicClusters_EB_h.isValid()) {
0246       edm::LogWarning("EBRecoSummary") << "basicClusters_EB_h not found";
0247     }
0248 
0249     bcEB_PV->Fill(PVCollection_h->size(), basicClusters_EB_h->size());
0250 
0251     //----------------- Basic Cluster Collection Ecal Endcal  ---------
0252 
0253     edm::Handle<edm::View<reco::CaloCluster>> basicClusters_EE_h;
0254     e.getByToken(basicClusterCollection_EE_, basicClusters_EE_h);
0255     if (!basicClusters_EE_h.isValid()) {
0256       edm::LogWarning("EERecoSummary") << "basicClusters_EE_h not found";
0257     }
0258 
0259     bcEE_PV->Fill(PVCollection_h->size(), basicClusters_EE_h->size());
0260   } else {
0261     //----------------- Basic Cluster Collection Ecal Barrel  ---------
0262     edm::Handle<edm::View<reco::CaloCluster>> basicClusters_h;
0263     e.getByToken(basicClusterCollection_, basicClusters_h);
0264     if (!basicClusters_h.isValid()) {
0265       edm::LogWarning("EBRecoSummary") << "basicClusters_h not found";
0266     }
0267 
0268     int nBarrel(0);
0269     int nEndcap(0);
0270     for (edm::View<reco::CaloCluster>::const_iterator bcItr(basicClusters_h->begin()); bcItr != basicClusters_h->end();
0271          ++bcItr) {
0272       if (bcItr->caloID().detector(reco::CaloID::DET_ECAL_BARREL))
0273         ++nBarrel;
0274       if (bcItr->caloID().detector(reco::CaloID::DET_ECAL_ENDCAP))
0275         ++nEndcap;
0276     }
0277 
0278     bcEB_PV->Fill(PVCollection_h->size(), nBarrel);
0279     bcEE_PV->Fill(PVCollection_h->size(), nEndcap);
0280   }
0281 
0282   //----------------- Reconstructed Hit Ecal barrel
0283 
0284   edm::Handle<EcalRecHitCollection> RecHitsEB;
0285   e.getByToken(RecHitCollection_EB_, RecHitsEB);
0286   if (!RecHitsEB.isValid()) {
0287     edm::LogWarning("EBRecoSummary") << "RecHitsEB not found";
0288   }
0289 
0290   //----------------- Reconstructed Hit Ecal Endcap
0291 
0292   edm::Handle<EcalRecHitCollection> RecHitsEE;
0293   e.getByToken(RecHitCollection_EE_, RecHitsEE);
0294   if (!RecHitsEE.isValid()) {
0295     edm::LogWarning("EBRecoSummary") << "RecHitsEB not found";
0296   }
0297   //----------------- Super Cluster Collection Ecal Endcap  ---------
0298 
0299   edm::Handle<reco::SuperClusterCollection> superClusters_EE_h;
0300   e.getByToken(superClusterCollection_EE_, superClusters_EE_h);
0301   if (!superClusters_EE_h.isValid()) {
0302     edm::LogWarning("EERecoSummary") << "superClusters_EE_h not found";
0303   }
0304 
0305   //--------- Fill Isolation -----------------
0306 
0307   if (electronCollection_h.isValid()) {
0308     for (reco::GsfElectronCollection::const_iterator recoElectron = electronCollection_h->begin();
0309          recoElectron != electronCollection_h->end();
0310          recoElectron++) {
0311       double IsoEcal = recoElectron->dr03EcalRecHitSumEt();  /// recoElectron->et()
0312       emIso_PV->Fill(PVCollection_h->size(), IsoEcal);
0313       emIso->Fill(IsoEcal);
0314     }
0315   }
0316 
0317   // fill super clusters EE
0318   scEE_PV->Fill(PVCollection_h->size(), superClusters_EE_h->size());
0319 
0320   for (reco::SuperClusterCollection::const_iterator itSC = superClusters_EE_h->begin();
0321        itSC != superClusters_EE_h->end();
0322        ++itSC) {
0323     double scEE_Et = itSC->energy() * sin(2. * atan(exp(-itSC->position().eta())));
0324     //    double scEE_E=itSC->energy();
0325 
0326     // fill super cluster endcap eta/phi
0327     //     scEta_EE->Fill(itSC->position().eta());
0328     //     scPhi_EE->Fill(itSC->position().phi());
0329 
0330     // get sigma eta_eta etc
0331 
0332     CaloTopology const *p_topology = &es.getData(caloTop);  // get calo topology
0333     const EcalRecHitCollection *eeRecHits = RecHitsEE.product();
0334 
0335     reco::BasicCluster const &seedCluster(*itSC->seed());
0336     const auto &cov = EcalClusterTools::localCovariances(seedCluster, eeRecHits, p_topology);
0337     float sigmaIetaIeta = std::sqrt(cov[0]);
0338     float sigmaIetaIphi = cov[1];
0339 
0340     float e3x3 = EcalClusterTools::e3x3(seedCluster, eeRecHits, p_topology);
0341     float r9 = e3x3 / itSC->rawEnergy();
0342 
0343     r9_EE->Fill(r9);
0344     scSigmaIetaIeta_EE->Fill(sigmaIetaIeta);
0345     scSigmaIetaIphi_EE->Fill(sigmaIetaIphi);
0346 
0347     // std::cout  << " sigmaIetaIeta: " << sigmaIetaIeta << std::endl;
0348     scEtEE_PV->Fill(PVCollection_h->size(), scEE_Et);
0349     scHitEtEE->Fill(scEE_Et);  // super cluster Et historam
0350     //    scHitE_EE->Fill(scEE_E); //super cluster energy histogram
0351 
0352   }  // sc-EE loop
0353 
0354   //----------------- Super Cluster Collection Ecal Barrel  ---------
0355 
0356   edm::Handle<reco::SuperClusterCollection> superClusters_EB_h;
0357   e.getByToken(superClusterCollection_EB_, superClusters_EB_h);
0358   if (!superClusters_EB_h.isValid()) {
0359     edm::LogWarning("EBRecoSummary") << "superClusters_EB_h not found";
0360   }
0361   scEB_PV->Fill(PVCollection_h->size(), superClusters_EB_h->size());
0362 
0363   for (reco::SuperClusterCollection::const_iterator itSC = superClusters_EB_h->begin();
0364        itSC != superClusters_EB_h->end();
0365        ++itSC) {
0366     double scEB_Et = itSC->energy() * sin(2. * atan(exp(-itSC->position().eta())));  // super cluster transverse energy
0367     //    double scEB_E= itSC->energy(); // super cluster energy
0368 
0369     // fill super cluster Barrel eta/phi
0370     //     scEta_EB->Fill(itSC->position().eta()); //super cluster eta
0371     //     scPhi_EB->Fill(itSC->position().phi()); // super cluster phi
0372 
0373     // sigma ietaieta etc
0374 
0375     CaloTopology const *p_topology = &es.getData(caloTop);  // get calo topology
0376     const EcalRecHitCollection *ebRecHits = RecHitsEB.product();
0377 
0378     reco::BasicCluster const &seedCluster(*itSC->seed());
0379     const auto &cov = EcalClusterTools::localCovariances(seedCluster, ebRecHits, p_topology);
0380     float sigmaIetaIeta = std::sqrt(cov[0]);
0381     float sigmaIetaIphi = cov[1];
0382 
0383     float e3x3 = EcalClusterTools::e3x3(seedCluster, ebRecHits, p_topology);
0384     float r9 = e3x3 / itSC->rawEnergy();
0385 
0386     r9_EB->Fill(r9);
0387     scSigmaIetaIeta_EB->Fill(sigmaIetaIeta);
0388     scSigmaIetaIphi_EB->Fill(sigmaIetaIphi);
0389 
0390     scEtEB_PV->Fill(PVCollection_h->size(), scEB_Et);
0391     scHitEtEB->Fill(scEB_Et);
0392     //    scHitE_EB->Fill(scEB_E);
0393   }  // sc-EB loop
0394 
0395   //-------------------Compute scalar sum of reconstructed hit Et
0396   double RecHitEt_EB = 0;
0397 
0398   for (EcalRecHitCollection::const_iterator itr = RecHitsEB->begin(); itr != RecHitsEB->end(); ++itr) {
0399     // RecHitEt_EB +=itr->energy();
0400 
0401     GlobalPoint const &position = geom->getGeometry(itr->detid())->getPosition();
0402     RecHitEt_EB += itr->energy() * sin(position.theta());
0403   }  // EB Rec Hit
0404 
0405   recHitEtEB->Fill(RecHitEt_EB);
0406   recHitEtEB_PV->Fill(PVCollection_h->size(), RecHitEt_EB);
0407 
0408   //-------------------Compute scalar sum of reconstructed hit Et
0409   double RecHitEt_EE = 0;
0410 
0411   for (EcalRecHitCollection::const_iterator itr = RecHitsEE->begin(); itr != RecHitsEE->end(); ++itr) {
0412     GlobalPoint const &position = geom->getGeometry(itr->detid())->getPosition();
0413     RecHitEt_EE += itr->energy() * sin(position.theta());
0414   }  // EB Rec Hit
0415 
0416   recHitEtEE->Fill(RecHitEt_EE);
0417   recHitEtEE_PV->Fill(PVCollection_h->size(), RecHitEt_EE);
0418 }
0419 
0420 DEFINE_FWK_MODULE(EcalPileUpDepMonitor);