HcalHits/src/HcalSimHitsClient.cc

0001 #include "FWCore/Framework/interface/MakerMacros.h"
0002 #include "Validation/HcalHits/interface/HcalSimHitsClient.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/Run.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008
0009 #include "DQMServices/Core/interface/DQMStore.h"
0010
0011 HcalSimHitsClient::HcalSimHitsClient(const edm::ParameterSet &iConfig)
0012     : dirName_(iConfig.getParameter<std::string>("DQMDirName")),
0013       verbose_(iConfig.getUntrackedParameter<bool>("Verbosity", false)),
0014       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {}
0015
0016 void HcalSimHitsClient::beginRun(edm::Run const &run, edm::EventSetup const &c) {
0017   const auto &pHRNDC = c.getData(tok_HRNDC_);
0018   hcons = &pHRNDC;
0019   maxDepthHB_ = hcons->getMaxDepth(0);
0020   maxDepthHE_ = hcons->getMaxDepth(1);
0021   maxDepthHF_ = hcons->getMaxDepth(2);
0022   maxDepthHO_ = hcons->getMaxDepth(3);
0023
0024   edm::LogVerbatim("HitsValidationHcal") << " Maximum Depths HB:" << maxDepthHB_ << " HE:" << maxDepthHE_
0025                                          << " HO:" << maxDepthHO_ << " HF:" << maxDepthHF_;
0026 }
0027
0028 void HcalSimHitsClient::dqmEndJob(DQMStore::IBooker &ib, DQMStore::IGetter &ig) { runClient_(ib, ig); }
0029
0030 void HcalSimHitsClient::runClient_(DQMStore::IBooker &ib, DQMStore::IGetter &ig) {
0031   ig.setCurrentFolder(dirName_);
0032
0033   if (verbose_)
0034     edm::LogVerbatim("HitsValidationHcal") << "runClient";
0035
0036   std::vector<MonitorElement *> hcalMEs;
0037
0038   std::vector<std::string> fullPathHLTFolders = ig.getSubdirs();
0039   for (unsigned int i = 0; i < fullPathHLTFolders.size(); i++) {
0040     if (verbose_)
0041       edm::LogVerbatim("HitsValidationHcal") << "fullPath: " << fullPathHLTFolders[i];
0042     ig.setCurrentFolder(fullPathHLTFolders[i]);
0043
0044     std::vector<std::string> fullSubPathHLTFolders = ig.getSubdirs();
0045     for (unsigned int j = 0; j < fullSubPathHLTFolders.size(); j++) {
0046       if (verbose_)
0047         edm::LogVerbatim("HitsValidationHcal") << "fullSub: " << fullSubPathHLTFolders[j];
0048
0049       if (strcmp(fullSubPathHLTFolders[j].c_str(), "HcalHitsV/SimHitsValidationHcal") == 0) {
0050         hcalMEs = ig.getContents(fullSubPathHLTFolders[j]);
0051         if (verbose_)
0052           edm::LogVerbatim("HitsValidationHcal") << "hltMES size : " << hcalMEs.size();
0053         if (!SimHitsEndjob(hcalMEs))
0054           edm::LogWarning("HitsValidationHcal") << "Error in SimhitEndjob!";
0055       }
0056     }
0057   }
0058 }
0059
0060 // called after entering the  directory
0061 // hcalMEs are within that directory
0062 int HcalSimHitsClient::SimHitsEndjob(const std::vector<MonitorElement *> &hcalMEs) {
0063   std::vector<std::string> divisions = getHistogramTypes();
0064   MonitorElement *Occupancy_map[nTime][divisions.size()];
0065   MonitorElement *Energy[nType1], *Time_weighteden[nType1];
0066   MonitorElement *HitEnergyvsieta[divisions.size()], *HitTimevsieta[divisions.size()];
0067
0068   std::string time[nTime] = {"25", "50", "100", "250"};
0069   std::string detdivision[nType1] = {"HB", "HE", "HF", "HO"};
0070   char name[40], name1[40], name2[40], name3[40], name4[40];
0071
0072   for (int k = 0; k < nType1; k++) {
0073     Energy[k] = nullptr;
0074     Time_weighteden[k] = nullptr;
0075     for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
0076       sprintf(name1, "Energy_%s", detdivision[k].c_str());
0077       sprintf(name2, "Time_Enweighted_%s", detdivision[k].c_str());
0078       if (strcmp(hcalMEs[ih]->getName().c_str(), name1) == 0) {
0079         Energy[k] = hcalMEs[ih];
0080       }
0081       if (strcmp(hcalMEs[ih]->getName().c_str(), name2) == 0) {
0082         Time_weighteden[k] = hcalMEs[ih];
0083       }
0084     }
0085   }
0086
0087   for (int i = 0; i < nTime; i++) {
0088     for (unsigned int j = 0; j < divisions.size(); j++) {
0089       for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
0090         sprintf(name, "HcalHitE%s%s", time[i].c_str(), divisions[j].c_str());
0091         if (strcmp(hcalMEs[ih]->getName().c_str(), name) == 0) {
0092           Occupancy_map[i][j] = hcalMEs[ih];
0093         }
0094       }
0095     }
0096   }
0097
0098   for (unsigned int k = 0; k < divisions.size(); k++) {
0099     HitEnergyvsieta[k] = nullptr;
0100     HitTimevsieta[k] = nullptr;
0101     for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
0102       sprintf(name3, "HcalHitEta%s", divisions[k].c_str());
0103       sprintf(name4, "HcalHitTimeAEta%s", divisions[k].c_str());
0104       if (strcmp(hcalMEs[ih]->getName().c_str(), name3) == 0) {
0105         HitEnergyvsieta[k] = hcalMEs[ih];
0106       }
0107       if (strcmp(hcalMEs[ih]->getName().c_str(), name4) == 0) {
0108         HitTimevsieta[k] = hcalMEs[ih];
0109       }
0110     }
0111   }
0112
0113   // mean energy
0114
0115   double nevent = Energy[0]->getEntries();
0116   if (verbose_)
0117     edm::LogVerbatim("HitsValidationHcal") << "nevent : " << nevent;
0118
0119   float cont[nTime][divisions.size()];
0120   float en[nType1], tme[nType1];
0121   float hitenergy[divisions.size()], hittime[divisions.size()];
0122   float fev = float(nevent);
0123
0124   for (int dettype = 0; dettype < nType1; dettype++) {
0125     int nx1 = Energy[dettype]->getNbinsX();
0126     for (int i = 0; i <= nx1; i++) {
0127       en[dettype] = Energy[dettype]->getBinContent(i) / fev;
0128       Energy[dettype]->setBinContent(i, en[dettype]);
0129     }
0130     int nx2 = Time_weighteden[dettype]->getNbinsX();
0131     for (int i = 0; i <= nx2; i++) {
0132       tme[dettype] = Time_weighteden[dettype]->getBinContent(i) / fev;
0133       Time_weighteden[dettype]->setBinContent(i, tme[dettype]);
0134     }
0135   }
0136
0137   for (unsigned int dettype = 0; dettype < divisions.size(); dettype++) {
0138     int nx1 = HitEnergyvsieta[dettype]->getNbinsX();
0139     for (int i = 0; i <= nx1; i++) {
0140       hitenergy[dettype] = HitEnergyvsieta[dettype]->getBinContent(i) / fev;
0141       HitEnergyvsieta[dettype]->setBinContent(i, hitenergy[dettype]);
0142     }
0143     int nx2 = HitTimevsieta[dettype]->getNbinsX();
0144     for (int i = 0; i <= nx2; i++) {
0145       hittime[dettype] = HitTimevsieta[dettype]->getBinContent(i) / fev;
0146       HitTimevsieta[dettype]->setBinContent(i, hittime[dettype]);
0147     }
0148   }
0149
0150   for (int itime = 0; itime < nTime; itime++) {
0151     for (unsigned int det = 0; det < divisions.size(); det++) {
0152       int ny = Occupancy_map[itime][det]->getNbinsY();
0153       int nx = Occupancy_map[itime][det]->getNbinsX();
0154       for (int i = 1; i < nx + 1; i++) {
0155         for (int j = 1; j < ny + 1; j++) {
0156           cont[itime][det] = Occupancy_map[itime][det]->getBinContent(i, j) / fev;
0157           Occupancy_map[itime][det]->setBinContent(i, j, cont[itime][det]);
0158         }
0159       }
0160     }
0161   }
0162
0163   return 1;
0164 }
0165
0166 std::vector<std::string> HcalSimHitsClient::getHistogramTypes() {
0167   int maxDepth = std::max(maxDepthHB_, maxDepthHE_);
0168   if (verbose_)
0169     edm::LogVerbatim("HitsValidationHcal") << "Max depth 1st step:: " << maxDepth;
0170   maxDepth = std::max(maxDepth, maxDepthHF_);
0171   if (verbose_)
0172     edm::LogVerbatim("HitsValidationHcal") << "Max depth 2nd step:: " << maxDepth;
0173   maxDepth = std::max(maxDepth, maxDepthHO_);
0174   if (verbose_)
0175     edm::LogVerbatim("HitsValidationHcal") << "Max depth 3rd step:: " << maxDepth;
0176   std::vector<std::string> divisions;
0177   char name1[20];
0178
0179   // first overall Hcal
0180   for (int depth = 0; depth < maxDepth; ++depth) {
0181     sprintf(name1, "HC%d", depth);
0182     divisions.push_back(std::string(name1));
0183   }
0184   // HB
0185   for (int depth = 0; depth < maxDepthHB_; ++depth) {
0186     sprintf(name1, "HB%d", depth);
0187     divisions.push_back(std::string(name1));
0188   }
0189   // HE
0190   for (int depth = 0; depth < maxDepthHE_; ++depth) {
0191     sprintf(name1, "HE%d+z", depth);
0192     divisions.push_back(std::string(name1));
0193     sprintf(name1, "HE%d-z", depth);
0194     divisions.push_back(std::string(name1));
0195   }
0196   // HO
0197   {
0198     int depth = maxDepthHO_;
0199     sprintf(name1, "HO%d", depth);
0200     divisions.push_back(std::string(name1));
0201   }
0202   // HF (first absorber, then different types of abnormal hits)
0203   std::string hfty1[4] = {"A", "W", "B", "J"};
0204   for (int k = 0; k < 4; ++k) {
0205     for (int depth = 0; depth < maxDepthHF_; ++depth) {
0206       sprintf(name1, "HF%s%d+z", hfty1[k].c_str(), depth);
0207       divisions.push_back(std::string(name1));
0208       sprintf(name1, "HF%s%d-z", hfty1[k].c_str(), depth);
0209       divisions.push_back(std::string(name1));
0210     }
0211   }
0212   return divisions;
0213 }
0214
0215 DEFINE_FWK_MODULE(HcalSimHitsClient);