File indexing completed on 2024-04-06 12:33:14
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;
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
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
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
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
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);
