Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:11:07

0001 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0002 #include "DQMServices/Core/interface/DQMStore.h"
0003 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0004 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
0005 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0006 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/DetId/interface/DetId.h"
0009 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0010 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0011 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0012 #include "DataFormats/Math/interface/Vector3D.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0015 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0016 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementCluster.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementTrack.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
0021 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0022 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0023 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/PluginManager/interface/ModuleDef.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0031 
0032 #include <algorithm>
0033 #include <cmath>
0034 #include <iostream>
0035 #include <string>
0036 #include <utility>
0037 #include <vector>
0038 
0039 class PFClusterValidation : public DQMEDAnalyzer {
0040 public:
0041   PFClusterValidation(edm::ParameterSet const& conf);
0042   ~PFClusterValidation() override;
0043   void analyze(edm::Event const& e, edm::EventSetup const& c) override;
0044   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0045 
0046 private:
0047   static constexpr double partR2 = 0.3 * 0.3;  // dr cutoff (squared)
0048   static double sumEnergy(edm::Handle<reco::PFClusterCollection> const& pfClusters, double eta, double phi);
0049 
0050   edm::EDGetTokenT<edm::HepMCProduct> hepMCTok_;
0051   edm::EDGetTokenT<reco::PFClusterCollection> pfClusterECALTok_;
0052   edm::EDGetTokenT<reco::PFClusterCollection> pfClusterHCALTok_;
0053   edm::EDGetTokenT<reco::PFClusterCollection> pfClusterHOTok_;
0054   edm::EDGetTokenT<reco::PFClusterCollection> pfClusterHFTok_;
0055 
0056   MonitorElement* emean_vs_eta_E_;
0057   MonitorElement* emean_vs_eta_H_;
0058   MonitorElement* emean_vs_eta_EH_;
0059 
0060   MonitorElement* emean_vs_eta_HF_;
0061   MonitorElement* emean_vs_eta_HO_;
0062   MonitorElement* emean_vs_eta_EHF_;
0063   MonitorElement* emean_vs_eta_EHFO_;
0064 
0065   MonitorElement* ratio_Esummed_ECAL_0_;
0066   MonitorElement* ratio_Esummed_HCAL_0_;
0067   MonitorElement* ratio_Esummed_HO_0_;
0068 
0069   MonitorElement* ratio_Esummed_ECAL_1_;
0070   MonitorElement* ratio_Esummed_HCAL_1_;
0071   MonitorElement* ratio_Esummed_HO_1_;
0072 
0073   MonitorElement* ratio_Esummed_ECAL_2_;
0074   MonitorElement* ratio_Esummed_HCAL_2_;
0075   MonitorElement* ratio_Esummed_HO_2_;
0076 
0077   MonitorElement* ratio_Esummed_ECAL_3_;
0078   MonitorElement* ratio_Esummed_HCAL_3_;
0079   MonitorElement* ratio_Esummed_HO_3_;
0080 
0081   MonitorElement* ratio_Esummed_ECAL_4_;
0082   MonitorElement* ratio_Esummed_HCAL_4_;
0083   MonitorElement* ratio_Esummed_HO_4_;
0084 
0085   MonitorElement* ratio_Esummed_HF_5_;
0086   MonitorElement* ratio_Esummed_HF_6_;
0087 
0088   MonitorElement* ratio_Esummed_ECAL_HCAL_0_;
0089   MonitorElement* ratio_Esummed_ECAL_HCAL_HO_0_;
0090   MonitorElement* ratio_Esummed_ECAL_HCAL_1_;
0091   MonitorElement* ratio_Esummed_ECAL_HCAL_HO_1_;
0092   MonitorElement* ratio_Esummed_ECAL_HCAL_2_;
0093   MonitorElement* ratio_Esummed_ECAL_HCAL_HO_2_;
0094   MonitorElement* ratio_Esummed_ECAL_HCAL_3_;
0095   MonitorElement* ratio_Esummed_ECAL_HCAL_HO_3_;
0096   MonitorElement* ratio_Esummed_ECAL_HCAL_4_;
0097   MonitorElement* ratio_Esummed_ECAL_HCAL_HO_4_;
0098 
0099   MonitorElement* egen_MC_;
0100 };
0101 
0102 PFClusterValidation::PFClusterValidation(const edm::ParameterSet& conf) {
0103   hepMCTok_ = consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"));
0104 
0105   pfClusterECALTok_ =
0106       consumes<reco::PFClusterCollection>(conf.getUntrackedParameter<edm::InputTag>("pflowClusterECAL"));
0107   pfClusterHCALTok_ =
0108       consumes<reco::PFClusterCollection>(conf.getUntrackedParameter<edm::InputTag>("pflowClusterHCAL"));
0109   pfClusterHOTok_ = consumes<reco::PFClusterCollection>(conf.getUntrackedParameter<edm::InputTag>("pflowClusterHO"));
0110   pfClusterHFTok_ = consumes<reco::PFClusterCollection>(conf.getUntrackedParameter<edm::InputTag>("pflowClusterHF"));
0111 }
0112 
0113 PFClusterValidation::~PFClusterValidation() {}
0114 
0115 void PFClusterValidation::bookHistograms(DQMStore::IBooker& ibooker,
0116                                          edm::Run const& irun,
0117                                          edm::EventSetup const& isetup) {
0118   constexpr auto size = 100;
0119   char histo[size];
0120 
0121   constexpr double etaBinsOffset[] = {
0122       -5.191, -4.889, -4.716, -4.538, -4.363, -4.191, -4.013, -3.839, -3.664, -3.489, -3.314, -3.139, -2.964, -2.853,
0123       -2.65,  -2.5,   -2.322, -2.172, -2.043, -1.93,  -1.83,  -1.74,  -1.653, -1.566, -1.479, -1.392, -1.305, -1.218,
0124       -1.131, -1.044, -0.957, -0.879, -0.783, -0.696, -0.609, -0.522, -0.435, -0.348, -0.261, -0.174, -0.087, 0,
0125       0.087,  0.174,  0.261,  0.348,  0.435,  0.522,  0.609,  0.696,  0.783,  0.879,  0.957,  1.044,  1.131,  1.218,
0126       1.305,  1.392,  1.479,  1.566,  1.653,  1.74,   1.83,   1.93,   2.043,  2.172,  2.322,  2.5,    2.65,   2.853,
0127       2.964,  3.139,  3.314,  3.489,  3.664,  3.839,  4.013,  4.191,  4.363,  4.538,  4.716,  4.889,  5.191};
0128   constexpr int etaBins = std::size(etaBinsOffset) - 1;
0129 
0130   ibooker.setCurrentFolder("ParticleFlow/PFClusterV");
0131 
0132   // These are the single pion scan histos
0133 
0134   strncpy(histo, "emean_vs_eta_E", size);
0135   emean_vs_eta_E_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0136   strncpy(histo, "emean_vs_eta_H", size);
0137   emean_vs_eta_H_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0138   strncpy(histo, "emean_vs_eta_EH", size);
0139   emean_vs_eta_EH_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0140 
0141   strncpy(histo, "emean_vs_eta_HF", size);
0142   emean_vs_eta_HF_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0143   strncpy(histo, "emean_vs_eta_HO", size);
0144   emean_vs_eta_HO_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0145 
0146   strncpy(histo, "emean_vs_eta_EHF", size);
0147   emean_vs_eta_EHF_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0148   strncpy(histo, "emean_vs_eta_EHFO", size);
0149   emean_vs_eta_EHFO_ = ibooker.bookProfile(histo, histo, etaBins, etaBinsOffset, -100., 2000., " ");
0150 
0151   // 1D histos
0152 
0153   strncpy(histo, "Ratio_Esummed_ECAL_0", size);
0154   ratio_Esummed_ECAL_0_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0155   strncpy(histo, "Ratio_Esummed_HCAL_0", size);
0156   ratio_Esummed_HCAL_0_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0157   strncpy(histo, "Ratio_Esummed_HO_0", size);
0158   ratio_Esummed_HO_0_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0159   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_0", size);
0160   ratio_Esummed_ECAL_HCAL_0_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0161   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_HO_0", size);
0162   ratio_Esummed_ECAL_HCAL_HO_0_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0163 
0164   strncpy(histo, "Ratio_Esummed_ECAL_1", size);
0165   ratio_Esummed_ECAL_1_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0166   strncpy(histo, "Ratio_Esummed_HCAL_1", size);
0167   ratio_Esummed_HCAL_1_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0168   strncpy(histo, "Ratio_Esummed_HO_1", size);
0169   ratio_Esummed_HO_1_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0170   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_1", size);
0171   ratio_Esummed_ECAL_HCAL_1_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0172   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_HO_1", size);
0173   ratio_Esummed_ECAL_HCAL_HO_1_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0174 
0175   strncpy(histo, "Ratio_Esummed_ECAL_2", size);
0176   ratio_Esummed_ECAL_2_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0177   strncpy(histo, "Ratio_Esummed_HCAL_2", size);
0178   ratio_Esummed_HCAL_2_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0179   strncpy(histo, "Ratio_Esummed_HO_2", size);
0180   ratio_Esummed_HO_2_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0181   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_2", size);
0182   ratio_Esummed_ECAL_HCAL_2_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0183   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_HO_2", size);
0184   ratio_Esummed_ECAL_HCAL_HO_2_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0185 
0186   strncpy(histo, "Ratio_Esummed_ECAL_3", size);
0187   ratio_Esummed_ECAL_3_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0188   strncpy(histo, "Ratio_Esummed_HCAL_3", size);
0189   ratio_Esummed_HCAL_3_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0190   strncpy(histo, "Ratio_Esummed_HO_3", size);
0191   ratio_Esummed_HO_3_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0192   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_3", size);
0193   ratio_Esummed_ECAL_HCAL_3_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0194   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_HO_3", size);
0195   ratio_Esummed_ECAL_HCAL_HO_3_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0196 
0197   strncpy(histo, "Ratio_Esummed_ECAL_4", size);
0198   ratio_Esummed_ECAL_4_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0199   strncpy(histo, "Ratio_Esummed_HCAL_4", size);
0200   ratio_Esummed_HCAL_4_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0201   strncpy(histo, "Ratio_Esummed_HO_4", size);
0202   ratio_Esummed_HO_4_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0203   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_4", size);
0204   ratio_Esummed_ECAL_HCAL_4_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0205   strncpy(histo, "Ratio_Esummed_ECAL_HCAL_HO_4", size);
0206   ratio_Esummed_ECAL_HCAL_HO_4_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0207 
0208   strncpy(histo, "Ratio_Esummed_HF_5", size);
0209   ratio_Esummed_HF_5_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0210   strncpy(histo, "Ratio_Esummed_HF_6", size);
0211   ratio_Esummed_HF_6_ = ibooker.book1D(histo, histo, 50, 0., 5.);
0212 
0213   strncpy(histo, "Egen_MC", size);
0214   egen_MC_ = ibooker.book1D(histo, histo, 50, 0, 50);
0215 }
0216 
0217 void PFClusterValidation::analyze(edm::Event const& event, edm::EventSetup const& c) {
0218   double eta_MC = 0.;
0219   double phi_MC = 0.;
0220   double energy_MC = 0.;
0221 
0222   edm::Handle<edm::HepMCProduct> hepMC;
0223   event.getByToken(hepMCTok_, hepMC);
0224   if (not hepMC.isValid()) {
0225     edm::LogWarning("PFClusterValidation") << "HepMCProduct not found";
0226     return;
0227   }
0228 
0229   // MC particle with highest pt is taken as a direction reference
0230   double maxPt = -99999.;
0231   const HepMC::GenEvent* myGenEvent = hepMC->GetEvent();
0232   for (auto p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p) {
0233     double phip = (*p)->momentum().phi();
0234     double etap = (*p)->momentum().eta();
0235     double pt = (*p)->momentum().perp();
0236     double energy = (*p)->momentum().e();
0237     if (pt > maxPt) {
0238       maxPt = pt;
0239       energy_MC = energy;
0240       phi_MC = phip;
0241       eta_MC = etap;
0242     }
0243   }
0244 
0245   egen_MC_->Fill(energy_MC);
0246 
0247   edm::Handle<reco::PFClusterCollection> pfClusterECAL;
0248   event.getByToken(pfClusterECALTok_, pfClusterECAL);
0249 
0250   edm::Handle<reco::PFClusterCollection> pfClusterHCAL;
0251   event.getByToken(pfClusterHCALTok_, pfClusterHCAL);
0252 
0253   edm::Handle<reco::PFClusterCollection> pfClusterHO;
0254   event.getByToken(pfClusterHOTok_, pfClusterHO);
0255 
0256   edm::Handle<reco::PFClusterCollection> pfClusterHF;
0257   event.getByToken(pfClusterHFTok_, pfClusterHF);
0258 
0259   // sum the energy in a dR cone for each subsystem
0260   const double Econe = sumEnergy(pfClusterECAL, eta_MC, phi_MC);
0261   const double Hcone = sumEnergy(pfClusterHCAL, eta_MC, phi_MC);
0262   const double HOcone = sumEnergy(pfClusterHO, eta_MC, phi_MC);
0263   const double HFcone = sumEnergy(pfClusterHF, eta_MC, phi_MC);
0264 
0265   if (energy_MC > 0.) {
0266     if (std::abs(eta_MC) < 0.5) {
0267       ratio_Esummed_ECAL_0_->Fill(Econe / energy_MC);
0268       ratio_Esummed_HCAL_0_->Fill(Hcone / energy_MC);
0269       ratio_Esummed_HO_0_->Fill(HOcone / energy_MC);
0270       ratio_Esummed_ECAL_HCAL_0_->Fill((Econe + Hcone) / energy_MC);
0271       ratio_Esummed_ECAL_HCAL_HO_0_->Fill((Econe + Hcone + HOcone) / energy_MC);
0272     } else if (std::abs(eta_MC) < 1.3 && std::abs(eta_MC) > 0.5) {
0273       ratio_Esummed_ECAL_1_->Fill(Econe / energy_MC);
0274       ratio_Esummed_HCAL_1_->Fill(Hcone / energy_MC);
0275       ratio_Esummed_HO_1_->Fill(HOcone / energy_MC);
0276       ratio_Esummed_ECAL_HCAL_1_->Fill((Econe + Hcone) / energy_MC);
0277       ratio_Esummed_ECAL_HCAL_HO_1_->Fill((Econe + Hcone + HOcone) / energy_MC);
0278     } else if (std::abs(eta_MC) < 2.1 && std::abs(eta_MC) > 1.3) {
0279       ratio_Esummed_ECAL_2_->Fill(Econe / energy_MC);
0280       ratio_Esummed_HCAL_2_->Fill(Hcone / energy_MC);
0281       ratio_Esummed_HO_2_->Fill(HOcone / energy_MC);
0282       ratio_Esummed_ECAL_HCAL_2_->Fill((Econe + Hcone) / energy_MC);
0283       ratio_Esummed_ECAL_HCAL_HO_2_->Fill((Econe + Hcone + HOcone) / energy_MC);
0284     } else if (std::abs(eta_MC) < 2.5 && std::abs(eta_MC) > 2.1) {
0285       ratio_Esummed_ECAL_3_->Fill(Econe / energy_MC);
0286       ratio_Esummed_HCAL_3_->Fill(Hcone / energy_MC);
0287       ratio_Esummed_HO_3_->Fill(HOcone / energy_MC);
0288       ratio_Esummed_ECAL_HCAL_3_->Fill((Econe + Hcone) / energy_MC);
0289       ratio_Esummed_ECAL_HCAL_HO_3_->Fill((Econe + Hcone + HOcone) / energy_MC);
0290     } else if (2.5 < std::abs(eta_MC) && std::abs(eta_MC) < 3.0) {
0291       ratio_Esummed_ECAL_4_->Fill(Econe / energy_MC);
0292       ratio_Esummed_HCAL_4_->Fill(Hcone / energy_MC);
0293       ratio_Esummed_HO_4_->Fill(HOcone / energy_MC);
0294       ratio_Esummed_ECAL_HCAL_4_->Fill((Econe + Hcone) / energy_MC);
0295       ratio_Esummed_ECAL_HCAL_HO_4_->Fill((Econe + Hcone + HOcone) / energy_MC);
0296     } else if (3.0 < std::abs(eta_MC) && std::abs(eta_MC) < 4.0) {
0297       ratio_Esummed_HF_5_->Fill(HFcone / energy_MC);
0298     } else if (4.0 < std::abs(eta_MC) && std::abs(eta_MC) < 5.0) {
0299       ratio_Esummed_HF_6_->Fill(HFcone / energy_MC);
0300     }
0301   }
0302 
0303   emean_vs_eta_E_->Fill(eta_MC, Econe);
0304   emean_vs_eta_H_->Fill(eta_MC, Hcone);
0305   emean_vs_eta_EH_->Fill(eta_MC, Econe + Hcone);
0306   emean_vs_eta_HF_->Fill(eta_MC, HFcone);
0307   emean_vs_eta_HO_->Fill(eta_MC, HOcone);
0308   emean_vs_eta_EHF_->Fill(eta_MC, Econe + Hcone + HFcone);
0309   emean_vs_eta_EHFO_->Fill(eta_MC, Econe + Hcone + HFcone + HOcone);
0310 }
0311 
0312 double PFClusterValidation::sumEnergy(edm::Handle<reco::PFClusterCollection> const& pfClusters,
0313                                       double eta,
0314                                       double phi) {
0315   if (not pfClusters.isValid())
0316     return 0.;
0317 
0318   double sum = 0.;
0319   for (auto pf = pfClusters->begin(); pf != pfClusters->end(); ++pf) {
0320     if (reco::deltaR2(eta, phi, pf->eta(), pf->phi()) < partR2) {
0321       sum += pf->energy();
0322     }
0323   }
0324 
0325   return sum;
0326 }
0327 
0328 #include "FWCore/Framework/interface/MakerMacros.h"
0329 DEFINE_FWK_MODULE(PFClusterValidation);