Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // consistency check
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   // conditional discriminator cut for efficiency histos
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   // discr. for computation
0158   vector<TH2F*> discrCfHistos = dDiscriminatorFC.getHistoVector();
0159   // discr no cut
0160   vector<TH2F*> discrNoCutHistos = discrNoCutEffic->getHistoVector();
0161   // discr no cut
0162   vector<TH2F*> discrCutHistos = discrCutEfficScan->getHistoVector();
0163 
0164   const int& dimHistos = discrCfHistos.size();  // they all have the same size
0165 
0166   // DISCR-CUT LOOP:
0167   // fill the histos for eff-pur computations by scanning the discriminatorFC histogram
0168 
0169   // better to loop over bins -> discrCut no longer needed
0170   const int& nBinsX = dDiscriminatorFC.nBinsX();
0171   const int& nBinsY = dDiscriminatorFC.nBinsY();
0172 
0173   // loop over flavours
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     // In Root histos, bin counting starts at 1 to nBins.
0183     // bin 0 is the underflow, and nBins+1 is the overflow.
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         // fill all jets into NoCut histo
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   // divide to get efficiency vs. discriminator cut from absolute numbers
0199   discrCutEfficScan->divide(*discrNoCutEffic);  // does: histos including discriminator cut / flat histo
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 /* = 0 */) {
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   // number of bins
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   // start
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   // end
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   // to have shorter names ......
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   // create histograms from base name and extension as given from user
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   // loop over eff. vs. discriminator cut b-histo and look in which bin the closest entry is;
0380   // use fact that eff decreases monotonously
0381 
0382   // any of the histos to be created can be taken here:
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++) {  // loop over the bins on the x-axis of the histograms to be filled
0388     const float& effBinWidthX = EffFlavVsXEff->getTH1F()->GetBinWidth(iBinX);
0389     const float& effMidX = EffFlavVsXEff->getTH1F()->GetBinCenter(iBinX);  // middle of efficiency bin
0390     const float& effLeftX = effMidX - 0.5 * effBinWidthX;                  // left edge of bin
0391     const float& effRightX = effMidX + 0.5 * effBinWidthX;                 // right edge of bin
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         // fill the histos
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>
[ 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!