EcalRecHits/src/EcalRecHitsValidation.cc

0001 /*
0002  * \file EcalRecHitsValidation.cc
0003  *
0004  * \author C. Rovelli
0005  *
0006  */
0007
0008 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
0009 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0012 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0013 #include <DataFormats/EcalDetId/interface/EBDetId.h>
0014 #include <DataFormats/EcalDetId/interface/EEDetId.h>
0015 #include <DataFormats/EcalDetId/interface/ESDetId.h>
0016 #include <Validation/EcalRecHits/interface/EcalRecHitsValidation.h>
0017 #include <string>
0018
0019 using namespace cms;
0020 using namespace edm;
0021 using namespace std;
0022
0023 EcalRecHitsValidation::EcalRecHitsValidation(const ParameterSet &ps)
0024     : pAgc(esConsumes()), pEcsToken(esConsumes()), pttMapToken(esConsumes()) {
0025   // ----------------------
0026   HepMCLabel = ps.getParameter<std::string>("moduleLabelMC");
0027   hitsProducer_ = ps.getParameter<std::string>("hitsProducer");
0028   EBrechitCollection_ = ps.getParameter<edm::InputTag>("EBrechitCollection");
0029   EErechitCollection_ = ps.getParameter<edm::InputTag>("EErechitCollection");
0030   ESrechitCollection_ = ps.getParameter<edm::InputTag>("ESrechitCollection");
0031   EBuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EBuncalibrechitCollection");
0032   EEuncalibrechitCollection_ = ps.getParameter<edm::InputTag>("EEuncalibrechitCollection");
0033   // switch to allow disabling endcaps for phase 2
0034   enableEndcaps_ = ps.getUntrackedParameter<bool>("enableEndcaps", true);
0035
0036   // fix for consumes
0037   HepMCLabel_Token_ = consumes<HepMCProduct>(ps.getParameter<std::string>("moduleLabelMC"));
0038   EBrechitCollection_Token_ = consumes<EBRecHitCollection>(ps.getParameter<edm::InputTag>("EBrechitCollection"));
0039   EBuncalibrechitCollection_Token_ =
0040       consumes<EBUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EBuncalibrechitCollection"));
0041
0042   EBHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0043       edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEB")));
0044   if (enableEndcaps_) {
0045     EErechitCollection_Token_ = consumes<EERecHitCollection>(ps.getParameter<edm::InputTag>("EErechitCollection"));
0046     ESrechitCollection_Token_ = consumes<ESRecHitCollection>(ps.getParameter<edm::InputTag>("ESrechitCollection"));
0047     EEuncalibrechitCollection_Token_ =
0048         consumes<EEUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EEuncalibrechitCollection"));
0049     EEHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0050         edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsEE")));
0051     ESHits_Token_ = consumes<CrossingFrame<PCaloHit>>(
0052         edm::InputTag(std::string("mix"), ps.getParameter<std::string>("hitsProducer") + std::string("EcalHitsES")));
0053   }
0054
0055   // ----------------------
0056   // DQM ROOT output
0057   outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
0058
0059   if (!outputFile_.empty()) {
0060     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will be saved to '" << outputFile_.c_str() << "'";
0061   } else {
0062     LogInfo("OutputInfo") << " Ecal RecHits Task histograms will NOT be saved";
0063   }
0064
0065   // ----------------------
0066   // verbosity switch
0067   verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
0068
0069   // ----------------------
0070   meGunEnergy_ = nullptr;
0071   meGunEta_ = nullptr;
0072   meGunPhi_ = nullptr;
0073   meEBRecHitSimHitRatio_ = nullptr;
0074   meEERecHitSimHitRatio_ = nullptr;
0075   meESRecHitSimHitRatio_ = nullptr;
0076   meEBRecHitSimHitRatio1011_ = nullptr;
0077   meEERecHitSimHitRatio1011_ = nullptr;
0078   meEBRecHitSimHitRatio12_ = nullptr;
0079   meEERecHitSimHitRatio12_ = nullptr;
0080   meEBRecHitSimHitRatio13_ = nullptr;
0081   meEERecHitSimHitRatio13_ = nullptr;
0082   meEBRecHitSimHitRatioGt35_ = nullptr;
0083   meEERecHitSimHitRatioGt35_ = nullptr;
0084   meEBUnRecHitSimHitRatio_ = nullptr;
0085   meEEUnRecHitSimHitRatio_ = nullptr;
0086   meEBUnRecHitSimHitRatioGt35_ = nullptr;
0087   meEEUnRecHitSimHitRatioGt35_ = nullptr;
0088   meEBe5x5_ = nullptr;
0089   meEBe5x5OverSimHits_ = nullptr;
0090   meEBe5x5OverGun_ = nullptr;
0091   meEEe5x5_ = nullptr;
0092   meEEe5x5OverSimHits_ = nullptr;
0093   meEEe5x5OverGun_ = nullptr;
0094
0095   meEBRecHitLog10Energy_ = nullptr;
0096   meEERecHitLog10Energy_ = nullptr;
0097   meESRecHitLog10Energy_ = nullptr;
0098   meEBRecHitLog10EnergyContr_ = nullptr;
0099   meEERecHitLog10EnergyContr_ = nullptr;
0100   meESRecHitLog10EnergyContr_ = nullptr;
0101   meEBRecHitLog10Energy5x5Contr_ = nullptr;
0102   meEERecHitLog10Energy5x5Contr_ = nullptr;
0103
0104   meEBRecHitsOccupancyFlag5_6_ = nullptr;
0105   meEBRecHitsOccupancyFlag8_9_ = nullptr;
0106   meEERecHitsOccupancyPlusFlag5_6_ = nullptr;
0107   meEERecHitsOccupancyMinusFlag5_6_ = nullptr;
0108   meEERecHitsOccupancyPlusFlag8_9_ = nullptr;
0109   meEERecHitsOccupancyMinusFlag8_9_ = nullptr;
0110
0111   meEBRecHitFlags_ = nullptr;
0112   meEBRecHitSimHitvsSimHitFlag5_6_ = nullptr;
0113   meEBRecHitSimHitFlag6_ = nullptr;
0114   meEBRecHitSimHitFlag7_ = nullptr;
0115   meEB5x5RecHitSimHitvsSimHitFlag8_ = nullptr;
0116
0117   meEERecHitFlags_ = nullptr;
0118   meEERecHitSimHitvsSimHitFlag5_6_ = nullptr;
0119   meEERecHitSimHitFlag6_ = nullptr;
0120   meEERecHitSimHitFlag7_ = nullptr;
0121 }
0122
0123 EcalRecHitsValidation::~EcalRecHitsValidation() {}
0124
0125 void EcalRecHitsValidation::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &, edm::EventSetup const &) {
0126   std::string histo;
0127
0128   ibooker.setCurrentFolder("EcalRecHitsV/EcalRecHitsTask");
0129
0130   histo = "EcalRecHitsTask Gun Momentum";
0131   meGunEnergy_ = ibooker.book1D(histo.c_str(), histo.c_str(), 100, 0., 1000.);
0132
0133   histo = "EcalRecHitsTask Gun Eta";
0134   meGunEta_ = ibooker.book1D(histo.c_str(), histo.c_str(), 700, -3.5, 3.5);
0135
0136   histo = "EcalRecHitsTask Gun Phi";
0137   meGunPhi_ = ibooker.book1D(histo.c_str(), histo.c_str(), 360, 0., 360.);
0138
0139   histo = "EcalRecHitsTask Barrel RecSimHit Ratio";
0140   meEBRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0141
0142   histo = "EcalRecHitsTask Barrel RecSimHit Ratio gt 3p5 GeV";
0143   meEBRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0144
0145   histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio";
0146   meEBUnRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0147
0148   histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=10 11";
0149   meEBRecHitSimHitRatio1011_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0150
0151   histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=12";
0152   meEBRecHitSimHitRatio12_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0153
0154   histo = "EcalRecHitsTask Barrel RecSimHit Ratio Channel Status=13";
0155   meEBRecHitSimHitRatio13_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0156
0157   histo = "EcalRecHitsTask Barrel Unc RecSimHit Ratio gt 3p5 GeV";
0158   meEBUnRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0159
0160   histo = "EcalRecHitsTask Barrel Rec E5x5";
0161   meEBe5x5_ = ibooker.book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
0162
0163   histo = "EcalRecHitsTask Barrel Rec E5x5 over Sim E5x5";
0164   meEBe5x5OverSimHits_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0165
0166   histo = "EcalRecHitsTask Barrel Rec E5x5 over gun energy";
0167   meEBe5x5OverGun_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0168
0169   meEBRecHitLog10Energy_ =
0170       ibooker.book1D("EcalRecHitsTask Barrel Log10 Energy", "EcalRecHitsTask Barrel Log10 Energy", 90, -5., 4.);
0171   meEBRecHitLog10EnergyContr_ = ibooker.bookProfile("EcalRecHits Barrel Log10En vs Hit Contribution",
0172                                                     "EcalRecHits Barrel Log10En vs Hit Contribution",
0173                                                     90,
0174                                                     -5.,
0175                                                     4.,
0176                                                     100,
0177                                                     0.,
0178                                                     1.);
0179   meEBRecHitLog10Energy5x5Contr_ = ibooker.bookProfile("EcalRecHits Barrel Log10En5x5 vs Hit Contribution",
0180                                                        "EcalRecHits Barrel Log10En5x5 vs Hit Contribution",
0181                                                        90,
0182                                                        -5.,
0183                                                        4.,
0184                                                        100,
0185                                                        0.,
0186                                                        1.);
0187
0188   histo = "EB Occupancy Flag=5 6";
0189   meEBRecHitsOccupancyFlag5_6_ = ibooker.book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
0190   histo = "EB Occupancy Flag=8 9";
0191   meEBRecHitsOccupancyFlag8_9_ = ibooker.book2D(histo, histo, 170, -85., 85., 360, 0., 360.);
0192
0193   histo = "EcalRecHitsTask Barrel Reco Flags";
0194   meEBRecHitFlags_ = ibooker.book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
0195   histo = "EcalRecHitsTask Barrel RecSimHit Ratio vs SimHit Flag=5 6";
0196   meEBRecHitSimHitvsSimHitFlag5_6_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400.);
0197   histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=6";
0198   meEBRecHitSimHitFlag6_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0199   histo = "EcalRecHitsTask Barrel RecSimHit Ratio Flag=7";
0200   meEBRecHitSimHitFlag7_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0201   histo = "EcalRecHitsTask Barrel 5x5 RecSimHit Ratio vs SimHit Flag=8";
0202   meEB5x5RecHitSimHitvsSimHitFlag8_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400.);
0203
0204   if (enableEndcaps_) {
0205     histo = "EcalRecHitsTask Preshower RecSimHit Ratio";
0206     meESRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0207
0208     histo = "EcalRecHitsTask Endcap RecSimHit Ratio";
0209     meEERecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0210
0211     histo = "EcalRecHitsTask Endcap RecSimHit Ratio gt 3p5 GeV";
0212     meEERecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0213
0214     histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio";
0215     meEEUnRecHitSimHitRatio_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0216
0217     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=10 11";
0218     meEERecHitSimHitRatio1011_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0219
0220     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=12";
0221     meEERecHitSimHitRatio12_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0222
0223     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Channel Status=13";
0224     meEERecHitSimHitRatio13_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0225
0226     histo = "EcalRecHitsTask Endcap Unc RecSimHit Ratio gt 3p5 GeV";
0227     meEEUnRecHitSimHitRatioGt35_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0228
0229     histo = "EcalRecHitsTask Endcap Rec E5x5";
0230     meEEe5x5_ = ibooker.book1D(histo.c_str(), histo.c_str(), 4000, 0., 400.);
0231
0232     histo = "EcalRecHitsTask Endcap Rec E5x5 over Sim E5x5";
0233     meEEe5x5OverSimHits_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0234
0235     histo = "EcalRecHitsTask Endcap Rec E5x5 over gun energy";
0236     meEEe5x5OverGun_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0.9, 1.1);
0237
0238     meEERecHitLog10Energy_ =
0239         ibooker.book1D("EcalRecHitsTask Endcap Log10 Energy", "EcalRecHitsTask Endcap Log10 Energy", 90, -5., 4.);
0240     meESRecHitLog10Energy_ =
0241         ibooker.book1D("EcalRecHitsTask Preshower Log10 Energy", "EcalRecHitsTask Preshower Log10 Energy", 90, -5., 4.);
0242
0243     meEERecHitLog10EnergyContr_ = ibooker.bookProfile("EcalRecHits Endcap Log10En vs Hit Contribution",
0244                                                       "EcalRecHits Endcap Log10En vs Hit Contribution",
0245                                                       90,
0246                                                       -5.,
0247                                                       4.,
0248                                                       100,
0249                                                       0.,
0250                                                       1.);
0251     meESRecHitLog10EnergyContr_ = ibooker.bookProfile("EcalRecHits Preshower Log10En vs Hit Contribution",
0252                                                       "EcalRecHits Preshower Log10En vs Hit Contribution",
0253                                                       90,
0254                                                       -5.,
0255                                                       4.,
0256                                                       100,
0257                                                       0.,
0258                                                       1.);
0259
0260     meEERecHitLog10Energy5x5Contr_ = ibooker.bookProfile("EcalRecHits Endcap Log10En5x5 vs Hit Contribution",
0261                                                          "EcalRecHits Endcap Log10En5x5 vs Hit Contribution",
0262                                                          90,
0263                                                          -5.,
0264                                                          4.,
0265                                                          100,
0266                                                          0.,
0267                                                          1.);
0268
0269     histo = "EE+ Occupancy Flag=5 6";
0270     meEERecHitsOccupancyPlusFlag5_6_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0271     histo = "EE- Occupancy Flag=5 6";
0272     meEERecHitsOccupancyMinusFlag5_6_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0273     histo = "EE+ Occupancy Flag=8 9";
0274     meEERecHitsOccupancyPlusFlag8_9_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0275     histo = "EE- Occupancy Flag=8 9";
0276     meEERecHitsOccupancyMinusFlag8_9_ = ibooker.book2D(histo, histo, 100, 0., 100., 100, 0., 100.);
0277
0278     histo = "EcalRecHitsTask Endcap Reco Flags";
0279     meEERecHitFlags_ = ibooker.book1D(histo.c_str(), histo.c_str(), 10, 0., 10.);
0280
0281     histo = "EcalRecHitsTask Endcap RecSimHit Ratio vs SimHit Flag=5 6";
0282     meEERecHitSimHitvsSimHitFlag5_6_ = ibooker.book2D(histo.c_str(), histo.c_str(), 80, 0., 2., 4000, 0., 400.);
0283
0284     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=6";
0285     meEERecHitSimHitFlag6_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0286
0287     histo = "EcalRecHitsTask Endcap RecSimHit Ratio Flag=7";
0288     meEERecHitSimHitFlag7_ = ibooker.book1D(histo.c_str(), histo.c_str(), 80, 0., 2.);
0289   }
0290 }
0291
0292 void EcalRecHitsValidation::analyze(const Event &e, const EventSetup &c) {
0293   // Temporary stuff
0294
0295   LogInfo("EcalRecHitsTask, EventInfo: ") << " Run = " << e.id().run() << " Event = " << e.id().event();
0296
0297   // ADC -> GeV Scale
0298   const EcalADCToGeVConstant *agc = &c.getData(pAgc);
0299   const double barrelADCtoGeV_ = agc->getEBValue();
0300   const double endcapADCtoGeV_ = agc->getEEValue();
0301
0302   Handle<HepMCProduct> MCEvt;
0303   bool skipMC = false;
0304   e.getByToken(HepMCLabel_Token_, MCEvt);
0305   if (!MCEvt.isValid()) {
0306     skipMC = true;
0307   }
0308
0309   edm::Handle<CrossingFrame<PCaloHit>> crossingFrame;
0310
0311   bool skipBarrel = false;
0312   const EBUncalibratedRecHitCollection *EBUncalibRecHit = nullptr;
0313   Handle<EBUncalibratedRecHitCollection> EcalUncalibRecHitEB;
0314   e.getByToken(EBuncalibrechitCollection_Token_, EcalUncalibRecHitEB);
0315   if (EcalUncalibRecHitEB.isValid()) {
0316     EBUncalibRecHit = EcalUncalibRecHitEB.product();
0317   } else {
0318     skipBarrel = true;
0319   }
0320
0321   bool skipEndcap = false;
0322   const EEUncalibratedRecHitCollection *EEUncalibRecHit = nullptr;
0323   Handle<EEUncalibratedRecHitCollection> EcalUncalibRecHitEE;
0324   if (enableEndcaps_) {
0325     e.getByToken(EEuncalibrechitCollection_Token_, EcalUncalibRecHitEE);
0326     if (EcalUncalibRecHitEE.isValid()) {
0327       EEUncalibRecHit = EcalUncalibRecHitEE.product();
0328     } else {
0329       skipEndcap = true;
0330     }
0331   } else
0332     skipEndcap = true;
0333
0334   const EBRecHitCollection *EBRecHit = nullptr;
0335   Handle<EBRecHitCollection> EcalRecHitEB;
0336   e.getByToken(EBrechitCollection_Token_, EcalRecHitEB);
0337   if (EcalRecHitEB.isValid()) {
0338     EBRecHit = EcalRecHitEB.product();
0339   } else {
0340     skipBarrel = true;
0341   }
0342
0343   const EERecHitCollection *EERecHit = nullptr;
0344   Handle<EERecHitCollection> EcalRecHitEE;
0345   if (enableEndcaps_) {
0346     e.getByToken(EErechitCollection_Token_, EcalRecHitEE);
0347     if (EcalRecHitEE.isValid()) {
0348       EERecHit = EcalRecHitEE.product();
0349     } else {
0350       skipEndcap = true;
0351     }
0352   }
0353
0354   bool skipPreshower = false;
0355   const ESRecHitCollection *ESRecHit = nullptr;
0356   Handle<ESRecHitCollection> EcalRecHitES;
0357   if (enableEndcaps_) {
0358     e.getByToken(ESrechitCollection_Token_, EcalRecHitES);
0359     if (EcalRecHitES.isValid()) {
0360       ESRecHit = EcalRecHitES.product();
0361     } else {
0362       skipPreshower = true;
0363     }
0364   } else
0365     skipPreshower = true;
0366
0367   // ----------------------
0368   // gun
0369   double eGun = 0.;
0370   if (!skipMC) {
0371     for (HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
0372          p != MCEvt->GetEvent()->particles_end();
0373          ++p) {
0374       double htheta = (*p)->momentum().theta();
0375       double heta = -99999.;
0376       if (tan(htheta * 0.5) > 0) {
0377         heta = -log(tan(htheta * 0.5));
0378       }
0379       double hphi = (*p)->momentum().phi();
0380       hphi = (hphi >= 0) ? hphi : hphi + 2 * M_PI;
0381       hphi = hphi / M_PI * 180.;
0382
0383       LogDebug("EventInfo") << "EcalRecHitsTask: Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n"
0384                             << "Energy = " << (*p)->momentum().e() << "\n"
0385                             << "Eta = " << heta << "\n"
0386                             << "Phi = " << hphi;
0387
0388       if ((*p)->momentum().e() > eGun)
0389         eGun = (*p)->momentum().e();
0390
0391       if (meGunEnergy_)
0392         meGunEnergy_->Fill((*p)->momentum().e());
0393       if (meGunEta_)
0394         meGunEta_->Fill(heta);
0395       if (meGunPhi_)
0396         meGunPhi_->Fill(hphi);
0397     }
0398   }
0399
0400   // -------------------------------------------------------------------
0401   // BARREL
0402
0403   if (!skipBarrel) {
0404     // 1) loop over simHits
0405     e.getByToken(EBHits_Token_, crossingFrame);
0406     const MixCollection<PCaloHit> barrelHits(crossingFrame.product());
0407
0408     MapType ebSimMap;
0409     MapType ebRecMap;
0410     const int ebcSize = 90;
0411     double ebcontr[ebcSize];
0412     double ebcontr25[ebcSize];
0413     for (int i = 0; i < ebcSize; i++) {
0414       ebcontr[i] = 0.0;
0415       ebcontr25[i] = 0.0;
0416     }
0417     double ebtotal = 0.;
0418
0419     for (auto const &iHit : barrelHits) {
0420       EBDetId ebid = EBDetId(iHit.id());
0421
0422       LogDebug("SimHitInfo, barrel") << "CaloHit " << iHit.getName() << " DetID = " << iHit.id() << "\n"
0423                                      << "Energy = " << iHit.energy() << " Time = " << iHit.time() << "\n"
0424                                      << "EBDetId = " << ebid.ieta() << " " << ebid.iphi();
0425
0426       uint32_t crystid = ebid.rawId();
0427       ebSimMap[crystid] += iHit.energy();
0428     }
0429
0430     // 2) loop over RecHits
0431     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EBUncalibRecHit->begin();
0432          uncalibRecHit != EBUncalibRecHit->end();
0433          ++uncalibRecHit) {
0434       EBDetId EBid = EBDetId(uncalibRecHit->id());
0435
0436       // Find corresponding recHit
0437       EcalRecHitCollection::const_iterator myRecHit = EBRecHit->find(EBid);
0438       if (myRecHit == EBRecHit->end())
0439         continue;
0440       ebRecMap[EBid.rawId()] += myRecHit->energy();
0441
0442       // Fill log10(Energy) stuff...
0443       ebtotal += myRecHit->energy();
0444       if (myRecHit->energy() > 0) {
0445         if (meEBRecHitLog10Energy_)
0446           meEBRecHitLog10Energy_->Fill(log10(myRecHit->energy()));
0447         int log10i = int((log10(myRecHit->energy()) + 5.) * 10.);
0448         if (log10i >= 0 and log10i < ebcSize)
0449           ebcontr[log10i] += myRecHit->energy();
0450       }
0451
0452       // comparison Rec/Sim hit
0453       if (ebSimMap[EBid.rawId()] != 0.) {
0454         double uncEnergy = uncalibRecHit->amplitude() * barrelADCtoGeV_;
0455         if (meEBUnRecHitSimHitRatio_) {
0456           meEBUnRecHitSimHitRatio_->Fill(uncEnergy / ebSimMap[EBid.rawId()]);
0457         }
0458         if (meEBUnRecHitSimHitRatioGt35_ && (myRecHit->energy() > 3.5)) {
0459           meEBUnRecHitSimHitRatioGt35_->Fill(uncEnergy / ebSimMap[EBid.rawId()]);
0460         }
0461       }
0462
0463       if (myRecHit != EBRecHit->end()) {
0464         if (ebSimMap[EBid.rawId()] != 0.) {
0465           if (meEBRecHitSimHitRatio_) {
0466             meEBRecHitSimHitRatio_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0467           }
0468           if (meEBRecHitSimHitRatioGt35_ && (myRecHit->energy() > 3.5)) {
0469             meEBRecHitSimHitRatioGt35_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0470           }
0471           uint16_t sc = 0;
0472
0473           c.getData(pEcsToken);
0474           auto pEcs = c.getHandle(pEcsToken);
0475           const EcalChannelStatus *ecs = nullptr;
0476           if (pEcs.isValid())
0477             ecs = &c.getData(pEcsToken);
0478           ;
0479           if (ecs != nullptr) {
0480             EcalChannelStatusMap::const_iterator csmi = ecs->find(EBid.rawId());
0481             EcalChannelStatusCode csc = 0;
0482             if (csmi != ecs->end())
0483               csc = *csmi;
0484             sc = csc.getStatusCode();
0485           }
0486
0487           if (meEBRecHitSimHitRatio1011_ != nullptr && (sc == 10 || sc == 11)) {
0488             meEBRecHitSimHitRatio1011_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0489           }
0490           if (meEBRecHitSimHitRatio12_ != nullptr && sc == 12) {
0491             meEBRecHitSimHitRatio12_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0492           }
0493
0494           auto pttMap = c.getHandle(pttMapToken);
0495           const EcalTrigTowerConstituentsMap *ttMap = nullptr;
0496           if (pttMap.isValid())
0497             ttMap = &c.getData(pttMapToken);
0498           double ttSimEnergy = 0;
0499           if (ttMap != nullptr) {
0500             EcalTrigTowerDetId ttDetId = EBid.tower();
0501             std::vector<DetId> vid = ttMap->constituentsOf(ttDetId);
0502             for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); dit++) {
0503               EBDetId ttEBid = EBDetId(*dit);
0504               ttSimEnergy += ebSimMap[ttEBid.rawId()];
0505             }
0506             if (!vid.empty())
0507               ttSimEnergy = ttSimEnergy / vid.size();
0508           }
0509           if (meEBRecHitSimHitRatio13_ != nullptr && sc == 13 && ttSimEnergy != 0)
0510             meEBRecHitSimHitRatio13_->Fill(myRecHit->energy() / ttSimEnergy);
0511
0512           int flag = myRecHit->recoFlag();
0513           if (meEBRecHitFlags_ != nullptr)
0514             meEBRecHitFlags_->Fill(flag);
0515           if (meEBRecHitSimHitvsSimHitFlag5_6_ &&
0516               (flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered))
0517             meEBRecHitSimHitvsSimHitFlag5_6_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()], ebSimMap[EBid.rawId()]);
0518           if (meEBRecHitSimHitFlag6_ && (flag == EcalRecHit::kLeadingEdgeRecovered))
0519             meEBRecHitSimHitFlag6_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0520           if (meEBRecHitSimHitFlag7_ && (flag == EcalRecHit::kNeighboursRecovered))
0521             meEBRecHitSimHitFlag6_->Fill(myRecHit->energy() / ebSimMap[EBid.rawId()]);
0522           if (meEB5x5RecHitSimHitvsSimHitFlag8_ && (flag == EcalRecHit::kTowerRecovered) && ttSimEnergy != 0)
0523             meEB5x5RecHitSimHitvsSimHitFlag8_->Fill(myRecHit->energy() / ttSimEnergy, ttSimEnergy);
0524
0525           if (meEBRecHitsOccupancyFlag5_6_ &&
0526               ((flag == EcalRecHit::kSaturated) || (flag == EcalRecHit::kLeadingEdgeRecovered)))
0527             meEBRecHitsOccupancyFlag5_6_->Fill(EBid.ieta(), EBid.iphi());
0528           if (meEBRecHitsOccupancyFlag8_9_ && ((flag == EcalRecHit::kTowerRecovered) || (flag == EcalRecHit::kDead)))
0529             meEBRecHitsOccupancyFlag8_9_->Fill(EBid.ieta(), EBid.iphi());
0530         }
0531       } else
0532         continue;
0533     }  // loop over the UncalibratedRecHitCollection
0534
0535     // RecHits matrix
0536     uint32_t ebcenterid = getUnitWithMaxEnergy(ebRecMap);
0537     EBDetId myEBid(ebcenterid);
0538     int bx = myEBid.ietaAbs();
0539     int by = myEBid.iphi();
0540     int bz = myEBid.zside();
0541     findBarrelMatrix(5, 5, bx, by, bz, ebRecMap);
0542     double e5x5rec = 0.;
0543     double e5x5sim = 0.;
0544     for (unsigned int i = 0; i < crystalMatrix.size(); i++) {
0545       e5x5rec += ebRecMap[crystalMatrix[i]];
0546       e5x5sim += ebSimMap[crystalMatrix[i]];
0547       if (ebRecMap[crystalMatrix[i]] > 0) {
0548         int log10i25 = int((log10(ebRecMap[crystalMatrix[i]]) + 5.) * 10.);
0549         if (log10i25 >= 0 && log10i25 < ebcSize)
0550           ebcontr25[log10i25] += ebRecMap[crystalMatrix[i]];
0551       }
0552     }
0553
0554     if (meEBe5x5_)
0555       meEBe5x5_->Fill(e5x5rec);
0556     if (e5x5sim > 0. && meEBe5x5OverSimHits_)
0557       meEBe5x5OverSimHits_->Fill(e5x5rec / e5x5sim);
0558     if (eGun > 0. && meEBe5x5OverGun_)
0559       meEBe5x5OverGun_->Fill(e5x5rec / eGun);
0560
0561     if (meEBRecHitLog10EnergyContr_ && ebtotal != 0) {
0562       for (int i = 0; i < ebcSize; i++) {
0563         meEBRecHitLog10EnergyContr_->Fill(-5. + (float(i) + 0.5) / 10., ebcontr[i] / ebtotal);
0564       }
0565     }
0566
0567     if (meEBRecHitLog10Energy5x5Contr_ && e5x5rec != 0) {
0568       for (int i = 0; i < ebcSize; i++) {
0569         meEBRecHitLog10Energy5x5Contr_->Fill(-5. + (float(i) + 0.5) / 10., ebcontr25[i] / e5x5rec);
0570       }
0571     }
0572   }
0573
0574   // -------------------------------------------------------------------
0575   // ENDCAP
0576
0577   if (!skipEndcap) {
0578     // 1) loop over simHits
0579     e.getByToken(EEHits_Token_, crossingFrame);
0580     const MixCollection<PCaloHit> endcapHits(crossingFrame.product());
0581
0582     MapType eeSimMap;
0583     MapType eeRecMap;
0584     const int eecSize = 90;
0585     double eecontr[eecSize];
0586     double eecontr25[eecSize];
0587     for (int i = 0; i < eecSize; i++) {
0588       eecontr[i] = 0.0;
0589       eecontr25[i] = 0.0;
0590     }
0591     double eetotal = 0.;
0592
0593     for (auto const &iHit : endcapHits) {
0594       EEDetId eeid(iHit.id());
0595
0596       LogDebug("Endcap, HitInfo") << " CaloHit " << iHit.getName() << " DetID = " << iHit.id() << "\n"
0597                                   << "Energy = " << iHit.energy() << " Time = " << iHit.time() << "\n"
0598                                   << "EEDetId side " << eeid.zside() << " = " << eeid.ix() << " " << eeid.iy();
0599
0600       uint32_t crystid = eeid.rawId();
0601       eeSimMap[crystid] += iHit.energy();
0602     }
0603
0604     // 2) loop over RecHits
0605     for (EcalUncalibratedRecHitCollection::const_iterator uncalibRecHit = EEUncalibRecHit->begin();
0606          uncalibRecHit != EEUncalibRecHit->end();
0607          ++uncalibRecHit) {
0608       EEDetId EEid = EEDetId(uncalibRecHit->id());
0609
0610       // Find corresponding recHit
0611       EcalRecHitCollection::const_iterator myRecHit = EERecHit->find(EEid);
0612       if (myRecHit == EERecHit->end())
0613         continue;
0614       eeRecMap[EEid.rawId()] += myRecHit->energy();
0615
0616       // Fill log10(Energy) stuff...
0617       eetotal += myRecHit->energy();
0618       if (myRecHit->energy() > 0) {
0619         if (meEERecHitLog10Energy_)
0620           meEERecHitLog10Energy_->Fill(log10(myRecHit->energy()));
0621         int log10i = int((log10(myRecHit->energy()) + 5.) * 10.);
0622         if (log10i >= 0 and log10i < eecSize)
0623           eecontr[log10i] += myRecHit->energy();
0624       }
0625
0626       // comparison Rec/Sim hit
0627       if (eeSimMap[EEid.rawId()] != 0.) {
0628         double uncEnergy = uncalibRecHit->amplitude() * endcapADCtoGeV_;
0629         if (meEEUnRecHitSimHitRatio_) {
0630           meEEUnRecHitSimHitRatio_->Fill(uncEnergy / eeSimMap[EEid.rawId()]);
0631         }
0632         if (meEEUnRecHitSimHitRatioGt35_ && (myRecHit->energy() > 3.5)) {
0633           meEEUnRecHitSimHitRatioGt35_->Fill(uncEnergy / eeSimMap[EEid.rawId()]);
0634         }
0635       }
0636
0637       if (myRecHit != EERecHit->end()) {
0638         if (eeSimMap[EEid.rawId()] != 0.) {
0639           if (meEERecHitSimHitRatio_) {
0640             meEERecHitSimHitRatio_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0641           }
0642           if (meEERecHitSimHitRatioGt35_ && (myRecHit->energy() > 3.5)) {
0643             meEERecHitSimHitRatioGt35_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0644           }
0645
0646           auto pEcs = c.getHandle(pEcsToken);
0647           const EcalChannelStatus *ecs = nullptr;
0648           if (pEcs.isValid())
0649             ecs = &c.getData(pEcsToken);
0650           if (ecs != nullptr) {
0651             EcalChannelStatusMap::const_iterator csmi = ecs->find(EEid.rawId());
0652             EcalChannelStatusCode csc = 0;
0653             if (csmi != ecs->end())
0654               csc = *csmi;
0655             uint16_t sc = csc.getStatusCode();
0656             if (meEERecHitSimHitRatio1011_ != nullptr && (sc == 10 || sc == 11)) {
0657               meEERecHitSimHitRatio1011_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0658             }
0659             if (meEERecHitSimHitRatio12_ != nullptr && sc == 12) {
0660               meEERecHitSimHitRatio12_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0661             }
0662             if (meEERecHitSimHitRatio13_ != nullptr && sc == 13) {
0663               meEERecHitSimHitRatio13_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0664             }
0665           }
0666
0667           int flag = myRecHit->recoFlag();
0668           if (meEERecHitFlags_ != nullptr)
0669             meEERecHitFlags_->Fill(flag);
0670           if (meEERecHitSimHitvsSimHitFlag5_6_ &&
0671               (flag == EcalRecHit::kSaturated || flag == EcalRecHit::kLeadingEdgeRecovered))
0672             meEERecHitSimHitvsSimHitFlag5_6_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()], eeSimMap[EEid.rawId()]);
0673           if (meEERecHitSimHitFlag6_ && (flag == EcalRecHit::kLeadingEdgeRecovered))
0674             meEERecHitSimHitFlag6_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0675           if (meEERecHitSimHitFlag7_ && (flag == EcalRecHit::kNeighboursRecovered))
0676             meEERecHitSimHitFlag6_->Fill(myRecHit->energy() / eeSimMap[EEid.rawId()]);
0677
0678           if (EEid.zside() > 0) {
0679             if (meEERecHitsOccupancyPlusFlag5_6_ &&
0680                 ((flag == EcalRecHit::kSaturated) || (flag == EcalRecHit::kLeadingEdgeRecovered)))
0681               meEERecHitsOccupancyPlusFlag5_6_->Fill(EEid.ix(), EEid.iy());
0682             if (meEERecHitsOccupancyPlusFlag8_9_ &&
0683                 ((flag == EcalRecHit::kTowerRecovered) || (flag == EcalRecHit::kDead)))
0684               meEERecHitsOccupancyPlusFlag8_9_->Fill(EEid.ix(), EEid.iy());
0685           }
0686           if (EEid.zside() < 0) {
0687             if (meEERecHitsOccupancyMinusFlag5_6_ &&
0688                 ((flag == EcalRecHit::kSaturated) || (flag == EcalRecHit::kLeadingEdgeRecovered)))
0689               meEERecHitsOccupancyMinusFlag5_6_->Fill(EEid.ix(), EEid.iy());
0690             if (meEERecHitsOccupancyMinusFlag8_9_ &&
0691                 ((flag == EcalRecHit::kTowerRecovered) || (flag == EcalRecHit::kDead)))
0692               meEERecHitsOccupancyMinusFlag8_9_->Fill(EEid.ix(), EEid.iy());
0693           }
0694         }
0695       } else
0696         continue;
0697     }  // loop over the UncalibratedechitCollection
0698
0699     // RecHits matrix
0700     uint32_t eecenterid = getUnitWithMaxEnergy(eeRecMap);
0701     EEDetId myEEid(eecenterid);
0702     int bx = myEEid.ix();
0703     int by = myEEid.iy();
0704     int bz = myEEid.zside();
0705     findEndcapMatrix(5, 5, bx, by, bz, eeRecMap);
0706     double e5x5rec = 0.;
0707     double e5x5sim = 0.;
0708     for (unsigned int i = 0; i < crystalMatrix.size(); i++) {
0709       e5x5rec += eeRecMap[crystalMatrix[i]];
0710       e5x5sim += eeSimMap[crystalMatrix[i]];
0711       if (eeRecMap[crystalMatrix[i]] > 0) {
0712         int log10i25 = int((log10(eeRecMap[crystalMatrix[i]]) + 5.) * 10.);
0713         if (log10i25 >= 0 && log10i25 < eecSize)
0714           eecontr25[log10i25] += eeRecMap[crystalMatrix[i]];
0715       }
0716     }
0717
0718     if (meEEe5x5_)
0719       meEEe5x5_->Fill(e5x5rec);
0720     if (e5x5sim > 0. && meEEe5x5OverSimHits_)
0721       meEEe5x5OverSimHits_->Fill(e5x5rec / e5x5sim);
0722     if (eGun > 0. && meEEe5x5OverGun_)
0723       meEEe5x5OverGun_->Fill(e5x5rec / eGun);
0724
0725     if (meEERecHitLog10EnergyContr_ && eetotal != 0) {
0726       for (int i = 0; i < eecSize; i++) {
0727         meEERecHitLog10EnergyContr_->Fill(-5. + (float(i) + 0.5) / 10., eecontr[i] / eetotal);
0728       }
0729     }
0730
0731     if (meEERecHitLog10Energy5x5Contr_ && e5x5rec != 0) {
0732       for (int i = 0; i < eecSize; i++) {
0733         meEERecHitLog10Energy5x5Contr_->Fill(-5. + (float(i) + 0.5) / 10., eecontr25[i] / e5x5rec);
0734       }
0735     }
0736   }
0737
0738   // -------------------------------------------------------------------
0739   // PRESHOWER
0740
0741   if (!skipPreshower) {
0742     // 1) loop over simHits
0743     e.getByToken(ESHits_Token_, crossingFrame);
0744     const MixCollection<PCaloHit> preshowerHits(crossingFrame.product());
0745
0746     MapType esSimMap;
0747     const int escSize = 90;
0748     double escontr[escSize];
0749     for (int i = 0; i < escSize; i++) {
0750       escontr[i] = 0.0;
0751     }
0752     double estotal = 0.;
0753
0754     for (auto const &iHit : preshowerHits) {
0755       ESDetId esid(iHit.id());
0756       LogDebug("Preshower, HitInfo") << " CaloHit " << iHit.getName() << " DetID = " << iHit.id() << "\n"
0757                                      << "Energy = " << iHit.energy() << " Time = " << iHit.time() << "\n"
0758                                      << "ESDetId strip " << esid.strip() << " = " << esid.six() << " " << esid.siy();
0759
0760       uint32_t crystid = esid.rawId();
0761       esSimMap[crystid] += iHit.energy();
0762     }
0763
0764     // 2) loop over RecHits
0765     for (EcalRecHitCollection::const_iterator recHit = ESRecHit->begin(); recHit != ESRecHit->end(); ++recHit) {
0766       ESDetId ESid = ESDetId(recHit->id());
0767       if (esSimMap[ESid.rawId()] != 0.) {
0768         // Fill log10(Energy) stuff...
0769         estotal += recHit->energy();
0770         if (recHit->energy() > 0) {
0771           if (meESRecHitLog10Energy_)
0772             meESRecHitLog10Energy_->Fill(log10(recHit->energy()));
0773           int log10i = int((log10(recHit->energy()) + 5.) * 10.);
0774           if (log10i >= 0 and log10i < escSize)
0775             escontr[log10i] += recHit->energy();
0776         }
0777
0778         if (meESRecHitSimHitRatio_) {
0779           meESRecHitSimHitRatio_->Fill(recHit->energy() / esSimMap[ESid.rawId()]);
0780         }
0781       } else
0782         continue;
0783     }  // loop over the RechitCollection
0784
0785     if (meESRecHitLog10EnergyContr_ && estotal != 0) {
0786       for (int i = 0; i < escSize; i++) {
0787         meESRecHitLog10EnergyContr_->Fill(-5. + (float(i) + 0.5) / 10., escontr[i] / estotal);
0788       }
0789     }
0790   }
0791 }
0792
0793 uint32_t EcalRecHitsValidation::getUnitWithMaxEnergy(MapType &themap) {
0794   // look for max
0795   uint32_t unitWithMaxEnergy = 0;
0796   float maxEnergy = 0.;
0797
0798   MapType::iterator iter;
0799   for (iter = themap.begin(); iter != themap.end(); iter++) {
0800     if (maxEnergy < (*iter).second) {
0801       maxEnergy = (*iter).second;
0802       unitWithMaxEnergy = (*iter).first;
0803     }
0804   }
0805
0806   return unitWithMaxEnergy;
0807 }
0808
0809 void EcalRecHitsValidation::findBarrelMatrix(
0810     int nCellInEta, int nCellInPhi, int CentralEta, int CentralPhi, int CentralZ, MapType &themap) {
0811   int goBackInEta = nCellInEta / 2;
0812   int goBackInPhi = nCellInPhi / 2;
0813   int matrixSize = nCellInEta * nCellInPhi;
0814   crystalMatrix.clear();
0815   crystalMatrix.resize(matrixSize);
0816
0817   int startEta = CentralZ * CentralEta - goBackInEta;
0818   int startPhi = CentralPhi - goBackInPhi;
0819
0820   int i = 0;
0821   for (int ieta = startEta; ieta < startEta + nCellInEta; ieta++) {
0822     for (int iphi = startPhi; iphi < startPhi + nCellInPhi; iphi++) {
0823       uint32_t index;
0824       if (abs(ieta) > 85 || abs(ieta) < 1) {
0825         continue;
0826       }
0827       if (iphi < 1) {
0828         index = EBDetId(ieta, iphi + 360).rawId();
0829       } else if (iphi > 360) {
0830         index = EBDetId(ieta, iphi - 360).rawId();
0831       } else {
0832         index = EBDetId(ieta, iphi).rawId();
0833       }
0834       crystalMatrix[i++] = index;
0835     }
0836   }
0837 }
0838
0839 void EcalRecHitsValidation::findEndcapMatrix(
0840     int nCellInX, int nCellInY, int CentralX, int CentralY, int CentralZ, MapType &themap) {
0841   int goBackInX = nCellInX / 2;
0842   int goBackInY = nCellInY / 2;
0843   crystalMatrix.clear();
0844
0845   int startX = CentralX - goBackInX;
0846   int startY = CentralY - goBackInY;
0847
0848   for (int ix = startX; ix < startX + nCellInX; ix++) {
0849     for (int iy = startY; iy < startY + nCellInY; iy++) {
0850       uint32_t index;
0851
0852       if (EEDetId::validDetId(ix, iy, CentralZ)) {
0853         index = EEDetId(ix, iy, CentralZ).rawId();
0854       } else {
0855         continue;
0856       }
0857       crystalMatrix.push_back(index);
0858     }
0859   }
0860 }