Hcal/src/HcalRecHitsDQMClient.cc

0001 #include "DQMOffline/Hcal/interface/HcalRecHitsDQMClient.h"
0002 #include "FWCore/Framework/interface/MakerMacros.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/Utilities/interface/Transition.h"
0008 #include "FWCore/Utilities/interface/ESInputTag.h"
0009
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011
0012 HcalRecHitsDQMClient::HcalRecHitsDQMClient(const edm::ParameterSet &iConfig)
0013     : conf_(iConfig),
0014       hcalDDDRecConstantsToken_{esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()},
0015       caloGeometryRunToken_{esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()} {
0016   outputFile_ = iConfig.getUntrackedParameter<std::string>("outputFile", "myfile.root");
0017   debug_ = false;
0018   verbose_ = false;
0019   dirName_ = iConfig.getParameter<std::string>("DQMDirName");
0020 }
0021
0022 HcalRecHitsDQMClient::~HcalRecHitsDQMClient() {}
0023
0024 void HcalRecHitsDQMClient::beginJob() {}
0025
0026 void HcalRecHitsDQMClient::beginRun(const edm::Run &run, const edm::EventSetup &iSetup) {
0027   HcalDDDRecConstants const &hcons = iSetup.getData(hcalDDDRecConstantsToken_);
0028   maxDepthHB_ = hcons.getMaxDepth(0);
0029   maxDepthHE_ = hcons.getMaxDepth(1);
0030   maxDepthHF_ = std::max(hcons.getMaxDepth(2), 1);
0031   maxDepthHO_ = hcons.getMaxDepth(3);
0032
0033   CaloGeometry const &geometry = iSetup.getData(caloGeometryRunToken_);
0034   const std::vector<DetId> &hbCells = geometry.getValidDetIds(DetId::Hcal, HcalBarrel);
0035   const std::vector<DetId> &heCells = geometry.getValidDetIds(DetId::Hcal, HcalEndcap);
0036   const std::vector<DetId> &hoCells = geometry.getValidDetIds(DetId::Hcal, HcalOuter);
0037   const std::vector<DetId> &hfCells = geometry.getValidDetIds(DetId::Hcal, HcalForward);
0038
0039   nChannels_[1] = hbCells.size();
0040   nChannels_[2] = heCells.size();
0041   nChannels_[3] = hoCells.size();
0042   nChannels_[4] = hfCells.size();
0043   nChannels_[0] = nChannels_[1] + nChannels_[2] + nChannels_[3] + nChannels_[4];
0044   // avoid divide by zero
0045   for (unsigned i = 0; i < 5; ++i) {
0046     if (nChannels_[i] == 0)
0047       nChannels_[i] = 1;
0048   }
0049
0050   // std::cout << "Channels HB:" << nChannels_[1] << " HE:" << nChannels_[2] <<
0051   // " HO:" << nChannels_[3] << " HF:" << nChannels_[4] << std::endl;
0052
0053   // We hardcode the HF depths because in the dual readout configuration,
0054   // rechits are not defined for depths 3&4
0055   maxDepthHF_ = (maxDepthHF_ > 2 ? 2 : maxDepthHF_);  // We reatin the dynamic possibility
0056                                                       // that HF might have 0 or 1 depths
0057
0058   maxDepthAll_ = (maxDepthHB_ + maxDepthHO_ > maxDepthHE_ ? maxDepthHB_ + maxDepthHO_ : maxDepthHE_);
0059   maxDepthAll_ = (maxDepthAll_ > maxDepthHF_ ? maxDepthAll_ : maxDepthHF_);
0060 }
0061
0062 void HcalRecHitsDQMClient::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0063   igetter.setCurrentFolder(dirName_);
0064
0065   if (verbose_)
0066     std::cout << "\nrunClient" << std::endl;
0067
0068   std::vector<MonitorElement *> hcalMEs;
0069
0070   // Since out folders are fixed to three, we can just go over these three
0071   // folders i.e., CaloTowersD/CaloTowersTask, HcalRecHitsD/HcalRecHitTask,
0072   // NoiseRatesV/NoiseRatesTask.
0073   std::vector<std::string> fullPathHLTFolders = igetter.getSubdirs();
0074   for (unsigned int i = 0; i < fullPathHLTFolders.size(); i++) {
0075     if (verbose_)
0076       std::cout << "\nfullPath: " << fullPathHLTFolders[i] << std::endl;
0077     igetter.setCurrentFolder(fullPathHLTFolders[i]);
0078
0079     std::vector<std::string> fullSubPathHLTFolders = igetter.getSubdirs();
0080     for (unsigned int j = 0; j < fullSubPathHLTFolders.size(); j++) {
0081       if (verbose_)
0082         std::cout << "fullSub: " << fullSubPathHLTFolders[j] << std::endl;
0083
0084       if (strcmp(fullSubPathHLTFolders[j].c_str(), "HcalRecHitsD/HcalRecHitTask") == 0) {
0085         hcalMEs = igetter.getContents(fullSubPathHLTFolders[j]);
0086         if (verbose_)
0087           std::cout << "hltMES size : " << hcalMEs.size() << std::endl;
0088         if (!HcalRecHitsEndjob(hcalMEs))
0089           std::cout << "\nError in HcalRecHitsEndjob!" << std::endl << std::endl;
0090       }
0091     }
0092   }
0093 }
0094
0095 // called after entering the HcalRecHitsD/HcalRecHitTask directory
0096 // hcalMEs are within that directory
0097 int HcalRecHitsDQMClient::HcalRecHitsEndjob(const std::vector<MonitorElement *> &hcalMEs) {
0098   MonitorElement *Nhf = nullptr;
0099
0100   // Search for emap histograms, and collect them into this vector
0101   // All subdtectors are plotted together in these histograms. We only need to
0102   // look for different depths
0103   std::vector<MonitorElement *> emap_depths;
0104
0105   // This vector is filled occupancy_maps identified by both subdetector and
0106   // depth
0107   std::vector<MonitorElement *> occupancy_maps;
0108   std::vector<std::string> occupancyID;
0109
0110   // This vector is filled with emean_vs_ieta histograms, they are divided by
0111   // both subdetector and depth
0112   std::vector<MonitorElement *> emean_vs_ieta;
0113
0114   // These are the only histograms filled in this module; however, the
0115   // histograms are created empty in HcalRecHitsAnalyzer occupancy_vs_ieta,
0116   // divided by both subdetector and depth
0117   std::vector<MonitorElement *> occupancy_vs_ieta;
0118   std::vector<std::string> occupancy_vs_ietaID;
0119
0120   // RecHit_StatusWord & RecHit_Aux_StatusWord
0121   // Divided by subdectector
0122   std::vector<MonitorElement *> RecHit_StatusWord;
0123   std::vector<float> RecHit_StatusWord_Channels;
0124   std::vector<MonitorElement *> RecHit_Aux_StatusWord;
0125   std::vector<float> RecHit_Aux_StatusWord_Channels;
0126
0127   for (unsigned int ih = 0; ih < hcalMEs.size(); ih++) {
0128     // N_HF is not special, it is just convient to get the total number of
0129     // events The number of entries in N_HF is equal to the number of events
0130     if (hcalMEs[ih]->getName() == "N_HF") {
0131       Nhf = hcalMEs[ih];
0132       continue;
0133     }
0134
0135     // ***********************
0136     // * We fill the various MonitorElement vectors by searching for a matching
0137     // substring
0138     // * The methods that are used are agnostic to the ordering of vectors
0139     // ***********************
0140
0141     if (hcalMEs[ih]->getName().find("emap_depth") != std::string::npos) {
0142       emap_depths.push_back(hcalMEs[ih]);
0143       continue;
0144     }
0145
0146     if (hcalMEs[ih]->getName().find("occupancy_map_H") != std::string::npos) {
0147       occupancy_maps.push_back(hcalMEs[ih]);
0148
0149       // Use occupancyID to save the subdetector and depth information
0150       // This will help preserve both indifference to vector ordering and
0151       // specific details of the detector topology The position in occupancyID
0152       // must correspond to the histogram position in occupancy_maps
0153
0154       // Save the string after "occupancy_map_"
0155
0156       std::string prefix = "occupancy_map_";
0157
0158       occupancyID.push_back(hcalMEs[ih]->getName().substr(prefix.size()));
0159
0160       continue;
0161     }
0162
0163     if (hcalMEs[ih]->getName().find("emean_vs_ieta_H") != std::string::npos) {
0164       emean_vs_ieta.push_back(hcalMEs[ih]);
0165       continue;
0166     }
0167
0168     if (hcalMEs[ih]->getName().find("occupancy_vs_ieta_H") != std::string::npos) {
0169       occupancy_vs_ieta.push_back(hcalMEs[ih]);
0170
0171       // Use occupancy_vs_ietaID to save the subdetector and depth information
0172       // This will help preserve both indifference to vector ordering and
0173       // specific details of the detector topology The position in occupancyID
0174       // must correspond to the histogram position in occupancy_vs_ieta
0175
0176       // Save the string after "occupancy_vs_ieta_"
0177
0178       std::string prefix = "occupancy_vs_ieta_";
0179
0180       occupancy_vs_ietaID.push_back(hcalMEs[ih]->getName().substr(prefix.size()));
0181
0182       continue;
0183     }
0184
0185     if (hcalMEs[ih]->getName().find("HcalRecHitTask_RecHit_StatusWord_H") != std::string::npos) {
0186       RecHit_StatusWord.push_back(hcalMEs[ih]);
0187
0188       if (hcalMEs[ih]->getName().find("HB") != std::string::npos) {
0189         RecHit_StatusWord_Channels.push_back((float)nChannels_[1]);
0190       } else if (hcalMEs[ih]->getName().find("HE") != std::string::npos) {
0191         RecHit_StatusWord_Channels.push_back((float)nChannels_[2]);
0192       } else if (hcalMEs[ih]->getName().find("H0") != std::string::npos) {
0193         RecHit_StatusWord_Channels.push_back((float)nChannels_[3]);
0194       } else if (hcalMEs[ih]->getName().find("HF") != std::string::npos) {
0195         RecHit_StatusWord_Channels.push_back((float)nChannels_[4]);
0196       } else {
0197         RecHit_StatusWord_Channels.push_back(1.);
0198       }
0199
0200       continue;
0201     }
0202
0203     if (hcalMEs[ih]->getName().find("HcalRecHitTask_RecHit_Aux_StatusWord_H") != std::string::npos) {
0204       RecHit_Aux_StatusWord.push_back(hcalMEs[ih]);
0205
0206       if (hcalMEs[ih]->getName().find("HB") != std::string::npos) {
0207         RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[1]);
0208       } else if (hcalMEs[ih]->getName().find("HE") != std::string::npos) {
0209         RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[2]);
0210       } else if (hcalMEs[ih]->getName().find("H0") != std::string::npos) {
0211         RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[3]);
0212       } else if (hcalMEs[ih]->getName().find("HF") != std::string::npos) {
0213         RecHit_Aux_StatusWord_Channels.push_back((float)nChannels_[4]);
0214       } else {
0215         RecHit_Aux_StatusWord_Channels.push_back(1.);
0216       }
0217
0218       continue;
0219     }
0220   }
0221
0222   // mean energies and occupancies evaluation
0223
0224   double nevtot = Nhf->getEntries();  // Use the number of entries in the Nhf histogram to
0225                                       // give the total number of events
0226
0227   if (verbose_)
0228     std::cout << "nevtot : " << nevtot << std::endl;
0229
0230   // emap histograms are scaled by the number of events
0231   float fev = float(nevtot);
0232   double scaleBynevtot = 1 / fev;
0233
0234   // In this and the following histogram vectors, recognize that the for-loop
0235   // index does not have to correspond to any particular depth
0236   for (unsigned int depthIdx = 0; depthIdx < emap_depths.size(); depthIdx++) {
0237     int nx = emap_depths[depthIdx]->getNbinsX();
0238     int ny = emap_depths[depthIdx]->getNbinsY();
0239
0240     float cnorm;
0241     float enorm;
0242
0243     for (int i = 1; i <= nx; i++) {
0244       for (int j = 1; j <= ny; j++) {
0245         cnorm = emap_depths[depthIdx]->getBinContent(i, j) * scaleBynevtot;
0246         enorm = emap_depths[depthIdx]->getBinError(i, j) * scaleBynevtot;
0247         emap_depths[depthIdx]->setBinContent(i, j, cnorm);
0248         emap_depths[depthIdx]->setBinError(i, j, enorm);
0249       }
0250     }
0251   }
0252
0253   // occupancy_maps & matched occupancy_vs_ieta
0254
0255   bool omatched = false;
0256
0257   for (unsigned int occupancyIdx = 0; occupancyIdx < occupancy_maps.size(); occupancyIdx++) {
0258     int nx = occupancy_maps[occupancyIdx]->getNbinsX();
0259     int ny = occupancy_maps[occupancyIdx]->getNbinsY();
0260
0261     float cnorm;
0262     float enorm;
0263
0264     unsigned int vsIetaIdx = occupancy_vs_ieta.size();
0265     omatched = false;
0266
0267     for (vsIetaIdx = 0; vsIetaIdx < occupancy_vs_ieta.size(); vsIetaIdx++) {
0268       if (occupancyID[occupancyIdx] == occupancy_vs_ietaID[vsIetaIdx]) {
0269         omatched = true;
0270         break;
0271       }
0272     }  // match occupancy_vs_ieta histogram
0273
0274     for (int i = 1; i <= nx; i++) {
0275       for (int j = 1; j <= ny; j++) {
0276         cnorm = occupancy_maps[occupancyIdx]->getBinContent(i, j) * scaleBynevtot;
0277         enorm = occupancy_maps[occupancyIdx]->getBinError(i, j) * scaleBynevtot;
0278         occupancy_maps[occupancyIdx]->setBinContent(i, j, cnorm);
0279         occupancy_maps[occupancyIdx]->setBinError(i, j, enorm);
0280       }
0281     }
0282
0283     // Fill occupancy_vs_ieta
0284
0285     if (omatched) {
0286       // We run over all of the ieta values
0287       for (int ieta = -41; ieta <= 41; ieta++) {
0288         float phi_factor = 1.;
0289         float sumphi = 0.;
0290         float sumphie = 0.;
0291
0292         if (ieta == 0)
0293           continue;  // ieta=0 is not defined
0294
0295         phi_factor = phifactor(ieta);
0296
0297         // the rechits occupancy map defines iphi as 0..71
0298         for (int iphi = 0; iphi <= 71; iphi++) {
0299           int binIeta = occupancy_maps[occupancyIdx]->getTH2F()->GetXaxis()->FindBin(float(ieta));
0300           int binIphi = occupancy_maps[occupancyIdx]->getTH2F()->GetYaxis()->FindBin(float(iphi));
0301
0302           float content = occupancy_maps[occupancyIdx]->getBinContent(binIeta, binIphi);
0303           float econtent = occupancy_maps[occupancyIdx]->getBinError(binIeta, binIphi);
0304
0305           sumphi += content;
0306           sumphie += econtent * econtent;
0307         }  // for loop over phi
0308
0309         int ietabin = occupancy_vs_ieta[vsIetaIdx]->getTH1F()->GetXaxis()->FindBin(float(ieta));
0310
0311         // fill occupancies vs ieta
0312         cnorm = sumphi / phi_factor;
0313         enorm = sqrt(sumphie) / phi_factor;
0314         occupancy_vs_ieta[vsIetaIdx]->setBinContent(ietabin, cnorm);
0315         occupancy_vs_ieta[vsIetaIdx]->setBinError(ietabin, enorm);
0316
0317       }  // Fill occupancy_vs_ieta
0318     }    // if omatched
0319   }
0320
0321   // Status Word
0322   // Normalized by number of events and by number of channels per subdetector as
0323   // well
0324
0325   for (unsigned int StatusWordIdx = 0; StatusWordIdx < RecHit_StatusWord.size(); StatusWordIdx++) {
0326     int nx = RecHit_StatusWord[StatusWordIdx]->getNbinsX();
0327
0328     float cnorm;
0329     float enorm;
0330
0331     for (int i = 1; i <= nx; i++) {
0332       cnorm = RecHit_StatusWord[StatusWordIdx]->getBinContent(i) * scaleBynevtot /
0333               RecHit_StatusWord_Channels[StatusWordIdx];
0334       enorm =
0335           RecHit_StatusWord[StatusWordIdx]->getBinError(i) * scaleBynevtot / RecHit_StatusWord_Channels[StatusWordIdx];
0336       RecHit_StatusWord[StatusWordIdx]->setBinContent(i, cnorm);
0337       RecHit_StatusWord[StatusWordIdx]->setBinError(i, enorm);
0338     }
0339   }
0340
0341   for (unsigned int AuxStatusWordIdx = 0; AuxStatusWordIdx < RecHit_Aux_StatusWord.size(); AuxStatusWordIdx++) {
0342     int nx = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getNbinsX();
0343
0344     float cnorm;
0345     float enorm;
0346
0347     for (int i = 1; i <= nx; i++) {
0348       cnorm = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getBinContent(i) * scaleBynevtot /
0349               RecHit_Aux_StatusWord_Channels[AuxStatusWordIdx];
0350       enorm = RecHit_Aux_StatusWord[AuxStatusWordIdx]->getBinError(i) * scaleBynevtot /
0351               RecHit_Aux_StatusWord_Channels[AuxStatusWordIdx];
0352       RecHit_Aux_StatusWord[AuxStatusWordIdx]->setBinContent(i, cnorm);
0353       RecHit_Aux_StatusWord[AuxStatusWordIdx]->setBinError(i, enorm);
0354     }
0355   }
0356
0357   return 1;
0358 }
0359
0360 float HcalRecHitsDQMClient::phifactor(int ieta) {
0361   float phi_factor_;
0362
0363   if (ieta >= -20 && ieta <= 20) {
0364     phi_factor_ = 72.;
0365   } else {
0366     if (ieta >= 40 || ieta <= -40) {
0367       phi_factor_ = 18.;
0368     } else {
0369       phi_factor_ = 36.;
0370     }
0371   }
0372
0373   return phi_factor_;
0374 }
0375
0376 DEFINE_FWK_MODULE(HcalRecHitsDQMClient);