HcalDigis/src/HcalDigisClient.cc

0001 #include <Validation/HcalDigis/interface/HcalDigisClient.h>
0002 // -*- C++ -*-
0003 //
0004 // Package:    HcalDigisClient
0005 // Class:      HcalDigisClient
0006 //
0007 /**\class HcalDigisClient HcalDigisClient.cc Validation/HcalDigisClient/src/HcalDigisClient.cc
0008
0009  Description: [one line class summary]
0010
0011  Implementation:
0012      [Notes on implementation]
0013  */
0014 //
0015 // Original Author:  Ali Fahim,22 R-013,+41227672649,
0016 //         Created:  Wed Mar 23 22:54:28 CET 2011
0017 //
0018 //
0019
0020 // system include files
0021
0022 HcalDigisClient::HcalDigisClient(const edm::ParameterSet& iConfig) {
0023   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "HcalDigisClient.root");
0024   dirName_ = iConfig.getParameter<std::string>("DQMDirName");
0025   msm_ = new std::map<std::string, MonitorElement*>();
0026 }
0027
0028 HcalDigisClient::~HcalDigisClient() { delete msm_; }
0029
0030 void HcalDigisClient::runClient(DQMStore::IBooker& ib, DQMStore::IGetter& ig) {
0031   ig.setCurrentFolder(dirName_);
0032   std::vector<MonitorElement*> hcalMEs;
0033   // Since out folders are fixed to three, we can just go over these three folders
0034   // i.e., CaloTowersV/CaloTowersTask, HcalRecHitsV/HcalRecHitTask, NoiseRatesV/NoiseRatesTask.
0035   std::vector<std::string> fullPathHLTFolders = ig.getSubdirs();
0036   for (unsigned int i = 0; i < fullPathHLTFolders.size(); i++) {
0037     ig.setCurrentFolder(fullPathHLTFolders[i]);
0038     std::vector<std::string> fullSubPathHLTFolders = ig.getSubdirs();
0039     for (unsigned int j = 0; j < fullSubPathHLTFolders.size(); j++) {
0040       if (strcmp(fullSubPathHLTFolders[j].c_str(), "HcalDigisV/HcalDigiTask") == 0) {
0041         hcalMEs = ig.getContents(fullSubPathHLTFolders[j]);
0042         ig.setCurrentFolder("HcalDigisV/HcalDigiTask");
0043         if (!HcalDigisEndjob(hcalMEs, "HB", ib))
0044           edm::LogError("HcalDigisClient") << "Error in HcalDigisEndjob! HB";
0045         if (!HcalDigisEndjob(hcalMEs, "HE", ib))
0046           edm::LogError("HcalDigisClient") << "Error in HcalDigisEndjob! HE";
0047         if (!HcalDigisEndjob(hcalMEs, "HO", ib))
0048           edm::LogError("HcalDigisClient") << "Error in HcalDigisEndjob! HO";
0049         if (!HcalDigisEndjob(hcalMEs, "HF", ib))
0050           edm::LogError("HcalDigisClient") << "Error in HcalDigisEndjob! HF";
0051       }
0052     }
0053   }
0054 }
0055
0056 int HcalDigisClient::HcalDigisEndjob(const std::vector<MonitorElement*>& hcalMEs,
0057                                      std::string subdet_,
0058                                      DQMStore::IBooker& ib) {
0059   using namespace std;
0060   string strtmp;
0061
0062   MonitorElement* nevtot(nullptr);
0063
0064   std::vector<MonitorElement*> ieta_iphi_occupancy_maps;
0065   std::vector<std::string> depthID;
0066
0067   edm::LogVerbatim("HcalDigisClient") << " Number of histos " << hcalMEs.size();
0068
0069   for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
0070     if (hcalMEs[ih]->getName() == "nevtot") {
0071       nevtot = hcalMEs[ih];
0072       continue;
0073     }
0074
0075     //We search the occupancy maps corresponding to this subdetector
0076     if ((hcalMEs[ih]->getName().find("HcalDigiTask_ieta_iphi_occupancy_map_depth") != std::string::npos) &&
0077         (hcalMEs[ih]->getName().find(subdet_) != std::string::npos)) {
0078       ieta_iphi_occupancy_maps.push_back(hcalMEs[ih]);
0079
0080       std::string start = "depth";
0081       std::string end = "_H";
0082
0083       int position = hcalMEs[ih]->getName().find(start) + start.length();
0084       int length = hcalMEs[ih]->getName().find(end) - position;
0085
0086       depthID.push_back(hcalMEs[ih]->getName().substr(position, length));
0087
0088       continue;
0089     }
0090   }
0091
0092   if (hcalMEs.empty()) {
0093     edm::LogError("HcalDigisClient") << "No nevtot or maps histo found...";
0094     return 0;
0095   }
0096   if (!nevtot) {
0097     edm::LogError("HcalDigisClient") << "No nevtot histoo found...";
0098     return 0;
0099   }
0100   if (ieta_iphi_occupancy_maps.empty()) {
0101     edm::LogError("HcalDigisClient") << "No maps histos found...";
0102     return 0;
0103   }
0104
0105   int ev = nevtot->getEntries();
0106
0107   if (ev <= 0) {
0108     edm::LogError("HcalDigisClient") << "normalization factor <= 0!";
0109     return 0;
0110   }
0111
0112   float fev = (float)nevtot->getEntries();
0113
0114   int depths = ieta_iphi_occupancy_maps.size();
0115
0116   HistLim ietaLim(85, -42.5, 42.5);
0117
0118   for (int depth = 1; depth <= depths; depth++) {
0119     strtmp = "HcalDigiTask_occupancy_vs_ieta_depth" + str(depth) + "_" + subdet_;
0120     book1D(ib, strtmp, ietaLim);
0121   }
0122
0123   std::vector<float> sumphi(depths, 0);
0124   std::vector<float> sumphie(depths, 0);
0125
0126   float phi_factor;
0127   float cnorm;
0128   float enorm;
0129
0130   for (int depth = 1; depth <= depths; depth++) {
0131     int nx = ieta_iphi_occupancy_maps[depth - 1]->getNbinsX();
0132     int ny = ieta_iphi_occupancy_maps[depth - 1]->getNbinsY();
0133
0134     for (int i = 1; i <= nx; i++) {
0135       for (int j = 1; j <= ny; j++) {
0136         // occupancies
0137         cnorm = ieta_iphi_occupancy_maps[depth - 1]->getBinContent(i, j) / fev;
0138         enorm = ieta_iphi_occupancy_maps[depth - 1]->getBinError(i, j) / fev;
0139         ieta_iphi_occupancy_maps[depth - 1]->setBinContent(i, j, cnorm);
0140         ieta_iphi_occupancy_maps[depth - 1]->setBinError(i, j, enorm);
0141
0142       }  //for loop over NbinsYU
0143     }  //for loop over NbinsX
0144   }  //for loop over the occupancy maps
0145
0146   for (int i = 1; i <= 82; i++) {
0147     int ieta = i - 42;  // -41 -1, 0 40
0148     if (ieta >= 0)
0149       ieta += 1;  // -41 -1, 1 41  - to make it detector-like
0150
0151     if (ieta >= -20 && ieta <= 20) {
0152       phi_factor = 72.;
0153     } else {
0154       if (ieta >= 40 || ieta <= -40)
0155         phi_factor = 18.;
0156       else
0157         phi_factor = 36.;
0158     }
0159
0160     //zero the sumphi and sumphie vector at the start of each ieta ring
0161     sumphi.assign(depths, 0);
0162     sumphie.assign(depths, 0);
0163
0164     for (int iphi = 1; iphi <= 72; iphi++) {
0165       for (int depth = 1; depth <= depths; depth++) {
0166         int binIeta = ieta_iphi_occupancy_maps[depth - 1]->getTH2F()->GetXaxis()->FindBin(ieta);
0167         int binIphi = ieta_iphi_occupancy_maps[depth - 1]->getTH2F()->GetYaxis()->FindBin(iphi);
0168
0169         float content = ieta_iphi_occupancy_maps[depth - 1]->getBinContent(binIeta, binIphi);
0170         float econtent = ieta_iphi_occupancy_maps[depth - 1]->getBinError(binIeta, binIphi);
0171
0172         sumphi[depth - 1] += content;
0173         sumphie[depth - 1] += econtent * econtent;
0174
0175       }  //for loop over depths
0176     }  //for loop over phi
0177
0178     //double deta = double(ieta);
0179
0180     // occupancies vs ieta
0181     for (int depth = 1; depth <= depths; depth++) {
0182       strtmp = "HcalDigiTask_occupancy_vs_ieta_depth" + depthID[depth - 1] + "_" + subdet_;
0183       MonitorElement* ME = msm_->find(strtmp)->second;
0184       int ietabin = ME->getTH1F()->GetXaxis()->FindBin(float(ieta));
0185
0186       if (sumphi[depth - 1] > 1.e-30) {
0187         cnorm = sumphi[depth - 1] / phi_factor;
0188         enorm = sqrt(sumphie[depth - 1]) / phi_factor;
0189         ME->setBinContent(ietabin, cnorm);
0190         ME->setBinError(ietabin, enorm);
0191       }
0192     }
0193   }  // end of i-loop
0194
0195   return 1;
0196 }
0197
0198 HcalDigisClient::MonitorElement* HcalDigisClient::monitor(std::string name) {
0199   if (!msm_->count(name))
0200     return nullptr;
0201   else
0202     return msm_->find(name)->second;
0203 }
0204
0205 std::string HcalDigisClient::str(int x) {
0206   std::stringstream out;
0207   out << x;
0208   return out.str();
0209 }
0210
0211 double HcalDigisClient::integralMETH2D(MonitorElement* ME, int i0, int i1, int j0, int j1) {
0212   double sum(0);
0213   for (int i = i0; i <= i1; i++) {
0214     for (int j = j0; j <= j1; j++) {
0215       sum += ME->getBinContent(i, j);
0216     }
0217   }
0218
0219   return sum;
0220 }
0221
0222 void HcalDigisClient::scaleMETH2D(MonitorElement* ME, double s) {
0223   int nx = ME->getNbinsX();
0224   int ny = ME->getNbinsY();
0225
0226   double content(0);
0227   double error(0);
0228   for (int i = 1; i <= nx; i++) {
0229     for (int j = 1; j <= ny; j++) {
0230       content = ME->getBinContent(i, j);
0231       error = ME->getBinError(i, j);
0232       content *= s;
0233       error *= s;
0234       ME->setBinContent(i, j, content);
0235       ME->setBinError(i, j, error);
0236     }
0237   }
0238 }
0239
0240 //define this as a plug-in
0241 DEFINE_FWK_MODULE(HcalDigisClient);