Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:41

0001 /// -*- C++ -*-
0002 //
0003 // Package:    HGCMissingRecHit
0004 // Class:      HGCMissingRecHit
0005 //
0006 /**\class HGCMissingRecHit HGCMissingRecHit.cc Validation/HGCalValidation/test/HGCMissingRecHit.cc
0007 
0008  Description: [one line class summary]
0009 
0010  Implementation:
0011     [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  "Sunanda Banerjee"
0015 //         Created:  Tue September 20 17:55:26 CDT 2022
0016 // $Id$
0017 //
0018 //
0019 
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021 
0022 #include "DataFormats/Common/interface/Handle.h"
0023 #include "DataFormats/DetId/interface/DetId.h"
0024 #include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
0025 #include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0027 #include "DataFormats/HGCDigi/interface/HGCDigiCollections.h"
0028 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0029 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0030 
0031 #include "FWCore/Framework/interface/Frameworkfwd.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0035 #include "FWCore/Framework/interface/ESHandle.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "FWCore/Utilities/interface/Exception.h"
0039 #include "FWCore/Utilities/interface/EDGetToken.h"
0040 #include "FWCore/Utilities/interface/transform.h"
0041 
0042 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0043 #include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h"
0044 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0045 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0046 
0047 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0048 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0049 
0050 #include <cmath>
0051 #include <iostream>
0052 #include <memory>
0053 #include <sstream>
0054 #include <string>
0055 #include <vector>
0056 
0057 #include <TH1D.h>
0058 
0059 class HGCMissingRecHit : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0060 public:
0061   explicit HGCMissingRecHit(const edm::ParameterSet &);
0062   ~HGCMissingRecHit() override = default;
0063   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0064 
0065 private:
0066   typedef std::tuple<float, GlobalPoint> HGCHitTuple;
0067 
0068   void beginJob() override;
0069   void endJob() override {}
0070   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0071   void analyze(edm::Event const &, edm::EventSetup const &) override;
0072   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0073   virtual void beginLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0074   virtual void endLuminosityBlock(edm::LuminosityBlock const &, edm::EventSetup const &) {}
0075   void analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0076                           int idet,
0077                           std::map<unsigned int, HGCHitTuple> &);
0078   template <class T1>
0079   void analyzeHGCalDigi(T1 const &theHits, int idet, std::map<unsigned int, HGCHitTuple> const &hitRefs);
0080   template <class T1>
0081   void analyzeHGCalRecHit(T1 const &theHits, int idet, std::map<unsigned int, HGCHitTuple> const &hitRefs);
0082 
0083 private:
0084   //HGC Geometry
0085   const std::vector<std::string> geometrySource_, detectors_;
0086   const std::vector<int> ietaExcludeBH_;
0087   const std::vector<edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord>> tok_hgcal_;
0088   const std::vector<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>> tok_hgcalg_;
0089   std::vector<const HGCalDDDConstants *> hgcCons_;
0090   std::vector<const HGCalGeometry *> hgcGeometry_;
0091 
0092   const edm::InputTag eeSimHitSource, fhSimHitSource, bhSimHitSource;
0093   const edm::EDGetTokenT<std::vector<PCaloHit>> eeSimHitToken_;
0094   const edm::EDGetTokenT<std::vector<PCaloHit>> fhSimHitToken_;
0095   const edm::EDGetTokenT<std::vector<PCaloHit>> bhSimHitToken_;
0096   const edm::InputTag eeDigiSource, fhDigiSource, bhDigiSource;
0097   const edm::EDGetTokenT<HGCalDigiCollection> eeDigiToken_;
0098   const edm::EDGetTokenT<HGCalDigiCollection> fhDigiToken_;
0099   const edm::EDGetTokenT<HGCalDigiCollection> bhDigiToken_;
0100   const edm::InputTag eeRecHitSource, fhRecHitSource, bhRecHitSource;
0101   const edm::EDGetTokenT<HGCeeRecHitCollection> eeRecHitToken_;
0102   const edm::EDGetTokenT<HGChefRecHitCollection> fhRecHitToken_;
0103   const edm::EDGetTokenT<HGChebRecHitCollection> bhRecHitToken_;
0104 
0105   std::vector<TH1D *> goodHitsDE_, missedHitsDE_, goodHitsDT_, missedHitsDT_;
0106   std::vector<TH1D *> goodHitsRE_, missedHitsRE_, goodHitsRT_, missedHitsRT_;
0107 };
0108 
0109 HGCMissingRecHit::HGCMissingRecHit(const edm::ParameterSet &cfg)
0110     : geometrySource_(cfg.getParameter<std::vector<std::string>>("geometrySource")),
0111       detectors_(cfg.getParameter<std::vector<std::string>>("detectors")),
0112       ietaExcludeBH_(cfg.getParameter<std::vector<int>>("ietaExcludeBH")),
0113       tok_hgcal_{
0114           edm::vector_transform(geometrySource_,
0115                                 [this](const std::string &name) {
0116                                   return esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0117                                       edm::ESInputTag{"", name});
0118                                 })},
0119       tok_hgcalg_{edm::vector_transform(
0120           geometrySource_,
0121           [this](const std::string &name) {
0122             return esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(edm::ESInputTag{"", name});
0123           })},
0124       eeSimHitSource(cfg.getParameter<edm::InputTag>("eeSimHitSource")),
0125       fhSimHitSource(cfg.getParameter<edm::InputTag>("fhSimHitSource")),
0126       bhSimHitSource(cfg.getParameter<edm::InputTag>("bhSimHitSource")),
0127       eeSimHitToken_(consumes<std::vector<PCaloHit>>(eeSimHitSource)),
0128       fhSimHitToken_(consumes<std::vector<PCaloHit>>(fhSimHitSource)),
0129       bhSimHitToken_(consumes<std::vector<PCaloHit>>(bhSimHitSource)),
0130       eeDigiSource(cfg.getParameter<edm::InputTag>("eeDigiSource")),
0131       fhDigiSource(cfg.getParameter<edm::InputTag>("fhDigiSource")),
0132       bhDigiSource(cfg.getParameter<edm::InputTag>("bhDigiSource")),
0133       eeDigiToken_(consumes<HGCalDigiCollection>(eeDigiSource)),
0134       fhDigiToken_(consumes<HGCalDigiCollection>(fhDigiSource)),
0135       bhDigiToken_(consumes<HGCalDigiCollection>(bhDigiSource)),
0136       eeRecHitSource(cfg.getParameter<edm::InputTag>("eeRecHitSource")),
0137       fhRecHitSource(cfg.getParameter<edm::InputTag>("fhRecHitSource")),
0138       bhRecHitSource(cfg.getParameter<edm::InputTag>("bhRecHitSource")),
0139       eeRecHitToken_(consumes<HGCeeRecHitCollection>(eeRecHitSource)),
0140       fhRecHitToken_(consumes<HGChefRecHitCollection>(fhRecHitSource)),
0141       bhRecHitToken_(consumes<HGChebRecHitCollection>(bhRecHitSource)) {
0142   usesResource(TFileService::kSharedResource);
0143 
0144   edm::LogVerbatim("HGCalValid") << "Use " << geometrySource_.size() << " Geometry sources";
0145   for (unsigned int k = 0; k < geometrySource_.size(); k++)
0146     edm::LogVerbatim("HGCalValid") << "  " << detectors_[k] << ":" << geometrySource_[k];
0147   edm::LogVerbatim("HGCalValid") << "SimHit labels: " << eeSimHitSource << "  " << fhSimHitSource << "  "
0148                                  << bhSimHitSource;
0149   edm::LogVerbatim("HGCalValid") << "Digi labels: " << eeDigiSource << "  " << fhDigiSource << "  " << bhDigiSource;
0150   edm::LogVerbatim("HGCalValid") << "RecHit labels: " << eeRecHitSource << "  " << fhRecHitSource << "  "
0151                                  << bhRecHitSource;
0152   edm::LogVerbatim("HGCalValid") << "Exclude the following " << ietaExcludeBH_.size() << " ieta values from BH plots";
0153   for (unsigned int k = 0; k < ietaExcludeBH_.size(); ++k)
0154     edm::LogVerbatim("HGCalValid") << " [" << k << "] " << ietaExcludeBH_[k];
0155 }
0156 
0157 void HGCMissingRecHit::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0158   std::vector<std::string> sources = {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"};
0159   std::vector<std::string> names = {"EE", "HE Silicon", "HE Scintillator"};
0160   std::vector<int> etas;
0161   edm::ParameterSetDescription desc;
0162   desc.add<std::vector<std::string>>("geometrySource", sources);
0163   desc.add<std::vector<std::string>>("detectors", names);
0164   desc.add<edm::InputTag>("eeSimHitSource", edm::InputTag("g4SimHits", "HGCHitsEE"));
0165   desc.add<edm::InputTag>("fhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEfront"));
0166   desc.add<edm::InputTag>("bhSimHitSource", edm::InputTag("g4SimHits", "HGCHitsHEback"));
0167   desc.add<edm::InputTag>("eeDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "EE"));
0168   desc.add<edm::InputTag>("fhDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "HEfront"));
0169   desc.add<edm::InputTag>("bhDigiSource", edm::InputTag("simHGCalUnsuppressedDigis", "HEback"));
0170   desc.add<edm::InputTag>("eeRecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0171   desc.add<edm::InputTag>("fhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEFRecHits"));
0172   desc.add<edm::InputTag>("bhRecHitSource", edm::InputTag("HGCalRecHit", "HGCHEBRecHits"));
0173   desc.add<std::vector<int>>("ietaExcludeBH", etas);
0174   descriptions.add("hgcMissingRecHit", desc);
0175 }
0176 
0177 void HGCMissingRecHit::beginJob() {
0178   //initiating fileservice
0179   edm::Service<TFileService> fs;
0180   for (unsigned int k = 0; k < detectors_.size(); ++k) {
0181     char name[50], title[100];
0182     sprintf(name, "GoodDE%s", geometrySource_[k].c_str());
0183     sprintf(title, "SimHit energy present among Digis in %s", detectors_[k].c_str());
0184     goodHitsDE_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 0.01));
0185     goodHitsDE_.back()->Sumw2();
0186     sprintf(name, "MissDE%s", geometrySource_[k].c_str());
0187     sprintf(title, "SimHit energy absent among Digis in %s", detectors_[k].c_str());
0188     missedHitsDE_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 0.01));
0189     missedHitsDE_.back()->Sumw2();
0190     sprintf(name, "GoodDT%s", geometrySource_[k].c_str());
0191     sprintf(title, "|#eta| of SimHits present among Digis in %s", detectors_[k].c_str());
0192     goodHitsDT_.emplace_back(fs->make<TH1D>(name, title, 320, 2.7, 3.1));
0193     goodHitsDT_.back()->Sumw2();
0194     sprintf(name, "MissDT%s", geometrySource_[k].c_str());
0195     sprintf(title, "|#eta| of SimHits absent among Digis in %s", detectors_[k].c_str());
0196     missedHitsDT_.emplace_back(fs->make<TH1D>(name, title, 320, 2.7, 3.1));
0197     missedHitsDT_.back()->Sumw2();
0198     sprintf(name, "GoodRE%s", geometrySource_[k].c_str());
0199     sprintf(title, "SimHit energy present among RecHits in %s", detectors_[k].c_str());
0200     goodHitsRE_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 0.01));
0201     goodHitsRE_.back()->Sumw2();
0202     sprintf(name, "MissRE%s", geometrySource_[k].c_str());
0203     sprintf(title, "SimHit energy absent among RecHits in %s", detectors_[k].c_str());
0204     missedHitsRE_.emplace_back(fs->make<TH1D>(name, title, 1000, 0, 0.01));
0205     missedHitsRE_.back()->Sumw2();
0206     sprintf(name, "GoodRT%s", geometrySource_[k].c_str());
0207     sprintf(title, "|#eta| of SimHits present among RecHits in %s", detectors_[k].c_str());
0208     goodHitsRT_.emplace_back(fs->make<TH1D>(name, title, 320, 2.7, 3.1));
0209     goodHitsRT_.back()->Sumw2();
0210     sprintf(name, "MissRT%s", geometrySource_[k].c_str());
0211     sprintf(title, "|#eta| of SimHits absent among RecHits in %s", detectors_[k].c_str());
0212     missedHitsRT_.emplace_back(fs->make<TH1D>(name, title, 320, 2.7, 3.1));
0213     missedHitsRT_.back()->Sumw2();
0214   }
0215 }
0216 
0217 void HGCMissingRecHit::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
0218   //initiating hgc Geometry
0219   for (size_t i = 0; i < geometrySource_.size(); i++) {
0220     edm::LogVerbatim("HGCalValid") << "Tries to initialize HGCalGeometry and HGCalDDDConstants for " << i;
0221     const edm::ESHandle<HGCalDDDConstants> &hgcCons = iSetup.getHandle(tok_hgcal_[i]);
0222     if (hgcCons.isValid()) {
0223       hgcCons_.push_back(hgcCons.product());
0224     } else {
0225       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalDDDConstants for " << geometrySource_[i] << std::endl;
0226     }
0227     const edm::ESHandle<HGCalGeometry> &hgcGeom = iSetup.getHandle(tok_hgcalg_[i]);
0228     if (hgcGeom.isValid()) {
0229       hgcGeometry_.push_back(hgcGeom.product());
0230     } else {
0231       edm::LogWarning("HGCalValid") << "Cannot initiate HGCalGeometry for " << geometrySource_[i] << std::endl;
0232     }
0233   }
0234 }
0235 
0236 void HGCMissingRecHit::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0237   std::map<unsigned int, HGCHitTuple> eeHitRefs, fhHitRefs, bhHitRefs;
0238 
0239   //Accesing ee simhits
0240   const edm::Handle<std::vector<PCaloHit>> &eeSimHits = iEvent.getHandle(eeSimHitToken_);
0241   if (eeSimHits.isValid()) {
0242     analyzeHGCalSimHit(eeSimHits, 0, eeHitRefs);
0243     for (std::map<unsigned int, HGCHitTuple>::iterator itr = eeHitRefs.begin(); itr != eeHitRefs.end(); ++itr) {
0244       int idx = std::distance(eeHitRefs.begin(), itr);
0245       edm::LogVerbatim("HGCalValid") << "EEHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0246                                      << std::get<0>(itr->second) << "; Position = " << std::get<1>(itr->second) << ")";
0247     }
0248   } else {
0249     edm::LogWarning("HGCalValid") << "No EE SimHit Found " << std::endl;
0250   }
0251 
0252   //Accesing fh simhits
0253   const edm::Handle<std::vector<PCaloHit>> &fhSimHits = iEvent.getHandle(fhSimHitToken_);
0254   if (fhSimHits.isValid()) {
0255     analyzeHGCalSimHit(fhSimHits, 1, fhHitRefs);
0256     for (std::map<unsigned int, HGCHitTuple>::iterator itr = fhHitRefs.begin(); itr != fhHitRefs.end(); ++itr) {
0257       int idx = std::distance(fhHitRefs.begin(), itr);
0258       edm::LogVerbatim("HGCalValid") << "FHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0259                                      << std::get<0>(itr->second) << "; Position = " << std::get<1>(itr->second) << ")";
0260     }
0261   } else {
0262     edm::LogWarning("HGCalValid") << "No FH SimHit Found " << std::endl;
0263   }
0264 
0265   //Accessing bh simhits
0266   const edm::Handle<std::vector<PCaloHit>> &bhSimHits = iEvent.getHandle(bhSimHitToken_);
0267   if (bhSimHits.isValid()) {
0268     analyzeHGCalSimHit(bhSimHits, 2, bhHitRefs);
0269     for (std::map<unsigned int, HGCHitTuple>::iterator itr = bhHitRefs.begin(); itr != bhHitRefs.end(); ++itr) {
0270       int idx = std::distance(bhHitRefs.begin(), itr);
0271       edm::LogVerbatim("HGCalValid") << "BHHit[" << idx << "] " << std::hex << itr->first << std::dec << "; Energy "
0272                                      << std::get<0>(itr->second) << "; Position = (" << std::get<1>(itr->second) << ")";
0273     }
0274   } else {
0275     edm::LogWarning("HGCalValid") << "No BH SimHit Found " << std::endl;
0276   }
0277 
0278   //accessing EE Digi information
0279   const edm::Handle<HGCalDigiCollection> &eeDigi = iEvent.getHandle(eeDigiToken_);
0280   if (eeDigi.isValid()) {
0281     const HGCalDigiCollection *theHits = (eeDigi.product());
0282     analyzeHGCalDigi(theHits, 0, eeHitRefs);
0283   } else {
0284     edm::LogWarning("HGCalValid") << "No EE Digi Found " << std::endl;
0285   }
0286 
0287   //accessing FH Digi information
0288   const edm::Handle<HGCalDigiCollection> &fhDigi = iEvent.getHandle(fhDigiToken_);
0289   if (fhDigi.isValid()) {
0290     const HGCalDigiCollection *theHits = (fhDigi.product());
0291     analyzeHGCalDigi(theHits, 1, fhHitRefs);
0292   } else {
0293     edm::LogWarning("HGCalValid") << "No FH Digi Found " << std::endl;
0294   }
0295 
0296   //accessing BH Digi information
0297   const edm::Handle<HGCalDigiCollection> &bhDigi = iEvent.getHandle(bhDigiToken_);
0298   if (bhDigi.isValid()) {
0299     const HGCalDigiCollection *theHits = (bhDigi.product());
0300     analyzeHGCalDigi(theHits, 2, bhHitRefs);
0301   } else {
0302     edm::LogWarning("HGCalValid") << "No BH Digi Found " << std::endl;
0303   }
0304 
0305   //accessing EE Rechit information
0306   const edm::Handle<HGCeeRecHitCollection> &eeRecHit = iEvent.getHandle(eeRecHitToken_);
0307   if (eeRecHit.isValid()) {
0308     const HGCeeRecHitCollection *theHits = (eeRecHit.product());
0309     analyzeHGCalRecHit(theHits, 0, eeHitRefs);
0310   } else {
0311     edm::LogWarning("HGCalValid") << "No EE RecHit Found " << std::endl;
0312   }
0313 
0314   //accessing FH Rechit information
0315   const edm::Handle<HGChefRecHitCollection> &fhRecHit = iEvent.getHandle(fhRecHitToken_);
0316   if (fhRecHit.isValid()) {
0317     const HGChefRecHitCollection *theHits = (fhRecHit.product());
0318     analyzeHGCalRecHit(theHits, 1, fhHitRefs);
0319   } else {
0320     edm::LogWarning("HGCalValid") << "No FH RecHit Found " << std::endl;
0321   }
0322 
0323   //accessing BH Rechit information
0324   const edm::Handle<HGChebRecHitCollection> &bhRecHit = iEvent.getHandle(bhRecHitToken_);
0325   if (bhRecHit.isValid()) {
0326     const HGChebRecHitCollection *theHits = (bhRecHit.product());
0327     analyzeHGCalRecHit(theHits, 2, bhHitRefs);
0328   } else {
0329     edm::LogWarning("HGCalValid") << "No BH RecHit Found " << std::endl;
0330   }
0331 }
0332 
0333 void HGCMissingRecHit::analyzeHGCalSimHit(edm::Handle<std::vector<PCaloHit>> const &simHits,
0334                                           int idet,
0335                                           std::map<unsigned int, HGCHitTuple> &hitRefs) {
0336   for (auto const &simHit : *simHits) {
0337     DetId id(simHit.id());
0338     bool ok(true), valid2(false);
0339     GlobalPoint p;
0340     bool valid1 = (hgcGeometry_[idet]->topology()).valid(id);
0341     std::ostringstream st1;
0342     if (id.det() == DetId::HGCalHSc) {
0343       st1 << HGCScintillatorDetId(id);
0344     } else if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi)) {
0345       st1 << HGCSiliconDetId(id);
0346     } else {
0347       st1 << "Not a Standard One";
0348     }
0349     if ((hgcCons_[idet]->waferHexagon8()) && ((id.det() == DetId::HGCalEE) || (id.det() == DetId::HGCalHSi))) {
0350       p = hgcGeometry_[idet]->getPosition(id);
0351       HGCSiliconDetId hid(id);
0352       valid2 = hgcCons_[idet]->isValidHex8(hid.layer(), hid.waferU(), hid.waferV(), hid.cellU(), hid.cellV(), true);
0353     } else if ((hgcCons_[idet]->tileTrapezoid()) && (id.det() == DetId::HGCalHSc)) {
0354       p = hgcGeometry_[idet]->getPosition(id);
0355       HGCScintillatorDetId hid(id);
0356       valid2 = hgcCons_[idet]->isValidTrap(hid.zside(), hid.layer(), hid.ring(), hid.iphi());
0357       edm::LogVerbatim("HGCalGeom") << "Scint " << HGCScintillatorDetId(id) << " position (" << p.x() << ", " << p.y()
0358                                     << ", " << p.z() << ") " << valid1 << ":" << valid2;
0359     } else {
0360       // This is an invalid cell
0361       ok = false;
0362       edm::LogVerbatim("HGCalError") << "Hit " << std::hex << id.rawId() << std::dec << " " << st1.str()
0363                                      << " in the wrong collection for detector " << idet << ":" << geometrySource_[idet]
0364                                      << " ***** ERROR *****";
0365     }
0366 
0367     edm::LogVerbatim("HGCalValid") << "SimHit: " << std::hex << id.rawId() << std::dec << " " << st1.str() << " Flag "
0368                                    << ok;
0369 
0370     if (ok) {
0371       float energy = simHit.energy();
0372 
0373       float energySum(energy);
0374       if (hitRefs.count(id) != 0)
0375         energySum += std::get<0>(hitRefs[id]);
0376       hitRefs[id] = std::make_tuple(energySum, p);
0377       edm::LogVerbatim("HGCalValid") << "Position = " << p << " Energy " << simHit.energy() << ":" << energySum;
0378     }
0379     if ((!valid1) || (!valid2))
0380       edm::LogVerbatim("HGCalMiss") << "Invalid SimHit " << st1.str() << " position = " << p << " perp " << p.perp()
0381                                     << " Validity flags " << valid1 << ":" << valid2;
0382   }
0383 }
0384 
0385 template <class T1>
0386 void HGCMissingRecHit::analyzeHGCalDigi(T1 const &theHits,
0387                                         int idet,
0388                                         std::map<unsigned int, HGCHitTuple> const &hitRefs) {
0389   std::vector<unsigned int> ids;
0390   for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0391     ids.emplace_back((it->id().rawId()));
0392     if (!(hgcGeometry_[idet]->topology()).valid(it->id())) {
0393       std::ostringstream st1;
0394       if ((it->id()).det() == DetId::HGCalHSc)
0395         st1 << HGCScintillatorDetId(it->id());
0396       else
0397         st1 << HGCSiliconDetId(it->id());
0398       edm::LogVerbatim("HGCalError") << "Invalid Hit " << st1.str();
0399     }
0400   }
0401   for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) {
0402     double eta = std::get<1>(it->second).eta();
0403     auto itr = std::find(ids.begin(), ids.end(), it->first);
0404     if (itr == ids.end()) {
0405       bool ok = (hgcGeometry_[idet]->topology()).valid(DetId(it->first));
0406       missedHitsDE_[idet]->Fill(std::get<0>(it->second));
0407       missedHitsDT_[idet]->Fill(eta);
0408       std::ostringstream st1;
0409       if (DetId(it->first).det() == DetId::HGCalHSc)
0410         st1 << HGCScintillatorDetId(it->first);
0411       else
0412         st1 << HGCSiliconDetId(it->first);
0413       edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str()
0414                                     << " SimHit (E = " << std::get<0>(it->second)
0415                                     << ", Position = " << std::get<1>(it->second) << ") Valid " << ok
0416                                     << " is missing in the Digi collection";
0417     } else {
0418       goodHitsDE_[idet]->Fill(std::get<0>(it->second));
0419       goodHitsDT_[idet]->Fill(eta);
0420     }
0421   }
0422 }
0423 
0424 template <class T1>
0425 void HGCMissingRecHit::analyzeHGCalRecHit(T1 const &theHits,
0426                                           int idet,
0427                                           std::map<unsigned int, HGCHitTuple> const &hitRefs) {
0428   std::vector<unsigned int> ids;
0429   for (auto it = theHits->begin(); it != theHits->end(); ++it) {
0430     ids.emplace_back((it->id().rawId()));
0431     if (!(hgcGeometry_[idet]->topology()).valid(it->id())) {
0432       std::ostringstream st1;
0433       if ((it->id()).det() == DetId::HGCalHSc)
0434         st1 << HGCScintillatorDetId(it->id());
0435       else
0436         st1 << HGCSiliconDetId(it->id());
0437       edm::LogVerbatim("HGCalError") << "Invalid Hit " << st1.str();
0438     }
0439   }
0440   for (auto it = hitRefs.begin(); it != hitRefs.end(); ++it) {
0441     double eta = std::get<1>(it->second).eta();
0442     auto itr = std::find(ids.begin(), ids.end(), it->first);
0443     if (itr == ids.end()) {
0444       bool ok = (hgcGeometry_[idet]->topology()).valid(DetId(it->first));
0445       missedHitsRE_[idet]->Fill(std::get<0>(it->second));
0446       missedHitsRT_[idet]->Fill(eta);
0447       std::ostringstream st1;
0448       if (DetId(it->first).det() == DetId::HGCalHSc)
0449         st1 << HGCScintillatorDetId(it->first);
0450       else
0451         st1 << HGCSiliconDetId(it->first);
0452       edm::LogVerbatim("HGCalMiss") << "Hit: " << std::hex << (it->first) << std::dec << " " << st1.str()
0453                                     << " SimHit (E = " << std::get<0>(it->second)
0454                                     << ", Position = " << std::get<1>(it->second) << ") Valid " << ok
0455                                     << " is missing in the RecHit collection";
0456     } else {
0457       goodHitsRE_[idet]->Fill(std::get<0>(it->second));
0458       goodHitsRT_[idet]->Fill(eta);
0459     }
0460   }
0461 }
0462 
0463 //define this as a plug-in
0464 DEFINE_FWK_MODULE(HGCMissingRecHit);