Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-22 23:03:56

0001 // system include files
0002 #include <cmath>
0003 #include <fstream>
0004 #include <iostream>
0005 #include <map>
0006 #include <string>
0007 #include <vector>
0008 
0009 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0010 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
0011 #include "DataFormats/DetId/interface/DetId.h"
0012 
0013 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015 
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0023 
0024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0025 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0026 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0027 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0028 
0029 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0030 #include "TVector3.h"
0031 
0032 class HGCalRecHitValidation : public DQMEDAnalyzer {
0033 public:
0034   struct energysum {
0035     energysum() { e15 = e25 = e50 = e100 = e250 = e1000 = 0.0; }
0036     double e15, e25, e50, e100, e250, e1000;
0037   };
0038 
0039   struct HitsInfo {
0040     HitsInfo() {
0041       x = y = z = time = energy = phi = eta = 0.0;
0042       layer = 0;
0043     }
0044     float x, y, z, time, energy, phi, eta;
0045     float layer;
0046   };
0047 
0048   explicit HGCalRecHitValidation(const edm::ParameterSet&);
0049   ~HGCalRecHitValidation() override = default;
0050 
0051   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0052   void dqmBeginRun(const edm::Run&, const edm::EventSetup&) override;
0053   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0054   void analyze(const edm::Event&, const edm::EventSetup&) override;
0055 
0056 private:
0057   template <class T1, class T2>
0058   void recHitValidation(DetId& detId, int layer, const T1* geom, T2 it);
0059   void fillHitsInfo();
0060   void fillHitsInfo(HitsInfo& hits);
0061   void fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer);
0062 
0063   // ----------member data ---------------------------
0064   const std::string nameDetector_;
0065   const int verbosity_;
0066   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geometry_token_;
0067   const edm::ESGetToken<HGCalGeometry, IdealGeometryRecord> geometry_beginRun_token_;
0068   const edm::EDGetTokenT<HGCRecHitCollection> recHitSource_;
0069   unsigned int layers_;
0070   int firstLayer_;
0071   std::map<int, int> OccupancyMap_plus;
0072   std::map<int, int> OccupancyMap_minus;
0073 
0074   std::vector<MonitorElement*> EtaPhi_Plus_;
0075   std::vector<MonitorElement*> EtaPhi_Minus_;
0076   std::vector<MonitorElement*> energy_;
0077   std::vector<MonitorElement*> HitOccupancy_Plus_;
0078   std::vector<MonitorElement*> HitOccupancy_Minus_;
0079   MonitorElement* MeanHitOccupancy_Plus_;
0080   MonitorElement* MeanHitOccupancy_Minus_;
0081 };
0082 
0083 HGCalRecHitValidation::HGCalRecHitValidation(const edm::ParameterSet& iConfig)
0084     : nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
0085       verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0086       geometry_token_(esConsumes(edm::ESInputTag("", nameDetector_))),
0087       geometry_beginRun_token_(esConsumes<HGCalGeometry, IdealGeometryRecord, edm::Transition::BeginRun>(
0088           edm::ESInputTag("", nameDetector_))),
0089       recHitSource_(consumes(iConfig.getParameter<edm::InputTag>("RecHitSource"))),
0090       firstLayer_(1) {}
0091 
0092 void HGCalRecHitValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0093   edm::ParameterSetDescription desc;
0094   desc.add<std::string>("DetectorName", "HGCalEESensitive");
0095   desc.add<edm::InputTag>("RecHitSource", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
0096   desc.addUntracked<int>("Verbosity", 0);
0097   descriptions.add("hgcalRecHitValidationEE", desc);
0098 }
0099 
0100 void HGCalRecHitValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0101   OccupancyMap_plus.clear();
0102   OccupancyMap_minus.clear();
0103 
0104   bool ok(true);
0105   unsigned int ntot(0), nused(0);
0106   const edm::ESHandle<HGCalGeometry>& geomHandle = iSetup.getHandle(geometry_token_);
0107   if (!geomHandle.isValid()) {
0108     edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry "
0109                                         << "Object for " << nameDetector_;
0110   } else {
0111     const HGCalGeometry* geom0 = geomHandle.product();
0112     int geomType = ((geom0->topology().waferHexagon8()) ? 1 : ((geom0->topology().tileTrapezoid()) ? 2 : 0));
0113 
0114     const edm::Handle<HGCRecHitCollection>& theRecHitContainers = iEvent.getHandle(recHitSource_);
0115     if (theRecHitContainers.isValid()) {
0116       if (verbosity_ > 0)
0117         edm::LogVerbatim("HGCalValidation")
0118             << nameDetector_ << " with " << theRecHitContainers->size() << " element(s)";
0119       for (const auto& it : *(theRecHitContainers.product())) {
0120         ntot++;
0121         nused++;
0122         DetId detId = it.id();
0123         int layer = ((geomType == 0)
0124                          ? HGCalDetId(detId).layer()
0125                          : ((geomType == 1) ? HGCSiliconDetId(detId).layer() : HGCScintillatorDetId(detId).layer()));
0126         recHitValidation(detId, layer, geom0, &it);
0127       }
0128     } else {
0129       ok = false;
0130       edm::LogVerbatim("HGCalValidation") << "HGCRecHitCollection Handle "
0131                                           << "does not exist !!!";
0132     }
0133   }
0134   if (ok)
0135     fillHitsInfo();
0136   if (verbosity_ > 0)
0137     edm::LogVerbatim("HGCalValidation") << "Event " << iEvent.id().event() << " with " << ntot << " total and " << nused
0138                                         << " used recHits";
0139 }
0140 
0141 template <class T1, class T2>
0142 void HGCalRecHitValidation::recHitValidation(DetId& detId, int layer, const T1* geom, T2 it) {
0143   const GlobalPoint& global = geom->getPosition(detId);
0144   double energy = it->energy();
0145 
0146   float globalx = global.x();
0147   float globaly = global.y();
0148   float globalz = global.z();
0149 
0150   HitsInfo hinfo;
0151   hinfo.energy = energy;
0152   hinfo.x = globalx;
0153   hinfo.y = globaly;
0154   hinfo.z = globalz;
0155   hinfo.layer = layer - firstLayer_;
0156   hinfo.phi = global.phi();
0157   hinfo.eta = global.eta();
0158 
0159   if (verbosity_ > 1)
0160     edm::LogVerbatim("HGCalValidation") << "--------------------------   gx = " << globalx << " gy = " << globaly
0161                                         << " gz = " << globalz << " phi = " << hinfo.phi << " eta = " << hinfo.eta
0162                                         << " lay = " << hinfo.layer;
0163 
0164   fillHitsInfo(hinfo);
0165 
0166   if (hinfo.eta > 0)
0167     fillOccupancyMap(OccupancyMap_plus, hinfo.layer);
0168   else
0169     fillOccupancyMap(OccupancyMap_minus, hinfo.layer);
0170 }
0171 
0172 void HGCalRecHitValidation::fillOccupancyMap(std::map<int, int>& OccupancyMap, int layer) {
0173   if (OccupancyMap.find(layer) != OccupancyMap.end())
0174     OccupancyMap[layer]++;
0175   else
0176     OccupancyMap[layer] = 1;
0177 }
0178 
0179 void HGCalRecHitValidation::fillHitsInfo() {
0180   for (auto const& itr : OccupancyMap_plus) {
0181     int layer = itr.first;
0182     int occupancy = itr.second;
0183     HitOccupancy_Plus_.at(layer)->Fill(occupancy);
0184   }
0185 
0186   for (auto const& itr : OccupancyMap_minus) {
0187     int layer = itr.first;
0188     int occupancy = itr.second;
0189     HitOccupancy_Minus_.at(layer)->Fill(occupancy);
0190   }
0191 }
0192 
0193 void HGCalRecHitValidation::fillHitsInfo(HitsInfo& hits) {
0194   unsigned int ilayer = hits.layer;
0195   energy_.at(ilayer)->Fill(hits.energy);
0196 
0197   EtaPhi_Plus_.at(ilayer)->Fill(hits.eta, hits.phi);
0198   EtaPhi_Minus_.at(ilayer)->Fill(hits.eta, hits.phi);
0199 }
0200 
0201 void HGCalRecHitValidation::dqmBeginRun(const edm::Run&, const edm::EventSetup& iSetup) {
0202   const edm::ESHandle<HGCalGeometry>& geomHandle = iSetup.getHandle(geometry_beginRun_token_);
0203   if (!geomHandle.isValid()) {
0204     edm::LogVerbatim("HGCalValidation") << "Cannot get valid HGCalGeometry "
0205                                         << "Object for " << nameDetector_;
0206   } else {
0207     const HGCalGeometry* geom = geomHandle.product();
0208     const HGCalDDDConstants& hgcons_ = geom->topology().dddConstants();
0209     layers_ = hgcons_.layers(true);
0210     firstLayer_ = hgcons_.firstLayer();
0211   }
0212 }
0213 
0214 void HGCalRecHitValidation::bookHistograms(DQMStore::IBooker& iB, edm::Run const&, edm::EventSetup const&) {
0215   iB.setCurrentFolder("HGCAL/HGCalRecHitsV/" + nameDetector_);
0216   std::ostringstream histoname;
0217   for (unsigned int il = 0; il < layers_; ++il) {
0218     int ilayer = firstLayer_ + static_cast<int>(il);
0219     auto istr1 = std::to_string(ilayer);
0220     while (istr1.size() < 2) {
0221       istr1.insert(0, "0");
0222     }
0223     histoname.str("");
0224     histoname << "HitOccupancy_Plus_layer_" << istr1;
0225     HitOccupancy_Plus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Plus", 100, 0, 10000));
0226     histoname.str("");
0227     histoname << "HitOccupancy_Minus_layer_" << istr1;
0228     HitOccupancy_Minus_.push_back(iB.book1D(histoname.str().c_str(), "RecHitOccupancy_Minus", 100, 0, 10000));
0229 
0230     histoname.str("");
0231     histoname << "EtaPhi_Plus_"
0232               << "layer_" << istr1;
0233     EtaPhi_Plus_.push_back(iB.book2D(histoname.str().c_str(), "Occupancy", 31, 1.45, 3.0, 72, -CLHEP::pi, CLHEP::pi));
0234     histoname.str("");
0235     histoname << "EtaPhi_Minus_"
0236               << "layer_" << istr1;
0237     EtaPhi_Minus_.push_back(
0238         iB.book2D(histoname.str().c_str(), "Occupancy", 31, -3.0, -1.45, 72, -CLHEP::pi, CLHEP::pi));
0239 
0240     histoname.str("");
0241     histoname << "energy_layer_" << istr1;
0242     energy_.push_back(iB.book1D(histoname.str().c_str(), "energy_", 500, 0, 1));
0243   }  //loop over layers ends here
0244 
0245   histoname.str("");
0246   histoname << "SUMOfRecHitOccupancy_Plus";
0247   MeanHitOccupancy_Plus_ =
0248       iB.book1D(histoname.str().c_str(), "SUMOfRecHitOccupancy_Plus", layers_, -0.5, layers_ - 0.5);
0249   histoname.str("");
0250   histoname << "SUMOfRecHitOccupancy_Minus";
0251   MeanHitOccupancy_Minus_ =
0252       iB.book1D(histoname.str().c_str(), "SUMOfRecHitOccupancy_Minus", layers_, -0.5, layers_ - 0.5);
0253 }
0254 
0255 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0256 
0257 #include "FWCore/Framework/interface/MakerMacros.h"
0258 
0259 //define this as a plug-in
0260 DEFINE_FWK_MODULE(HGCalRecHitValidation);