HGCalValidation/plugins/HGCalRecHitsClient.cc

0001 #include <unistd.h>
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "DataFormats/Math/interface/LorentzVector.h"
0008
0009 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011
0012 #include "FWCore/Framework/interface/Run.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018
0019 #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h"
0020 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0021
0022 class HGCalRecHitsClient : public DQMEDHarvester {
0023 private:
0024   //member data
0025   const std::string nameDetector_;
0026   const int verbosity_;
0027   const edm::ESGetToken<HGCalDDDConstants, IdealGeometryRecord> ddc_token_;
0028   unsigned int layers_;
0029
0030 public:
0031   explicit HGCalRecHitsClient(const edm::ParameterSet &);
0032   ~HGCalRecHitsClient() override = default;
0033
0034   void beginRun(const edm::Run &run, const edm::EventSetup &c) override;
0035   void dqmEndJob(DQMStore::IBooker &ib, DQMStore::IGetter &ig) override;
0036   virtual void runClient_(DQMStore::IBooker &ib, DQMStore::IGetter &ig);
0037
0038   int recHitsEndjob(const std::vector<MonitorElement *> &hgcalMEs);
0039 };
0040
0041 HGCalRecHitsClient::HGCalRecHitsClient(const edm::ParameterSet &iConfig)
0042     : nameDetector_(iConfig.getParameter<std::string>("DetectorName")),
0043       verbosity_(iConfig.getUntrackedParameter<int>("Verbosity", 0)),
0044       ddc_token_(esConsumes<HGCalDDDConstants, IdealGeometryRecord, edm::Transition::BeginRun>(
0045           edm::ESInputTag{"", nameDetector_})) {}
0046
0047 void HGCalRecHitsClient::beginRun(const edm::Run &run, const edm::EventSetup &iSetup) {
0048   const HGCalDDDConstants &hgcons_ = iSetup.getData(ddc_token_);
0049   layers_ = hgcons_.layers(true);
0050 }
0051
0052 void HGCalRecHitsClient::dqmEndJob(DQMStore::IBooker &ib, DQMStore::IGetter &ig) { runClient_(ib, ig); }
0053
0054 void HGCalRecHitsClient::runClient_(DQMStore::IBooker &ib, DQMStore::IGetter &ig) {
0055   ig.setCurrentFolder("/");
0056   if (verbosity_ > 0)
0057     edm::LogVerbatim("HGCalValidation") << "\nrunClient";
0058   std::vector<MonitorElement *> hgcalMEs;
0059   std::vector<std::string> fullDirPath = ig.getSubdirs();
0060
0061   for (unsigned int i = 0; i < fullDirPath.size(); i++) {
0062     if (verbosity_ > 0)
0063       edm::LogVerbatim("HGCalValidation") << "\nfullPath: " << fullDirPath.at(i);
0064     ig.setCurrentFolder(fullDirPath.at(i));
0065     std::vector<std::string> fullSubDirPath = ig.getSubdirs();
0066
0067     for (unsigned int j = 0; j < fullSubDirPath.size(); j++) {
0068       if (verbosity_ > 1)
0069         edm::LogVerbatim("HGCalValidation") << "fullSubPath: " << fullSubDirPath.at(j);
0070       std::string nameDirectory = "HGCAL/HGCalRecHitsV/" + nameDetector_;
0071       if (strcmp(fullSubDirPath.at(j).c_str(), nameDirectory.c_str()) == 0) {
0072         hgcalMEs = ig.getContents(fullSubDirPath.at(j));
0073         if (verbosity_ > 1)
0074           edm::LogVerbatim("HGCalValidation") << "hgcalMES size : " << hgcalMEs.size();
0075         if (!recHitsEndjob(hgcalMEs))
0076           edm::LogWarning("HGCalValidation") << "\nError in RecHitsEndjob!";
0077       }
0078     }
0079   }
0080 }
0081
0082 int HGCalRecHitsClient::recHitsEndjob(const std::vector<MonitorElement *> &hgcalMEs) {
0083   std::vector<MonitorElement *> energy_;
0084   std::vector<MonitorElement *> EtaPhi_Plus_;
0085   std::vector<MonitorElement *> EtaPhi_Minus_;
0086   std::vector<MonitorElement *> HitOccupancy_Plus_;
0087   std::vector<MonitorElement *> HitOccupancy_Minus_;
0088   std::vector<MonitorElement *> MeanHitOccupancy_Plus_;
0089   std::vector<MonitorElement *> MeanHitOccupancy_Minus_;
0090
0091   std::ostringstream name;
0092   double nevent;
0093   int nbinsx, nbinsy;
0094
0095   for (unsigned int ilayer = 0; ilayer < layers_; ilayer++) {
0096     name.str("");
0097     name << "energy_layer_" << ilayer;
0098     for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0099       if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0100         energy_.push_back(hgcalMEs[ih]);
0101       }
0102     }
0103
0104     //normalization
0105     nevent = energy_.at(ilayer)->getEntries();
0106     nbinsx = energy_.at(ilayer)->getNbinsX();
0107     for (int i = 1; i <= nbinsx; i++) {
0108       double binValue = energy_.at(ilayer)->getBinContent(i) / nevent;
0109       energy_.at(ilayer)->setBinContent(i, binValue);
0110     }
0111
0112     //EtaPhi 2d plots
0113     name.str("");
0114     name << "EtaPhi_Plus_"
0115          << "layer_" << ilayer;
0116     for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0117       if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0118         EtaPhi_Plus_.push_back(hgcalMEs[ih]);
0119       }
0120     }
0121
0122     name.str("");
0123     name << "EtaPhi_Minus_"
0124          << "layer_" << ilayer;
0125     for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0126       if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0127         EtaPhi_Minus_.push_back(hgcalMEs[ih]);
0128       }
0129     }
0130
0131     //normalization EtaPhi
0132     nevent = EtaPhi_Plus_.at(ilayer)->getEntries();
0133     nbinsx = EtaPhi_Plus_.at(ilayer)->getNbinsX();
0134     nbinsy = EtaPhi_Plus_.at(ilayer)->getNbinsY();
0135     for (int i = 1; i <= nbinsx; ++i) {
0136       for (int j = 1; j <= nbinsy; ++j) {
0137         double binValue = EtaPhi_Plus_.at(ilayer)->getBinContent(i, j) / nevent;
0138         EtaPhi_Plus_.at(ilayer)->setBinContent(i, j, binValue);
0139       }
0140     }
0141
0142     nevent = EtaPhi_Minus_.at(ilayer)->getEntries();
0143     nbinsx = EtaPhi_Minus_.at(ilayer)->getNbinsX();
0144     nbinsy = EtaPhi_Plus_.at(ilayer)->getNbinsY();
0145     for (int i = 1; i <= nbinsx; ++i) {
0146       for (int j = 1; j <= nbinsy; ++j) {
0147         double binValue = EtaPhi_Minus_.at(ilayer)->getBinContent(i, j) / nevent;
0148         EtaPhi_Minus_.at(ilayer)->setBinContent(i, j, binValue);
0149       }
0150     }
0151
0152     //HitOccupancy
0153     name.str("");
0154     name << "HitOccupancy_Plus_layer_" << ilayer;
0155     for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0156       if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0157         HitOccupancy_Plus_.push_back(hgcalMEs[ih]);
0158       }
0159     }
0160
0161     name.str("");
0162     name << "HitOccupancy_Minus_layer_" << ilayer;
0163     for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0164       if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0165         HitOccupancy_Minus_.push_back(hgcalMEs[ih]);
0166       }
0167     }
0168
0169     //normalization of hit occupancy histos
0170     nevent = HitOccupancy_Plus_.at(ilayer)->getEntries();
0171     nbinsx = HitOccupancy_Plus_.at(ilayer)->getNbinsX();
0172     for (int i = 1; i <= nbinsx; ++i) {
0173       double binValue = HitOccupancy_Plus_.at(ilayer)->getBinContent(i) / nevent;
0174       HitOccupancy_Plus_.at(ilayer)->setBinContent(i, binValue);
0175     }
0176
0177     nevent = HitOccupancy_Minus_.at(ilayer)->getEntries();
0178     nbinsx = HitOccupancy_Minus_.at(ilayer)->getNbinsX();
0179     for (int i = 1; i <= nbinsx; ++i) {
0180       double binValue = HitOccupancy_Minus_.at(ilayer)->getBinContent(i) / nevent;
0181       HitOccupancy_Minus_.at(ilayer)->setBinContent(i, binValue);
0182     }
0183
0184   }  //loop over layers
0185
0186   name.str("");
0187   name << "SUMOfRecHitOccupancy_Plus";
0188   for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0189     if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0190       MeanHitOccupancy_Plus_.push_back(hgcalMEs[ih]);
0191       unsigned int indx = MeanHitOccupancy_Plus_.size() - 1;
0192       for (int ilayer = 0; ilayer < static_cast<int>(layers_); ++ilayer) {
0193         double meanVal = HitOccupancy_Plus_.at(ilayer)->getMean();
0194         MeanHitOccupancy_Plus_[indx]->setBinContent(ilayer + 1, meanVal);
0195       }
0196       break;
0197     }
0198   }
0199
0200   name.str("");
0201   name << "SUMOfRecHitOccupancy_Plus";
0202   for (unsigned int ih = 0; ih < hgcalMEs.size(); ih++) {
0203     if (strcmp(hgcalMEs[ih]->getName().c_str(), name.str().c_str()) == 0) {
0204       MeanHitOccupancy_Minus_.push_back(hgcalMEs[ih]);
0205       unsigned indx = MeanHitOccupancy_Minus_.size() - 1;
0206       for (int ilayer = 0; ilayer < static_cast<int>(layers_); ++ilayer) {
0207         double meanVal = HitOccupancy_Minus_.at(ilayer)->getMean();
0208         MeanHitOccupancy_Minus_[indx]->setBinContent(ilayer + 1, meanVal);
0209       }
0210       break;
0211     }
0212   }
0213
0214   return 1;
0215 }
0216
0217 #include "FWCore/Framework/interface/MakerMacros.h"
0218
0219 DEFINE_FWK_MODULE(HGCalRecHitsClient);