Back to home page

Project CMSSW displayed by LXR

 
 

    


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

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