File indexing completed on 2025-03-13 02:31:41
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
0045 for (unsigned i = 0; i < 5; ++i) {
0046 if (nChannels_[i] == 0)
0047 nChannels_[i] = 1;
0048 }
0049
0050
0051
0052
0053
0054
0055 maxDepthHF_ = (maxDepthHF_ > 2 ? 2 : maxDepthHF_);
0056
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
0071
0072
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
0096
0097 int HcalRecHitsDQMClient::HcalRecHitsEndjob(const std::vector<MonitorElement *> &hcalMEs) {
0098 MonitorElement *Nhf = nullptr;
0099
0100
0101
0102
0103 std::vector<MonitorElement *> emap_depths;
0104
0105
0106
0107 std::vector<MonitorElement *> occupancy_maps;
0108 std::vector<std::string> occupancyID;
0109
0110
0111
0112 std::vector<MonitorElement *> emean_vs_ieta;
0113
0114
0115
0116
0117 std::vector<MonitorElement *> occupancy_vs_ieta;
0118 std::vector<std::string> occupancy_vs_ietaID;
0119
0120
0121
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
0129
0130 if (hcalMEs[ih]->getName() == "N_HF") {
0131 Nhf = hcalMEs[ih];
0132 continue;
0133 }
0134
0135
0136
0137
0138
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
0150
0151
0152
0153
0154
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
0172
0173
0174
0175
0176
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
0223 assert(Nhf);
0224 double nevtot = Nhf->getEntries();
0225
0226
0227 if (verbose_)
0228 std::cout << "nevtot : " << nevtot << std::endl;
0229
0230
0231 float fev = float(nevtot);
0232 double scaleBynevtot = 1 / fev;
0233
0234
0235
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
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 = 0;
0265 omatched = false;
0266
0267 for (; vsIetaIdx < occupancy_vs_ieta.size(); vsIetaIdx++) {
0268 if (occupancyID[occupancyIdx] == occupancy_vs_ietaID[vsIetaIdx]) {
0269 omatched = true;
0270 break;
0271 }
0272 }
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
0284
0285 if (omatched) {
0286
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;
0294
0295 phi_factor = phifactor(ieta);
0296
0297
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 }
0308
0309 int ietabin = occupancy_vs_ieta[vsIetaIdx]->getTH1F()->GetXaxis()->FindBin(float(ieta));
0310
0311
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 }
0318 }
0319 }
0320
0321
0322
0323
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);