File indexing completed on 2024-04-06 12:07:26
0001 #include "DQM/GEM/interface/GEMDQMEfficiencyClientBase.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/Utilities/interface/isFinite.h"
0005
0006 #include "TEfficiency.h"
0007 #include "TPRegexp.h"
0008 #include <regex>
0009
0010 GEMDQMEfficiencyClientBase::GEMDQMEfficiencyClientBase(const edm::ParameterSet& ps)
0011 : kConfidenceLevel_(ps.getUntrackedParameter<double>("confidenceLevel")),
0012 kLogCategory_(ps.getUntrackedParameter<std::string>("logCategory")) {}
0013
0014
0015
0016
0017
0018
0019 std::tuple<bool, std::string, std::string, bool> GEMDQMEfficiencyClientBase::parseEfficiencySourceName(
0020 std::string name) {
0021
0022
0023 const bool success = TPRegexp("\\w+(?:_match)?_GE\\d1-(P|M)[0-9\\-]*").MatchB(name);
0024 if (not success) {
0025 return std::make_tuple(success, "", "", false);
0026 }
0027
0028 const std::string numerator_pattern = "_match";
0029 const auto numerator_pattern_start = name.find(numerator_pattern);
0030 const bool is_numerator = numerator_pattern_start != std::string::npos;
0031 if (is_numerator) {
0032
0033
0034 name.erase(numerator_pattern_start, numerator_pattern.length());
0035 }
0036
0037
0038
0039 const unsigned long last_pos = name.find_last_of('_');
0040
0041
0042 const std::string var_name = name.substr(0, last_pos);
0043
0044
0045 const std::string gem_name = name.substr(last_pos + 1);
0046 return std::make_tuple(success, var_name, gem_name, is_numerator);
0047 }
0048
0049 GEMDetId GEMDQMEfficiencyClientBase::parseGEMLabel(const std::string gem_label, const std::string delimiter) {
0050
0051
0052
0053
0054 int region = 0;
0055 int station = 0;
0056 int layer = 0;
0057 int chamber = 0;
0058 int ieta = 0;
0059
0060 std::vector<std::string> tokens;
0061
0062
0063 const std::regex re_station{"GE\\d1"};
0064 const std::regex re_region{"(P|M)"};
0065 const std::regex re_layer{"L\\d"};
0066 const std::regex re_chamber_layer{"\\d+L\\d"};
0067 const std::regex re_ieta{"E\\d+"};
0068
0069 std::string::size_type last_pos = gem_label.find_first_not_of(delimiter, 0);
0070 std::string::size_type pos = gem_label.find_first_of(delimiter, last_pos);
0071 while ((pos != std::string::npos) or (last_pos != std::string::npos)) {
0072 const std::string token = gem_label.substr(last_pos, pos - last_pos);
0073
0074 if (std::regex_match(token, re_region)) {
0075 region = (token == "P") ? 1 : -1;
0076
0077 } else if (std::regex_match(token, re_station)) {
0078 station = std::stoi(token.substr(2, 1));
0079
0080 } else if (std::regex_match(token, re_layer)) {
0081 layer = std::stoi(token.substr(1));
0082
0083 } else if (std::regex_match(token, re_chamber_layer)) {
0084 const unsigned long layer_prefix_pos = token.find('L');
0085 chamber = std::stoi(token.substr(0, layer_prefix_pos));
0086 layer = std::stoi(token.substr(layer_prefix_pos + 1));
0087
0088 } else if (std::regex_match(token, re_ieta)) {
0089 ieta = std::stoi(token.substr(1));
0090
0091 } else {
0092 edm::LogError(kLogCategory_) << "unknown pattern: " << gem_label << " --> " << token;
0093 }
0094 }
0095
0096 const GEMDetId id{region, 1, station, layer, chamber, ieta};
0097 return id;
0098 }
0099
0100 std::map<std::string, GEMDQMEfficiencyClientBase::MEPair> GEMDQMEfficiencyClientBase::makeEfficiencySourcePair(
0101 DQMStore::IBooker& ibooker, DQMStore::IGetter& igetter, const std::string& folder, const std::string prefix) {
0102 ibooker.setCurrentFolder(folder);
0103 igetter.setCurrentFolder(folder);
0104
0105 std::map<std::string, MEPair> me_pairs;
0106
0107 for (const std::string& name : igetter.getMEs()) {
0108
0109
0110 if (name.rfind(prefix, 0) != 0) {
0111
0112 continue;
0113 }
0114
0115 const std::string fullpath = folder + "/" + name;
0116 const MonitorElement* me = igetter.get(fullpath);
0117 if (me == nullptr) {
0118 edm::LogError(kLogCategory_) << "failed to get " << fullpath;
0119 continue;
0120 }
0121
0122 const auto [parsing_success, var_name, gem_name, is_matched] = parseEfficiencySourceName(name);
0123 if (not parsing_success) {
0124
0125 continue;
0126 }
0127
0128 const std::string key = var_name + "_" + gem_name;
0129
0130 if (me_pairs.find(key) == me_pairs.end()) {
0131 me_pairs[key] = {nullptr, nullptr};
0132 }
0133
0134 if (is_matched)
0135 me_pairs[key].first = me;
0136 else
0137 me_pairs[key].second = me;
0138 }
0139
0140
0141 for (auto it = me_pairs.cbegin(); it != me_pairs.cend();) {
0142 auto [me_numerator, me_denominator] = (*it).second;
0143
0144 bool okay = true;
0145 if (me_numerator == nullptr) {
0146 okay = false;
0147
0148 } else if (me_denominator == nullptr) {
0149 okay = false;
0150
0151 } else if (me_numerator->kind() != me_denominator->kind()) {
0152 okay = false;
0153 }
0154
0155
0156 if (okay) {
0157 it++;
0158
0159 } else {
0160 it = me_pairs.erase(it);
0161 }
0162 }
0163
0164 return me_pairs;
0165 }
0166
0167 void GEMDQMEfficiencyClientBase::setBins(TH1F* dst_hist, const TAxis* src_axis) {
0168 const int nbins = src_axis->GetNbins();
0169 if (src_axis->IsVariableBinSize()) {
0170 std::vector<double> edges;
0171 edges.reserve(nbins + 1);
0172
0173 for (int bin = 1; bin <= nbins; bin++) {
0174 edges.push_back(src_axis->GetBinLowEdge(bin));
0175 }
0176 edges.push_back(src_axis->GetBinUpEdge(nbins));
0177
0178 dst_hist->SetBins(nbins, &edges[0]);
0179
0180 } else {
0181 const double xlow = src_axis->GetBinLowEdge(1);
0182 const double xup = src_axis->GetBinUpEdge(nbins);
0183
0184 dst_hist->SetBins(nbins, xlow, xup);
0185 }
0186
0187 for (int bin = 1; bin <= nbins; bin++) {
0188 const TString label{src_axis->GetBinLabel(bin)};
0189 if (label.Length() > 0) {
0190 dst_hist->GetXaxis()->SetBinLabel(bin, label);
0191 }
0192 }
0193 }
0194
0195
0196
0197
0198
0199
0200
0201 bool GEMDQMEfficiencyClientBase::checkConsistency(const TH1& pass, const TH1& total) {
0202 if (pass.GetDimension() != total.GetDimension()) {
0203 edm::LogError(kLogCategory_) << "numerator and denominator have different dimensions: " << pass.GetName() << " & "
0204 << total.GetName();
0205 return false;
0206 }
0207
0208 if (not TEfficiency::CheckBinning(pass, total)) {
0209 edm::LogError(kLogCategory_) << "numerator and denominator have different binning: " << pass.GetName() << " & "
0210 << total.GetName();
0211 return false;
0212 }
0213
0214 if (not TEfficiency::CheckEntries(pass, total)) {
0215 edm::LogError(kLogCategory_) << "numerator and denominator do not have consistent bin contents " << pass.GetName()
0216 << " & " << total.GetName();
0217 return false;
0218 }
0219
0220 return true;
0221 }
0222
0223
0224 TH1F* GEMDQMEfficiencyClientBase::makeEfficiency(const TH1F* h_numerator,
0225 const TH1F* h_denominator,
0226 const char* name,
0227 const char* title) {
0228 if (h_numerator == nullptr) {
0229 edm::LogError(kLogCategory_) << "numerator is nullptr";
0230 return nullptr;
0231 }
0232
0233 if (h_denominator == nullptr) {
0234 edm::LogError(kLogCategory_) << "denominator is nulpptr";
0235 return nullptr;
0236 }
0237
0238 if (not checkConsistency(*h_numerator, *h_denominator)) {
0239 return nullptr;
0240 }
0241
0242 if (name == nullptr) {
0243 name = Form("eff_%s", h_denominator->GetName());
0244 }
0245
0246 if (title == nullptr) {
0247 title = h_denominator->GetTitle();
0248 }
0249
0250 const TAxis* x_axis = h_denominator->GetXaxis();
0251
0252
0253 TH1F* h_eff = new TH1F();
0254 h_eff->SetName(name);
0255 h_eff->SetTitle(title);
0256 h_eff->GetXaxis()->SetTitle(x_axis->GetTitle());
0257 h_eff->GetYaxis()->SetTitle("Efficiency");
0258 setBins(h_eff, h_denominator->GetXaxis());
0259
0260
0261 const int nbins = x_axis->GetNbins();
0262 for (int bin = 1; bin <= nbins; bin++) {
0263 const double passed = h_numerator->GetBinContent(bin);
0264 const double total = h_denominator->GetBinContent(bin);
0265
0266 if (total < 1) {
0267 continue;
0268 }
0269
0270 const double efficiency = passed / total;
0271 const double lower_boundary = TEfficiency::ClopperPearson(total, passed, kConfidenceLevel_, false);
0272 const double upper_boundary = TEfficiency::ClopperPearson(total, passed, kConfidenceLevel_, true);
0273 const double error = std::max(efficiency - lower_boundary, upper_boundary - efficiency);
0274
0275 h_eff->SetBinContent(bin, efficiency);
0276 h_eff->SetBinError(bin, error);
0277 }
0278
0279 return h_eff;
0280 }
0281
0282
0283 TH2F* GEMDQMEfficiencyClientBase::makeEfficiency(const TH2F* h_numerator,
0284 const TH2F* h_denominator,
0285 const char* name,
0286 const char* title) {
0287 if (h_numerator == nullptr) {
0288 edm::LogError(kLogCategory_) << "numerator is nullptr";
0289 return nullptr;
0290 }
0291
0292 if (h_denominator == nullptr) {
0293 edm::LogError(kLogCategory_) << "denominator is nulpptr";
0294 return nullptr;
0295 }
0296
0297 if (not checkConsistency(*h_numerator, *h_denominator)) {
0298 return nullptr;
0299 }
0300
0301 if (name == nullptr) {
0302 name = Form("eff_%s", h_denominator->GetName());
0303 }
0304
0305 if (title == nullptr) {
0306 title = h_denominator->GetTitle();
0307 }
0308
0309 TEfficiency eff(*h_numerator, *h_denominator);
0310 auto h_eff = dynamic_cast<TH2F*>(eff.CreateHistogram());
0311 h_eff->SetName(name);
0312 h_eff->SetTitle(title);
0313
0314 return h_eff;
0315 }
0316
0317
0318 TH1F* GEMDQMEfficiencyClientBase::projectHistogram(const TH2F* h_2d, const unsigned int on_which_axis) {
0319 if ((on_which_axis != TH1::kXaxis) and (on_which_axis != TH1::kYaxis)) {
0320 edm::LogError(kLogCategory_) << "invalid choice: " << on_which_axis << "."
0321 << " choose from [TH1::kXaxis (=1), TH1::kYaxis (=2)]";
0322 return nullptr;
0323 }
0324
0325 const bool on_x_axis = (on_which_axis == TH1::kXaxis);
0326
0327
0328 const TAxis* src_proj_axis = on_x_axis ? h_2d->GetXaxis() : h_2d->GetYaxis();
0329
0330 const TAxis* src_accum_axis = on_x_axis ? h_2d->GetYaxis() : h_2d->GetXaxis();
0331
0332 const TString prefix = on_x_axis ? "_proj_on_x" : "_proj_on_y";
0333 const TString name = h_2d->GetName() + prefix;
0334 const TString title = h_2d->GetTitle();
0335
0336 TH1F* h_proj = new TH1F();
0337 h_proj->SetName(name);
0338 h_proj->SetTitle(title);
0339 h_proj->GetXaxis()->SetTitle(src_proj_axis->GetTitle());
0340 setBins(h_proj, src_proj_axis);
0341
0342 for (int proj_bin = 1; proj_bin <= src_proj_axis->GetNbins(); proj_bin++) {
0343 double cumsum = 0.0;
0344 for (int accum_bin = 1; accum_bin <= src_accum_axis->GetNbins(); accum_bin++) {
0345 if (on_x_axis) {
0346 cumsum += h_2d->GetBinContent(proj_bin, accum_bin);
0347 } else {
0348 cumsum += h_2d->GetBinContent(accum_bin, proj_bin);
0349 }
0350 }
0351 h_proj->SetBinContent(proj_bin, cumsum);
0352 }
0353 h_proj->Sumw2();
0354 return h_proj;
0355 }
0356
0357 void GEMDQMEfficiencyClientBase::bookEfficiencyAuto(DQMStore::IBooker& ibooker,
0358 DQMStore::IGetter& igetter,
0359 const std::string& folder) {
0360 const std::map<std::string, MEPair> me_pairs = makeEfficiencySourcePair(ibooker, igetter, folder);
0361
0362 for (auto& [key, value] : me_pairs) {
0363 const auto& [me_numerator, me_denominator] = value;
0364
0365 const MonitorElement::Kind me_kind = me_numerator->kind();
0366 if (me_kind == MonitorElement::Kind::TH1F) {
0367 TH1F* h_numerator = me_numerator->getTH1F();
0368 if (h_numerator == nullptr) {
0369 edm::LogError(kLogCategory_) << "failed to get TH1F from h_numerator " << key;
0370 continue;
0371 }
0372
0373 TH1F* h_denominator = me_denominator->getTH1F();
0374 if (h_denominator == nullptr) {
0375 edm::LogError(kLogCategory_) << "failed to get TH1F from h_denominator" << key;
0376 continue;
0377 }
0378
0379 if (TH1F* eff = makeEfficiency(h_numerator, h_denominator)) {
0380 ibooker.book1D(eff->GetName(), eff);
0381
0382 } else {
0383
0384 continue;
0385 }
0386
0387 } else if (me_kind == MonitorElement::Kind::TH2F) {
0388 TH2F* h_numerator = me_numerator->getTH2F();
0389 if (h_numerator == nullptr) {
0390 edm::LogError(kLogCategory_) << "failed to get TH1F from h_numerator " << key;
0391 continue;
0392 }
0393
0394 TH2F* h_denominator = me_denominator->getTH2F();
0395 if (h_denominator == nullptr) {
0396 edm::LogError(kLogCategory_) << "failed to get TH1F from h_denominator" << key;
0397 continue;
0398 }
0399
0400 if (TH2F* eff = makeEfficiency(h_numerator, h_denominator)) {
0401 ibooker.book2D(eff->GetName(), eff);
0402
0403 } else {
0404
0405 continue;
0406 }
0407
0408 } else {
0409 edm::LogError(kLogCategory_) << "got an unepxected MonitorElement::Kind "
0410 << "0x" << std::hex << static_cast<int>(me_kind);
0411 continue;
0412 }
0413
0414 }
0415 }