MuonGEMHits/src/MuonGEMBaseHarvestor.cc

0001 #include "Validation/MuonGEMHits/interface/MuonGEMBaseHarvestor.h"
0002
0003 #include "TEfficiency.h"
0004
0005 MuonGEMBaseHarvestor::MuonGEMBaseHarvestor(const edm::ParameterSet& pset, std::string log_category)
0006     : kLogCategory_(log_category) {}
0007
0008 TProfile* MuonGEMBaseHarvestor::computeEfficiency(
0009     const TH1F& passed, const TH1F& total, const char* name, const char* title, Double_t confidence_level) {
0010   const TAxis* total_x = total.GetXaxis();
0011
0012   TProfile* eff_profile = new TProfile(name, title, total_x->GetNbins(), total_x->GetXmin(), total_x->GetXmax());
0013   eff_profile->GetXaxis()->SetTitle(total_x->GetTitle());
0014   eff_profile->GetYaxis()->SetTitle("#epsilon");
0015
0016   for (Int_t bin = 1; bin <= total.GetXaxis()->GetNbins(); bin++) {
0017     Double_t num_passed = passed.GetBinContent(bin);
0018     Double_t num_total = total.GetBinContent(bin);
0019
0020     if (num_passed > num_total) {
0021       edm::LogError(kLogCategory_) << "# of passed events > # of total events. "
0022                                    << "These numbers are not consistent." << std::endl;
0023       continue;
0024     }
0025
0026     if (num_total < 1) {
0027       eff_profile->SetBinEntries(bin, 0);
0028       continue;
0029     }
0030
0031     Double_t efficiency = num_passed / num_total;
0032
0033     Double_t lower_bound = TEfficiency::ClopperPearson(num_total, num_passed, confidence_level, false);
0034     Double_t upper_bound = TEfficiency::ClopperPearson(num_total, num_passed, confidence_level, true);
0035
0036     Double_t width = std::max(efficiency - lower_bound, upper_bound - efficiency);
0037     Double_t error = std::hypot(efficiency, width);
0038
0039     eff_profile->SetBinContent(bin, efficiency);
0040     eff_profile->SetBinError(bin, error);
0041     eff_profile->SetBinEntries(bin, 1);
0042   }
0043
0044   return eff_profile;
0045 }
0046
0047 TH2F* MuonGEMBaseHarvestor::computeEfficiency(const TH2F& passed,
0048                                               const TH2F& total,
0049                                               const char* name,
0050                                               const char* title) {
0051   TEfficiency eff(passed, total);
0052   TH2F* eff_hist = dynamic_cast<TH2F*>(eff.CreateHistogram());
0053   eff_hist->SetName(name);
0054   eff_hist->SetTitle(title);
0055
0056   const TAxis* total_x = total.GetXaxis();
0057   TAxis* eff_hist_x = eff_hist->GetXaxis();
0058   eff_hist_x->SetTitle(total_x->GetTitle());
0059   for (Int_t bin = 1; bin <= total.GetNbinsX(); bin++) {
0060     const char* label = total_x->GetBinLabel(bin);
0061     eff_hist_x->SetBinLabel(bin, label);
0062   }
0063
0064   const TAxis* total_y = total.GetYaxis();
0065   TAxis* eff_hist_y = eff_hist->GetYaxis();
0066   eff_hist_y->SetTitle(total_y->GetTitle());
0067   for (Int_t bin = 1; bin <= total.GetNbinsY(); bin++) {
0068     const char* label = total_y->GetBinLabel(bin);
0069     eff_hist_y->SetBinLabel(bin, label);
0070   }
0071
0072   return eff_hist;
0073 }
0074
0075 void MuonGEMBaseHarvestor::bookEff1D(DQMStore::IBooker& booker,
0076                                      DQMStore::IGetter& getter,
0077                                      const TString& passed_path,
0078                                      const TString& total_path,
0079                                      const TString& folder,
0080                                      const TString& eff_name,
0081                                      const TString& eff_title) {
0082   TH1F* passed = getElement<TH1F>(getter, passed_path);
0083   if (passed == nullptr) {
0084     edm::LogError(kLogCategory_) << "failed to get " << passed_path << std::endl;
0085     return;
0086   }
0087
0088   TH1F* total = getElement<TH1F>(getter, total_path);
0089   if (total == nullptr) {
0090     edm::LogError(kLogCategory_) << "failed to get " << total_path << std::endl;
0091     return;
0092   }
0093
0094   TProfile* eff = computeEfficiency(*passed, *total, eff_name, eff_title);
0095
0096   booker.setCurrentFolder(folder.Data());
0097   booker.bookProfile(eff_name, eff);
0098 }
0099
0100 void MuonGEMBaseHarvestor::bookEff2D(DQMStore::IBooker& booker,
0101                                      DQMStore::IGetter& getter,
0102                                      const TString& passed_path,
0103                                      const TString& total_path,
0104                                      const TString& folder,
0105                                      const TString& eff_name,
0106                                      const TString& eff_title) {
0107   TH2F* passed = getElement<TH2F>(getter, passed_path);
0108   if (passed == nullptr) {
0109     edm::LogError(kLogCategory_) << "failed to get " << passed_path << std::endl;
0110     return;
0111   }
0112
0113   TH2F* total = getElement<TH2F>(getter, total_path);
0114   if (total == nullptr) {
0115     edm::LogError(kLogCategory_) << "failed to get " << total_path << std::endl;
0116     return;
0117   }
0118
0119   TH2F* eff = computeEfficiency(*passed, *total, eff_name, eff_title);
0120   booker.setCurrentFolder(folder.Data());
0121   booker.book2D(eff_name, eff);
0122 }