Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-26 22:39:25

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 // Returns a tuple of
0015 //   - a boolean indicating whether the parsing is successful or not
0016 //   - name of a variable used in the efficiency monitoring
0017 //   - GEM subdetector name like GE11-P-L1
0018 //   - a boolean indicating whether the name is a numerator name.
0019 std::tuple<bool, std::string, std::string, bool> GEMDQMEfficiencyClientBase::parseEfficiencySourceName(
0020     std::string name) {
0021   // NOTE This expression must be consistent with TODO
0022   // TODO use regex
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     // keep a delimiter between a variable name and a GEM name
0033     // e.g. 'muon_pt_matched_GE11-L1' --> 'muon_pt_GE11-L1'
0034     name.erase(numerator_pattern_start, numerator_pattern.length());
0035   }
0036   // find the position of the delimiter.
0037   // Because variable name can has "_", find the last one.
0038   // NOTE The GEM name must not contains "_"
0039   const unsigned long last_pos = name.find_last_of('_');
0040 
0041   // "muon_pt"
0042   const std::string var_name = name.substr(0, last_pos);
0043 
0044   // "GE11-L1"
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   // GE11-P
0051   // GE11-P-L1
0052   // GE11-P-E1
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   // static const?
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     // If name doesn't start with prefix
0109     // The default prefix is empty string.
0110     if (name.rfind(prefix, 0) != 0) {
0111       // TODO LogDebug
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       // TODO LogDebug
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   // remove invalid pairs
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     // anyways, move on to the next one
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 // Returns a boolean indicating whether the numerator and the denominator are
0196 // consistent.
0197 //
0198 // TEfficiency::CheckConsistency raises errors and leads to an exception.
0199 // So, the efficiency client will skip inconsitent two histograms.
0200 // https://github.com/root-project/root/blob/v6-24-06/hist/hist/src/TEfficiency.cxx#L1494-L1512
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 // MonitorElement doesn't support TGraphAsymmErrors
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   // create an empty TProfile for storing efficiencies and uncertainties.
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   // efficiency calculation
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 // FIXME TH2D::ProjectionX looks buggy
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   // on which axis is the histogram projected?
0328   const TAxis* src_proj_axis = on_x_axis ? h_2d->GetXaxis() : h_2d->GetYaxis();
0329   // along which axis do the entries accumulate?
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         // makeEfficiency will report the error.
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         // makeEfficiency will report the error.
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   }  // me_pairs
0415 }