File indexing completed on 2024-04-06 12:32:33
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
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
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
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 }
0113
0114
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
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
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 }
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);