HGCalValidation/plugins/HGCalSimHitsClient.cc

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