Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18:47 

0001 #include "HLTriggerOffline/Btag/interface/HLTBTagHarvestingAnalyzer.h"
0002 
0003 HLTBTagHarvestingAnalyzer::HLTBTagHarvestingAnalyzer(const edm::ParameterSet &iConfig) {
0004   // getParameter
0005   mainFolder_ = iConfig.getParameter<std::string>("mainFolder");
0006   hltPathNames_ = iConfig.getParameter<std::vector<std::string>>("HLTPathNames");
0007   edm::ParameterSet mc = iConfig.getParameter<edm::ParameterSet>("mcFlavours");
0008   m_mcLabels = mc.getParameterNamesForType<std::vector<unsigned int>>();
0009   m_histoName = iConfig.getParameter<std::vector<std::string>>("histoName");
0010   m_minTag = iConfig.getParameter<double>("minTag");
0011 
0012   HCALSpecialsNames[HEP17] = "HEP17";
0013   HCALSpecialsNames[HEP18] = "HEP18";
0014   HCALSpecialsNames[HEM17] = "HEM17";
0015 }
0016 
0017 HLTBTagHarvestingAnalyzer::~HLTBTagHarvestingAnalyzer() {}
0018 
0019 // ------------ method called once each job just after ending the event loop
0020 // ------------
0021 void HLTBTagHarvestingAnalyzer::dqmEndJob(DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter) {
0022   using namespace edm;
0023   Exception excp(errors::LogicError);
0024   std::string dqmFolder_hist;
0025 
0026   // for each hltPath and for each flavour, do the "b-tag efficiency vs jet pt"
0027   // and "b-tag efficiency vs mistag rate" plots
0028   for (unsigned int ind = 0; ind < hltPathNames_.size(); ind++) {
0029     dqmFolder_hist = Form("%s/Discriminator/%s", mainFolder_.c_str(), hltPathNames_[ind].c_str());
0030     std::string effDir = Form("%s/Discriminator/%s/efficiency", mainFolder_.c_str(), hltPathNames_[ind].c_str());
0031     std::string relationsDir = Form("%s/Discriminator/%s/HEP17_HEM17", mainFolder_.c_str(), hltPathNames_[ind].c_str());
0032     ibooker.setCurrentFolder(effDir);
0033     TH1 *den = nullptr;
0034     TH1 *num = nullptr;
0035     std::map<std::string, TH1F> effics;
0036     std::map<std::string, bool> efficsOK;
0037     std::map<std::string, std::map<HCALSpecials, TH1F>> efficsmod;
0038     std::map<std::string, std::map<HCALSpecials, bool>> efficsmodOK;
0039     for (unsigned int i = 0; i < m_mcLabels.size(); ++i) {
0040       bool isOK = false;
0041       std::string label = m_histoName.at(ind) + "__";  //"JetTag__";
0042       std::string flavour = m_mcLabels[i];
0043       label += flavour;
0044       isOK =
0045           GetNumDenumerators(ibooker, igetter, dqmFolder_hist + "/" + label, dqmFolder_hist + "/" + label, num, den, 0);
0046       if (isOK) {
0047         // do the 'b-tag efficiency vs discr' plot
0048         effics[flavour] = calculateEfficiency1D(ibooker, igetter, *num, *den, label + "_efficiency_vs_disc");
0049         efficsOK[flavour] = isOK;
0050       }
0051       // for modules (HEP17 etc.)
0052       for (const auto &j : HCALSpecialsNames) {
0053         ibooker.setCurrentFolder(dqmFolder_hist + "/" + j.second + "/efficiency");
0054         isOK = GetNumDenumerators(ibooker,
0055                                   igetter,
0056                                   dqmFolder_hist + "/" + j.second + "/" + label,
0057                                   dqmFolder_hist + "/" + j.second + "/" + label,
0058                                   num,
0059                                   den,
0060                                   0);
0061         if (isOK) {
0062           // do the 'b-tag efficiency vs discr' plot
0063           efficsmod[flavour][j.first] =
0064               calculateEfficiency1D(ibooker, igetter, *num, *den, label + "_efficiency_vs_disc");
0065           efficsmodOK[flavour][j.first] = isOK;
0066         }
0067       }
0068       ibooker.setCurrentFolder(effDir);
0069       label = m_histoName.at(ind) + "___";
0070       std::string labelEta = label;
0071       std::string labelPhi = label;
0072       label += flavour + "_disc_pT";
0073       labelEta += flavour + "_disc_eta";
0074       labelPhi += flavour + "_disc_phi";
0075       isOK =
0076           GetNumDenumerators(ibooker, igetter, dqmFolder_hist + "/" + label, dqmFolder_hist + "/" + label, num, den, 1);
0077       if (isOK) {
0078         // do the 'b-tag efficiency vs pT' plot
0079         TH1F eff = calculateEfficiency1D(ibooker, igetter, *num, *den, label + "_efficiency_vs_pT");
0080       }
0081       isOK = GetNumDenumerators(
0082           ibooker, igetter, dqmFolder_hist + "/" + labelEta, dqmFolder_hist + "/" + labelEta, num, den, 2);
0083       if (isOK) {
0084         // do the 'b-tag efficiency vs Eta' plot
0085         TH1F eff = calculateEfficiency1D(ibooker, igetter, *num, *den, labelEta + "_efficiency_vs_eta");
0086       }
0087       isOK = GetNumDenumerators(
0088           ibooker, igetter, dqmFolder_hist + "/" + labelPhi, dqmFolder_hist + "/" + labelPhi, num, den, 2);
0089       if (isOK) {
0090         // do the 'b-tag efficiency vs Phi' plot
0091         TH1F eff = calculateEfficiency1D(ibooker, igetter, *num, *den, labelPhi + "_efficiency_vs_phi");
0092       }
0093 
0094       /// save efficiency_vs_disc_HEP17 / efficiency_vs_disc_HEM17 plots
0095       ibooker.setCurrentFolder(relationsDir);
0096       if (efficsmodOK[flavour][HEP17] && efficsmodOK[flavour][HEM17])
0097         modulesrate(ibooker,
0098                     igetter,
0099                     &efficsmod[flavour][HEP17],
0100                     &efficsmod[flavour][HEM17],
0101                     m_histoName.at(ind) + "_" + flavour + "_HEP17_HEM17_effs_vs_disc_rate");
0102       ibooker.setCurrentFolder(effDir);
0103 
0104     }  /// for mc labels
0105 
0106     /// save mistagrate vs b-eff plots
0107     if (efficsOK["b"] && efficsOK["c"])
0108       mistagrate(ibooker, igetter, &effics["b"], &effics["c"], m_histoName.at(ind) + "_b_c_mistagrate");
0109     if (efficsOK["b"] && efficsOK["light"])
0110       mistagrate(ibooker, igetter, &effics["b"], &effics["light"], m_histoName.at(ind) + "_b_light_mistagrate");
0111     if (efficsOK["b"] && efficsOK["g"])
0112       mistagrate(ibooker, igetter, &effics["b"], &effics["g"], m_histoName.at(ind) + "_b_g_mistagrate");
0113 
0114     /// save mistagrate vs b-eff plots for modules (HEP17 etc.)
0115     for (const auto &j : HCALSpecialsNames) {
0116       ibooker.setCurrentFolder(dqmFolder_hist + "/" + j.second + "/efficiency");
0117       if (efficsmodOK["b"][j.first] && efficsmodOK["c"][j.first])
0118         mistagrate(ibooker,
0119                    igetter,
0120                    &efficsmod["b"][j.first],
0121                    &efficsmod["c"][j.first],
0122                    m_histoName.at(ind) + "_b_c_mistagrate");
0123       if (efficsmodOK["b"][j.first] && efficsmodOK["light"][j.first])
0124         mistagrate(ibooker,
0125                    igetter,
0126                    &efficsmod["b"][j.first],
0127                    &efficsmod["light"][j.first],
0128                    m_histoName.at(ind) + "_b_light_mistagrate");
0129       if (efficsmodOK["b"][j.first] && efficsmodOK["g"][j.first])
0130         mistagrate(ibooker,
0131                    igetter,
0132                    &efficsmod["b"][j.first],
0133                    &efficsmod["g"][j.first],
0134                    m_histoName.at(ind) + "_b_g_mistagrate");
0135     }
0136 
0137     /// save mistagrate_HEP17 / mistagrate_HEM17 plots
0138     ibooker.setCurrentFolder(relationsDir);
0139     bool isOK = false;
0140     isOK = GetNumDenumerators(ibooker,
0141                               igetter,
0142                               dqmFolder_hist + "/HEP17/efficiency/" + m_histoName.at(ind) + "_b_c_mistagrate",
0143                               dqmFolder_hist + "/HEM17/efficiency/" + m_histoName.at(ind) + "_b_c_mistagrate",
0144                               num,
0145                               den,
0146                               3);
0147     if (isOK)
0148       modulesrate(ibooker, igetter, (TH1F *)num, (TH1F *)den, m_histoName.at(ind) + "_HEP17_HEM17_b_c_mistagrate");
0149     isOK = GetNumDenumerators(ibooker,
0150                               igetter,
0151                               dqmFolder_hist + "/HEP17/efficiency/" + m_histoName.at(ind) + "_b_light_mistagrate",
0152                               dqmFolder_hist + "/HEM17/efficiency/" + m_histoName.at(ind) + "_b_light_mistagrate",
0153                               num,
0154                               den,
0155                               3);
0156     if (isOK)
0157       modulesrate(ibooker, igetter, (TH1F *)num, (TH1F *)den, m_histoName.at(ind) + "_HEP17_HEM17_b_light_mistagrate");
0158     isOK = GetNumDenumerators(ibooker,
0159                               igetter,
0160                               dqmFolder_hist + "/HEP17/efficiency/" + m_histoName.at(ind) + "_b_g_mistagrate",
0161                               dqmFolder_hist + "/HEM17/efficiency/" + m_histoName.at(ind) + "_b_g_mistagrate",
0162                               num,
0163                               den,
0164                               3);
0165     if (isOK)
0166       modulesrate(ibooker, igetter, (TH1F *)num, (TH1F *)den, m_histoName.at(ind) + "_HEP17_HEM17_b_g_mistagrate");
0167   }  /// for triggers
0168 }
0169 
0170 bool HLTBTagHarvestingAnalyzer::GetNumDenumerators(DQMStore::IBooker &ibooker,
0171                                                    DQMStore::IGetter &igetter,
0172                                                    std::string num,
0173                                                    std::string den,
0174                                                    TH1 *&ptrnum,
0175                                                    TH1 *&ptrden,
0176                                                    int type) {
0177   using namespace edm;
0178   /*
0179      possible types:
0180      type =0 for eff_vs_discriminator
0181      type =1 for eff_vs_pT
0182      type =2 for eff_vs_eta or eff_vs_phi
0183      type =3 for HEP17 / HEM17 mistagrate relation
0184    */
0185   MonitorElement *denME = nullptr;
0186   MonitorElement *numME = nullptr;
0187   denME = igetter.get(den);
0188   numME = igetter.get(num);
0189   Exception excp(errors::LogicError);
0190 
0191   if (denME == nullptr || numME == nullptr) {
0192     excp << "Plots not found:\n";
0193     if (denME == nullptr)
0194       excp << den << "\n";
0195     if (numME == nullptr)
0196       excp << num << "\n";
0197     excp.raise();
0198   }
0199 
0200   if (type == 0)  // efficiency_vs_discr: fill "ptrnum" with the cumulative function of
0201                   // the DQM plots contained in "num" and "ptrden" with a flat function
0202   {
0203     TH1 *numH1 = numME->getTH1();
0204     TH1 *denH1 = denME->getTH1();
0205     ptrden = (TH1 *)denH1->Clone("denominator");
0206     ptrnum = (TH1 *)numH1->Clone("numerator");
0207 
0208     ptrnum->SetBinContent(1, numH1->Integral());
0209     ptrden->SetBinContent(1, numH1->Integral());
0210     for (int j = 2; j <= numH1->GetNbinsX(); j++) {
0211       ptrnum->SetBinContent(j, numH1->Integral() - numH1->Integral(1, j - 1));
0212       ptrden->SetBinContent(j, numH1->Integral());
0213     }
0214   }
0215 
0216   if (type == 1)  // efficiency_vs_pT: fill "ptrden" with projection of the plots
0217                   // contained in "den" and fill "ptrnum" with projection of the
0218                   // plots contained in "num", having btag>m_minTag
0219   {
0220     TH2F *numH2 = numME->getTH2F();
0221     TH2F *denH2 = denME->getTH2F();
0222 
0223     /// numerator preparing
0224     TCutG *cutg_num = new TCutG("cutg_num", 4);
0225     cutg_num->SetPoint(0, m_minTag, 0);
0226     cutg_num->SetPoint(1, m_minTag, 9999);
0227     cutg_num->SetPoint(2, 1.1, 9999);
0228     cutg_num->SetPoint(3, 1.1, 0);
0229     ptrnum = numH2->ProjectionY("numerator", 0, -1, "[cutg_num]");
0230 
0231     /// denominator preparing
0232     TCutG *cutg_den = new TCutG("cutg_den", 4);
0233     cutg_den->SetPoint(0, -10.1, 0);
0234     cutg_den->SetPoint(1, -10.1, 9999);
0235     cutg_den->SetPoint(2, 1.1, 9999);
0236     cutg_den->SetPoint(3, 1.1, 0);
0237     ptrden = denH2->ProjectionY("denumerator", 0, -1, "[cutg_den]");
0238     delete cutg_num;
0239     delete cutg_den;
0240   }
0241 
0242   if (type == 2)  // efficiency_vs_eta: fill "ptrden" with projection of the
0243                   // plots contained in "den" and fill "ptrnum" with projection
0244                   // of the plots contained in "num", having btag>m_minTag
0245   {
0246     TH2F *numH2 = numME->getTH2F();
0247     TH2F *denH2 = denME->getTH2F();
0248 
0249     /// numerator preparing
0250     TCutG *cutg_num = new TCutG("cutg_num", 4);
0251     cutg_num->SetPoint(0, m_minTag, -10);
0252     cutg_num->SetPoint(1, m_minTag, 10);
0253     cutg_num->SetPoint(2, 1.1, 10);
0254     cutg_num->SetPoint(3, 1.1, -10);
0255     ptrnum = numH2->ProjectionY("numerator", 0, -1, "[cutg_num]");
0256 
0257     /// denominator preparing
0258     TCutG *cutg_den = new TCutG("cutg_den", 4);
0259     cutg_den->SetPoint(0, -10.1, -10);
0260     cutg_den->SetPoint(1, -10.1, 10);
0261     cutg_den->SetPoint(2, 1.1, 10);
0262     cutg_den->SetPoint(3, 1.1, -10);
0263     ptrden = denH2->ProjectionY("denumerator", 0, -1, "[cutg_den]");
0264     delete cutg_num;
0265     delete cutg_den;
0266   }
0267 
0268   if (type == 3)  // mistagrate HEP17 / HEM17 relation: fill "ptrnum" with HEP17
0269                   // mistagrate and "ptrden" with HEM17 mistagrate
0270   {
0271     ptrden = denME->getTH1();
0272     ptrnum = numME->getTH1();
0273   }
0274   return true;
0275 }
0276 
0277 void HLTBTagHarvestingAnalyzer::mistagrate(
0278     DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter, TH1F *num, TH1F *den, std::string effName) {
0279   // do the efficiency_vs_mistag_rate plot
0280   TH1F *eff;
0281   eff = new TH1F(effName.c_str(), effName.c_str(), 100, 0, 1);
0282   eff->SetTitle(effName.c_str());
0283   eff->SetXTitle("b-effficiency");
0284   eff->SetYTitle("mistag rate");
0285   eff->SetOption("E");
0286   eff->SetLineColor(2);
0287   eff->SetLineWidth(2);
0288   eff->SetMarkerStyle(20);
0289   eff->SetMarkerSize(0.8);
0290   eff->GetYaxis()->SetRangeUser(0.001, 1.001);
0291   eff->GetXaxis()->SetRangeUser(-0.001, 1.001);
0292   eff->SetStats(kFALSE);
0293 
0294   // for each bin in the discr -> find efficiency and mistagrate -> put them in
0295   // a plot
0296   for (int i = 1; i <= num->GetNbinsX(); i++) {
0297     double beff = num->GetBinContent(i);
0298     double miseff = den->GetBinContent(i);
0299     double miseffErr = den->GetBinError(i);
0300     int binX = eff->GetXaxis()->FindBin(beff);
0301     if (eff->GetBinContent(binX) != 0)
0302       continue;
0303     eff->SetBinContent(binX, miseff);
0304     eff->SetBinError(binX, miseffErr);
0305   }
0306   MonitorElement *me;
0307   me = ibooker.book1D(effName, eff);
0308   me->setEfficiencyFlag();
0309 
0310   delete eff;
0311   return;
0312 }
0313 
0314 void HLTBTagHarvestingAnalyzer::modulesrate(
0315     DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter, TH1F *num, TH1F *den, std::string effName) {
0316   // do the eff_vs_disc_HEP17 / eff_vs_disc_HEM17 plot
0317   TH1F *eff = new TH1F(*num);
0318   // eff = new TH1F(effName.c_str(),effName.c_str(),100,0,1);
0319   eff->Divide(den);
0320   eff->SetTitle(effName.c_str());
0321   eff->SetXTitle(num->GetXaxis()->GetTitle());
0322   eff->SetYTitle("");
0323   eff->SetOption("E");
0324   eff->SetLineColor(2);
0325   eff->SetLineWidth(2);
0326   eff->SetMarkerStyle(20);
0327   eff->SetMarkerSize(0.8);
0328   eff->GetYaxis()->SetRangeUser(0.001, 2.001);
0329   // eff->GetXaxis()->SetRangeUser(-0.001,1.001);
0330   eff->SetStats(kFALSE);
0331 
0332   MonitorElement *me;
0333   me = ibooker.book1D(effName, eff);
0334   me->setEfficiencyFlag();
0335 
0336   delete eff;
0337   return;
0338 }
0339 
0340 TH1F HLTBTagHarvestingAnalyzer::calculateEfficiency1D(
0341     DQMStore::IBooker &ibooker, DQMStore::IGetter &igetter, TH1 &num, TH1 &den, std::string effName) {
0342   // calculate the efficiency as num/den ratio
0343   TH1F eff;
0344   if (num.GetXaxis()->GetXbins()->GetSize() == 0) {
0345     eff = TH1F(effName.c_str(),
0346                effName.c_str(),
0347                num.GetXaxis()->GetNbins(),
0348                num.GetXaxis()->GetXmin(),
0349                num.GetXaxis()->GetXmax());
0350   } else {
0351     eff = TH1F(effName.c_str(), effName.c_str(), num.GetXaxis()->GetNbins(), num.GetXaxis()->GetXbins()->GetArray());
0352   }
0353   eff.SetTitle(effName.c_str());
0354   eff.SetXTitle(num.GetXaxis()->GetTitle());
0355   eff.SetYTitle("Efficiency");
0356   eff.SetOption("PE");
0357   eff.SetLineColor(2);
0358   eff.SetLineWidth(2);
0359   eff.SetMarkerStyle(20);
0360   eff.SetMarkerSize(0.8);
0361   eff.GetYaxis()->SetRangeUser(-0.001, 1.001);
0362   for (int i = 1; i <= num.GetNbinsX(); i++) {
0363     double d, n, err;
0364     d = den.GetBinContent(i);
0365     n = num.GetBinContent(i);
0366     double e;
0367     if (d != 0) {
0368       e = n / d;
0369       err = std::max(e - TEfficiency::ClopperPearson(d, n, 0.683, false),
0370                      TEfficiency::ClopperPearson(d, n, 0.683, true) - e);
0371       // err =  sqrt(e*(1-e)/d); //from binomial standard deviation
0372     } else {
0373       e = 0;
0374       err = 0;
0375     }
0376     eff.SetBinContent(i, e);
0377     eff.SetBinError(i, err);
0378   }
0379 
0380   MonitorElement *me;
0381   me = ibooker.book1D(effName, &eff);
0382   me->setEfficiencyFlag();
0383 
0384   return eff;
0385 }
0386 
0387 // define this as a plug-in
0388 DEFINE_FWK_MODULE(HLTBTagHarvestingAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!