File indexing completed on 2024-04-06 12:09:48
0001 #include "DQMOffline/RecoB/interface/EffPurFromHistos2D.h"
0002 #include "DQMOffline/RecoB/interface/Tools.h"
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "TStyle.h"
0005 #include "TCanvas.h"
0006
0007 #include <iostream>
0008 #include <cmath>
0009
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012
0013 using namespace std;
0014 using namespace RecoBTag;
0015
0016 EffPurFromHistos2D::EffPurFromHistos2D(const std::string& ext,
0017 TH2F* h_d,
0018 TH2F* h_u,
0019 TH2F* h_s,
0020 TH2F* h_c,
0021 TH2F* h_b,
0022 TH2F* h_g,
0023 TH2F* h_ni,
0024 TH2F* h_dus,
0025 TH2F* h_dusg,
0026 TH2F* h_pu,
0027 const std::string& label,
0028 unsigned int mc,
0029 int nBinX,
0030 double startOX,
0031 double endOX)
0032 : fromDiscriminatorDistr(false),
0033 mcPlots_(mc),
0034 doCTagPlots_(false),
0035 label_(label),
0036 histoExtension(ext),
0037 effVersusDiscr_d(h_d),
0038 effVersusDiscr_u(h_u),
0039 effVersusDiscr_s(h_s),
0040 effVersusDiscr_c(h_c),
0041 effVersusDiscr_b(h_b),
0042 effVersusDiscr_g(h_g),
0043 effVersusDiscr_ni(h_ni),
0044 effVersusDiscr_dus(h_dus),
0045 effVersusDiscr_dusg(h_dusg),
0046 effVersusDiscr_pu(h_pu),
0047 nBinOutputX(nBinX),
0048 startOutputX(startOX),
0049 endOutputX(endOX) {
0050
0051 check();
0052 }
0053
0054 EffPurFromHistos2D::EffPurFromHistos2D(const FlavourHistograms2D<double, double>& dDiscriminatorFC,
0055 const std::string& label,
0056 unsigned int mc,
0057 DQMStore::IBooker& ibook,
0058 int nBinX,
0059 double startOX,
0060 double endOX)
0061 : fromDiscriminatorDistr(true),
0062 mcPlots_(mc),
0063 doCTagPlots_(false),
0064 label_(label),
0065 nBinOutputX(nBinX),
0066 startOutputX(startOX),
0067 endOutputX(endOX) {
0068 histoExtension = "_" + dDiscriminatorFC.baseNameTitle();
0069
0070 discrNoCutEffic =
0071 std::make_unique<FlavourHistograms2D<double, double>>("totalEntries" + histoExtension,
0072 "Total Entries: " + dDiscriminatorFC.baseNameDescription(),
0073 dDiscriminatorFC.nBinsX(),
0074 dDiscriminatorFC.lowerBoundX(),
0075 dDiscriminatorFC.upperBoundX(),
0076 dDiscriminatorFC.nBinsY(),
0077 dDiscriminatorFC.lowerBoundY(),
0078 dDiscriminatorFC.upperBoundY(),
0079 false,
0080 label,
0081 mcPlots_,
0082 false,
0083 ibook);
0084
0085
0086 discrCutEfficScan = std::make_unique<FlavourHistograms2D<double, double>>(
0087 "effVsDiscrCut" + histoExtension,
0088 "Eff. vs Disc. Cut: " + dDiscriminatorFC.baseNameDescription(),
0089 dDiscriminatorFC.nBinsX(),
0090 dDiscriminatorFC.lowerBoundX(),
0091 dDiscriminatorFC.upperBoundX(),
0092 dDiscriminatorFC.nBinsY(),
0093 dDiscriminatorFC.lowerBoundY(),
0094 dDiscriminatorFC.upperBoundY(),
0095 false,
0096 label,
0097 mcPlots_,
0098 false,
0099 ibook);
0100 discrCutEfficScan->SetMinimum(1E-4);
0101
0102 if (mcPlots_) {
0103 if (mcPlots_ > 2) {
0104 effVersusDiscr_d = discrCutEfficScan->histo_d();
0105 effVersusDiscr_u = discrCutEfficScan->histo_u();
0106 effVersusDiscr_s = discrCutEfficScan->histo_s();
0107 effVersusDiscr_g = discrCutEfficScan->histo_g();
0108 effVersusDiscr_dus = discrCutEfficScan->histo_dus();
0109 } else {
0110 effVersusDiscr_d = nullptr;
0111 effVersusDiscr_u = nullptr;
0112 effVersusDiscr_s = nullptr;
0113 effVersusDiscr_g = nullptr;
0114 effVersusDiscr_dus = nullptr;
0115 }
0116 effVersusDiscr_c = discrCutEfficScan->histo_c();
0117 effVersusDiscr_b = discrCutEfficScan->histo_b();
0118 effVersusDiscr_ni = discrCutEfficScan->histo_ni();
0119 effVersusDiscr_dusg = discrCutEfficScan->histo_dusg();
0120 effVersusDiscr_pu = discrCutEfficScan->histo_pu();
0121
0122 if (mcPlots_ > 2) {
0123 effVersusDiscr_d->SetXTitle("Discriminant");
0124 effVersusDiscr_d->GetXaxis()->SetTitleOffset(0.75);
0125 effVersusDiscr_u->SetXTitle("Discriminant");
0126 effVersusDiscr_u->GetXaxis()->SetTitleOffset(0.75);
0127 effVersusDiscr_s->SetXTitle("Discriminant");
0128 effVersusDiscr_s->GetXaxis()->SetTitleOffset(0.75);
0129 effVersusDiscr_g->SetXTitle("Discriminant");
0130 effVersusDiscr_g->GetXaxis()->SetTitleOffset(0.75);
0131 effVersusDiscr_dus->SetXTitle("Discriminant");
0132 effVersusDiscr_dus->GetXaxis()->SetTitleOffset(0.75);
0133 }
0134 effVersusDiscr_c->SetXTitle("Discriminant");
0135 effVersusDiscr_c->GetXaxis()->SetTitleOffset(0.75);
0136 effVersusDiscr_b->SetXTitle("Discriminant");
0137 effVersusDiscr_b->GetXaxis()->SetTitleOffset(0.75);
0138 effVersusDiscr_ni->SetXTitle("Discriminant");
0139 effVersusDiscr_ni->GetXaxis()->SetTitleOffset(0.75);
0140 effVersusDiscr_dusg->SetXTitle("Discriminant");
0141 effVersusDiscr_dusg->GetXaxis()->SetTitleOffset(0.75);
0142 effVersusDiscr_pu->SetXTitle("Discriminant");
0143 effVersusDiscr_pu->GetXaxis()->SetTitleOffset(0.75);
0144 } else {
0145 effVersusDiscr_d = nullptr;
0146 effVersusDiscr_u = nullptr;
0147 effVersusDiscr_s = nullptr;
0148 effVersusDiscr_c = nullptr;
0149 effVersusDiscr_b = nullptr;
0150 effVersusDiscr_g = nullptr;
0151 effVersusDiscr_ni = nullptr;
0152 effVersusDiscr_dus = nullptr;
0153 effVersusDiscr_dusg = nullptr;
0154 effVersusDiscr_pu = nullptr;
0155 }
0156
0157
0158 vector<TH2F*> discrCfHistos = dDiscriminatorFC.getHistoVector();
0159
0160 vector<TH2F*> discrNoCutHistos = discrNoCutEffic->getHistoVector();
0161
0162 vector<TH2F*> discrCutHistos = discrCutEfficScan->getHistoVector();
0163
0164 const int& dimHistos = discrCfHistos.size();
0165
0166
0167
0168
0169
0170 const int& nBinsX = dDiscriminatorFC.nBinsX();
0171 const int& nBinsY = dDiscriminatorFC.nBinsY();
0172
0173
0174 for (int iFlav = 0; iFlav < dimHistos; iFlav++) {
0175 if (discrCfHistos[iFlav] == nullptr)
0176 continue;
0177 discrNoCutHistos[iFlav]->SetXTitle("Discriminant A");
0178 discrNoCutHistos[iFlav]->GetXaxis()->SetTitleOffset(0.75);
0179 discrNoCutHistos[iFlav]->SetYTitle("Discriminant B");
0180 discrNoCutHistos[iFlav]->GetYaxis()->SetTitleOffset(0.75);
0181
0182
0183
0184 const double& nJetsFlav = discrCfHistos[iFlav]->GetEntries();
0185
0186 for (int iDiscrX = nBinsX; iDiscrX > 0; --iDiscrX) {
0187 for (int iDiscrY = nBinsY; iDiscrY > 0; --iDiscrY) {
0188
0189 discrNoCutHistos[iFlav]->SetBinContent(iDiscrX, iDiscrY, nJetsFlav);
0190 discrNoCutHistos[iFlav]->SetBinError(iDiscrX, iDiscrY, sqrt(nJetsFlav));
0191 const double& sum = nJetsFlav - discrCfHistos[iFlav]->Integral(0, iDiscrX - 1, 0, iDiscrY - 1);
0192 discrCutHistos[iFlav]->SetBinContent(iDiscrX, iDiscrY, sum);
0193 discrCutHistos[iFlav]->SetBinError(iDiscrX, iDiscrY, sqrt(sum));
0194 }
0195 }
0196 }
0197
0198
0199 discrCutEfficScan->divide(*discrNoCutEffic);
0200 discrCutEfficScan->setEfficiencyFlag();
0201 }
0202
0203 EffPurFromHistos2D::~EffPurFromHistos2D() {}
0204
0205 void EffPurFromHistos2D::epsPlot(const std::string& name) {
0206 if (fromDiscriminatorDistr) {
0207 discrNoCutEffic->epsPlot(name);
0208 discrCutEfficScan->epsPlot(name);
0209 }
0210 plot(name, ".eps");
0211 }
0212
0213 void EffPurFromHistos2D::psPlot(const std::string& name) { plot(name, ".ps"); }
0214
0215 void EffPurFromHistos2D::plot(const std::string& name, const std::string& ext) {
0216 std::string hX = "";
0217 std::string Title = "";
0218 if (!doCTagPlots_) {
0219 hX = "FlavEffVsBEff";
0220 Title = "b";
0221 } else {
0222 hX = "FlavEffVsCEff";
0223 Title = "c";
0224 }
0225 TCanvas tc((hX + histoExtension).c_str(),
0226 ("Flavour misidentification vs. " + Title + "-tagging efficiency " + histoExtension).c_str());
0227 plot(&tc);
0228 tc.Print((name + hX + histoExtension + ext).c_str());
0229 }
0230
0231 void EffPurFromHistos2D::plot(TPad* plotCanvas ) {
0232 setTDRStyle()->cd();
0233
0234 if (plotCanvas)
0235 plotCanvas->cd();
0236
0237 gPad->UseCurrentStyle();
0238 gPad->SetFillColor(0);
0239 gPad->SetLogy(1);
0240 gPad->SetGridx(1);
0241 gPad->SetGridy(1);
0242 }
0243
0244 void EffPurFromHistos2D::check() {
0245
0246 int nBinsX_d = 0;
0247 int nBinsX_u = 0;
0248 int nBinsX_s = 0;
0249 int nBinsX_g = 0;
0250 int nBinsX_dus = 0;
0251 int nBinsY_d = 0;
0252 int nBinsY_u = 0;
0253 int nBinsY_s = 0;
0254 int nBinsY_g = 0;
0255 int nBinsY_dus = 0;
0256 if (mcPlots_ > 2) {
0257 nBinsX_d = effVersusDiscr_d->GetNbinsX();
0258 nBinsX_u = effVersusDiscr_u->GetNbinsX();
0259 nBinsX_s = effVersusDiscr_s->GetNbinsX();
0260 nBinsX_g = effVersusDiscr_g->GetNbinsX();
0261 nBinsX_dus = effVersusDiscr_dus->GetNbinsX();
0262 nBinsY_d = effVersusDiscr_d->GetNbinsY();
0263 nBinsY_u = effVersusDiscr_u->GetNbinsY();
0264 nBinsY_s = effVersusDiscr_s->GetNbinsY();
0265 nBinsY_g = effVersusDiscr_g->GetNbinsY();
0266 nBinsY_dus = effVersusDiscr_dus->GetNbinsY();
0267 }
0268 const int& nBinsX_c = effVersusDiscr_c->GetNbinsX();
0269 const int& nBinsX_b = effVersusDiscr_b->GetNbinsX();
0270 const int& nBinsX_ni = effVersusDiscr_ni->GetNbinsX();
0271 const int& nBinsX_dusg = effVersusDiscr_dusg->GetNbinsX();
0272 const int& nBinsX_pu = effVersusDiscr_pu->GetNbinsX();
0273 const int& nBinsY_c = effVersusDiscr_c->GetNbinsY();
0274 const int& nBinsY_b = effVersusDiscr_b->GetNbinsY();
0275 const int& nBinsY_ni = effVersusDiscr_ni->GetNbinsY();
0276 const int& nBinsY_dusg = effVersusDiscr_dusg->GetNbinsY();
0277 const int& nBinsY_pu = effVersusDiscr_pu->GetNbinsY();
0278
0279 const bool& lNBinsX =
0280 ((nBinsX_d == nBinsX_u && nBinsX_d == nBinsX_s && nBinsX_d == nBinsX_c && nBinsX_d == nBinsX_b &&
0281 nBinsX_d == nBinsX_g && nBinsX_d == nBinsX_ni && nBinsX_d == nBinsX_dus && nBinsX_d == nBinsX_dusg) ||
0282 (nBinsX_c == nBinsX_b && nBinsX_c == nBinsX_dusg && nBinsX_c == nBinsX_ni && nBinsX_c == nBinsX_pu));
0283
0284 const bool& lNBinsY =
0285 ((nBinsY_d == nBinsY_u && nBinsY_d == nBinsY_s && nBinsY_d == nBinsY_c && nBinsY_d == nBinsY_b &&
0286 nBinsY_d == nBinsY_g && nBinsY_d == nBinsY_ni && nBinsY_d == nBinsY_dus && nBinsY_d == nBinsY_dusg) ||
0287 (nBinsY_c == nBinsY_b && nBinsY_c == nBinsY_dusg && nBinsY_c == nBinsY_ni && nBinsY_c == nBinsY_pu));
0288
0289 if (!lNBinsX || !lNBinsY) {
0290 throw cms::Exception("Configuration") << "Input histograms do not all have the same number of bins!\n";
0291 }
0292
0293
0294 float sBin_d = 0;
0295 float sBin_u = 0;
0296 float sBin_s = 0;
0297 float sBin_g = 0;
0298 float sBin_dus = 0;
0299 if (mcPlots_ > 2) {
0300 sBin_d = effVersusDiscr_d->GetBinCenter(1);
0301 sBin_u = effVersusDiscr_u->GetBinCenter(1);
0302 sBin_s = effVersusDiscr_s->GetBinCenter(1);
0303 sBin_g = effVersusDiscr_g->GetBinCenter(1);
0304 sBin_dus = effVersusDiscr_dus->GetBinCenter(1);
0305 }
0306 const float& sBin_c = effVersusDiscr_c->GetBinCenter(1);
0307 const float& sBin_b = effVersusDiscr_b->GetBinCenter(1);
0308 const float& sBin_ni = effVersusDiscr_ni->GetBinCenter(1);
0309 const float& sBin_dusg = effVersusDiscr_dusg->GetBinCenter(1);
0310 const float& sBin_pu = effVersusDiscr_pu->GetBinCenter(1);
0311
0312 const bool& lSBin = ((sBin_d == sBin_u && sBin_d == sBin_s && sBin_d == sBin_c && sBin_d == sBin_b &&
0313 sBin_d == sBin_g && sBin_d == sBin_ni && sBin_d == sBin_dus && sBin_d == sBin_dusg) ||
0314 (sBin_c == sBin_b && sBin_c == sBin_dusg && sBin_c == sBin_ni && sBin_c == sBin_pu));
0315
0316 if (!lSBin) {
0317 throw cms::Exception("Configuration")
0318 << "EffPurFromHistos::check() : Input histograms do not all have the same start bin!\n";
0319 }
0320
0321
0322 float eBin_d = 0;
0323 float eBin_u = 0;
0324 float eBin_s = 0;
0325 float eBin_g = 0;
0326 float eBin_dus = 0;
0327 const int& binEnd = effVersusDiscr_b->GetBin(nBinsX_b - 1, nBinsY_b - 1);
0328 if (mcPlots_ > 2) {
0329 eBin_d = effVersusDiscr_d->GetBinCenter(binEnd);
0330 eBin_u = effVersusDiscr_u->GetBinCenter(binEnd);
0331 eBin_s = effVersusDiscr_s->GetBinCenter(binEnd);
0332 eBin_g = effVersusDiscr_g->GetBinCenter(binEnd);
0333 eBin_dus = effVersusDiscr_dus->GetBinCenter(binEnd);
0334 }
0335 const float& eBin_c = effVersusDiscr_c->GetBinCenter(binEnd);
0336 const float& eBin_b = effVersusDiscr_b->GetBinCenter(binEnd);
0337 const float& eBin_ni = effVersusDiscr_ni->GetBinCenter(binEnd);
0338 const float& eBin_dusg = effVersusDiscr_dusg->GetBinCenter(binEnd);
0339 const float& eBin_pu = effVersusDiscr_pu->GetBinCenter(binEnd);
0340
0341 const bool& lEBin = ((eBin_d == eBin_u && eBin_d == eBin_s && eBin_d == eBin_c && eBin_d == eBin_b &&
0342 eBin_d == eBin_g && eBin_d == eBin_ni && eBin_d == eBin_dus && eBin_d == eBin_dusg) ||
0343 (eBin_c == eBin_b && eBin_c == eBin_dusg && eBin_c == eBin_ni && eBin_c == eBin_pu));
0344
0345 if (!lEBin) {
0346 throw cms::Exception("Configuration")
0347 << "EffPurFromHistos::check() : Input histograms do not all have the same end bin!\n";
0348 }
0349 }
0350
0351 void EffPurFromHistos2D::compute(DQMStore::IBooker& ibook, vector<double> fixedEff) {
0352 if (!mcPlots_ || fixedEff.empty()) {
0353 return;
0354 }
0355
0356
0357 const std::string& hE = histoExtension;
0358 std::string hX = "DUSG_vs_B_eff_at_fixedCeff_";
0359 if (!doCTagPlots_)
0360 hX = "DUSG_vs_C_eff_at_fixedBeff_";
0361
0362
0363 HistoProviderDQM prov("Btag", label_, ibook);
0364
0365 for (unsigned int ieff = 0; ieff < fixedEff.size(); ieff++) {
0366 std::string fixedEfficiency = std::to_string(fixedEff[ieff]);
0367 fixedEfficiency.replace(1, 1, "_");
0368 X_vs_Y_eff_at_fixedZeff.push_back(
0369 (prov.book1D(hX + fixedEfficiency + hE, hX + fixedEfficiency + hE, nBinOutputX, startOutputX, endOutputX)));
0370 X_vs_Y_eff_at_fixedZeff[ieff]->setEfficiencyFlag();
0371
0372 X_vs_Y_eff_at_fixedZeff[ieff]->setAxisTitle("Light mistag");
0373 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetYTitle("B mistag");
0374 if (!doCTagPlots_)
0375 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetYTitle("C mistag");
0376 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->GetXaxis()->SetTitleOffset(0.75);
0377 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->GetYaxis()->SetTitleOffset(0.75);
0378 }
0379
0380
0381
0382
0383 MonitorElement* EffFlavVsXEff = X_vs_Y_eff_at_fixedZeff[0];
0384
0385 const int& nBinX = EffFlavVsXEff->getTH1F()->GetNbinsX();
0386
0387 for (int iBinX = 1; iBinX <= nBinX; iBinX++) {
0388 const float& effBinWidthX = EffFlavVsXEff->getTH1F()->GetBinWidth(iBinX);
0389 const float& effMidX = EffFlavVsXEff->getTH1F()->GetBinCenter(iBinX);
0390 const float& effLeftX = effMidX - 0.5 * effBinWidthX;
0391 const float& effRightX = effMidX + 0.5 * effBinWidthX;
0392
0393 vector<int> binClosest;
0394 if (doCTagPlots_)
0395 binClosest =
0396 findBinClosestYValueAtFixedZ(effVersusDiscr_dusg, effMidX, effLeftX, effRightX, effVersusDiscr_c, fixedEff);
0397 else
0398 binClosest =
0399 findBinClosestYValueAtFixedZ(effVersusDiscr_dusg, effMidX, effLeftX, effRightX, effVersusDiscr_b, fixedEff);
0400
0401 for (unsigned int ieff = 0; ieff < binClosest.size(); ieff++) {
0402 const bool& binFound = (binClosest[ieff] > 0);
0403
0404 if (binFound) {
0405
0406 if (doCTagPlots_) {
0407 X_vs_Y_eff_at_fixedZeff[ieff]->Fill(effMidX, effVersusDiscr_b->GetBinContent(binClosest[ieff]));
0408 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetBinError(iBinX, effVersusDiscr_b->GetBinError(binClosest[ieff]));
0409 } else {
0410 X_vs_Y_eff_at_fixedZeff[ieff]->Fill(effMidX, effVersusDiscr_c->GetBinContent(binClosest[ieff]));
0411 X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetBinError(iBinX, effVersusDiscr_c->GetBinError(binClosest[ieff]));
0412 }
0413 }
0414 }
0415 }
0416 }
0417
0418 #include <typeinfo>