RecoB/src/BTagDifferentialPlot.cc

0001 #include "DQMOffline/RecoB/interface/BTagDifferentialPlot.h"
0002 #include "DQMOffline/RecoB/interface/EffPurFromHistos.h"
0003 #include "DQMOffline/RecoB/interface/Tools.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "TF1.h"
0006
0007 #include "DQMServices/Core/interface/DQMStore.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009
0010 #include "DQMOffline/RecoB/interface/HistoProviderDQM.h"
0011
0012 #include <algorithm>
0013 #include <iostream>
0014 #include <sstream>
0015 using namespace std;
0016
0017 BTagDifferentialPlot::BTagDifferentialPlot(double bEff,
0018                                            const ConstVarType& constVariable,
0019                                            const std::string& tagName,
0020                                            unsigned int mc)
0021     : fixedBEfficiency(bEff),
0022       noProcessing(false),
0023       processed(false),
0024       constVar(constVariable),
0025       constVariableName(""),
0026       diffVariableName(""),
0027       constVariableValue(999.9, 999.9),
0028       commonName("MisidForBEff_" + tagName + "_"),
0029       theDifferentialHistoB_d(nullptr),
0030       theDifferentialHistoB_u(nullptr),
0031       theDifferentialHistoB_s(nullptr),
0032       theDifferentialHistoB_c(nullptr),
0033       theDifferentialHistoB_b(nullptr),
0034       theDifferentialHistoB_g(nullptr),
0035       theDifferentialHistoB_ni(nullptr),
0036       theDifferentialHistoB_dus(nullptr),
0037       theDifferentialHistoB_dusg(nullptr),
0038       theDifferentialHistoB_pu(nullptr),
0039       mcPlots_(mc) {}
0040
0041 BTagDifferentialPlot::~BTagDifferentialPlot() {}
0042
0043 void BTagDifferentialPlot::plot(TCanvas& thePlotCanvas) {
0044   //   thePlotCanvas = new TCanvas( commonName,
0045   //                commonName,
0046   //                btppXCanvas, btppYCanvas);
0047   //
0048   //   if (!btppTitle) gStyle->SetOptTitle(0);
0049
0050   if (!processed)
0051     return;
0052   //fixme:
0053   bool btppNI = false;
0054   bool btppColour = true;
0055
0056   thePlotCanvas.SetFillColor(0);
0057   thePlotCanvas.cd(1);
0058   gPad->SetLogy(1);
0059   gPad->SetGridx(1);
0060   gPad->SetGridy(1);
0061
0062   //  int col_b;
0063   int col_c;
0064   int col_g;
0065   int col_dus;
0066   int col_ni;
0067
0068   //  int mStyle_b;
0069   int mStyle_c;
0070   int mStyle_g;
0071   int mStyle_dus;
0072   int mStyle_ni;
0073
0074   // marker size(same for all)
0075   float mSize = 1.5;
0076
0077   if (btppColour) {
0078     //    col_b    = 2;
0079     col_c = 6;
0080     col_g = 3;
0081     col_dus = 4;
0082     col_ni = 5;
0083     //    mStyle_b   = 20;
0084     mStyle_c = 20;
0085     mStyle_g = 20;
0086     mStyle_dus = 20;
0087     mStyle_ni = 20;
0088   } else {
0089     //    col_b    = 1;
0090     col_c = 1;
0091     col_g = 1;
0092     col_dus = 1;
0093     col_ni = 1;
0094     //    mStyle_b   = 12;
0095     mStyle_c = 22;
0096     mStyle_g = 29;
0097     mStyle_dus = 20;
0098     mStyle_ni = 27;
0099   }
0100
0101   // for the moment: plot b(to see what the constant b-efficiency is), c, g, uds
0102   // b in red
0103   // No, do not plot b(because only visible for the soft leptons)
0104   // theDifferentialHistoB_b   -> GetXaxis()->SetTitle(diffVariableName);
0105   // theDifferentialHistoB_b   -> GetYaxis()->SetTitle("non b-jet efficiency");
0106   // theDifferentialHistoB_b   -> GetYaxis()->SetTitleOffset(1.25);
0107   // theDifferentialHistoB_b   -> SetMaximum(0.4);
0108   // theDifferentialHistoB_b   -> SetMinimum(1.e-4);
0109   // theDifferentialHistoB_b   -> SetMarkerColor(col_b);
0110   // theDifferentialHistoB_b   -> SetLineColor (col_b);
0111   // theDifferentialHistoB_b   -> SetMarkerSize(mSize);
0112   // theDifferentialHistoB_b   -> SetMarkerStyle(mStyle_b);
0113   // theDifferentialHistoB_b   -> SetStats(false);
0114   // theDifferentialHistoB_b   -> Draw("pe");
0115   // c in magenta
0116   theDifferentialHistoB_c->getTH1F()->SetMaximum(0.4);
0117   theDifferentialHistoB_c->getTH1F()->SetMinimum(1.e-4);
0118   theDifferentialHistoB_c->getTH1F()->SetMarkerColor(col_c);
0119   theDifferentialHistoB_c->getTH1F()->SetLineColor(col_c);
0120   theDifferentialHistoB_c->getTH1F()->SetMarkerSize(mSize);
0121   theDifferentialHistoB_c->getTH1F()->SetMarkerStyle(mStyle_c);
0122   theDifferentialHistoB_c->getTH1F()->SetStats(false);
0123   //  theDifferentialHistoB_c   -> Draw("peSame");
0124   theDifferentialHistoB_c->getTH1F()->Draw("pe");
0125   if (mcPlots_ > 2) {
0126     // uds in blue
0127     theDifferentialHistoB_dus->getTH1F()->SetMarkerColor(col_dus);
0128     theDifferentialHistoB_dus->getTH1F()->SetLineColor(col_dus);
0129     theDifferentialHistoB_dus->getTH1F()->SetMarkerSize(mSize);
0130     theDifferentialHistoB_dus->getTH1F()->SetMarkerStyle(mStyle_dus);
0131     theDifferentialHistoB_dus->getTH1F()->SetStats(false);
0132     theDifferentialHistoB_dus->getTH1F()->Draw("peSame");
0133     // g in green
0134     // only uds not to confuse
0135     theDifferentialHistoB_g->getTH1F()->SetMarkerColor(col_g);
0136     theDifferentialHistoB_g->getTH1F()->SetLineColor(col_g);
0137     theDifferentialHistoB_g->getTH1F()->SetMarkerSize(mSize);
0138     theDifferentialHistoB_g->getTH1F()->SetMarkerStyle(mStyle_g);
0139     theDifferentialHistoB_g->getTH1F()->SetStats(false);
0140     theDifferentialHistoB_g->getTH1F()->Draw("peSame");
0141   }
0142   // NI if wanted
0143   if (btppNI) {
0144     theDifferentialHistoB_ni->getTH1F()->SetMarkerColor(col_ni);
0145     theDifferentialHistoB_ni->getTH1F()->SetLineColor(col_ni);
0146     theDifferentialHistoB_ni->getTH1F()->SetMarkerSize(mSize);
0147     theDifferentialHistoB_ni->getTH1F()->SetMarkerStyle(mStyle_ni);
0148     theDifferentialHistoB_ni->getTH1F()->SetStats(false);
0149     theDifferentialHistoB_ni->getTH1F()->Draw("peSame");
0150   }
0151 }
0152
0153 void BTagDifferentialPlot::epsPlot(const std::string& name) { plot(name, ".eps"); }
0154
0155 void BTagDifferentialPlot::psPlot(const std::string& name) { plot(name, ".ps"); }
0156
0157 void BTagDifferentialPlot::plot(const std::string& name, const std::string& ext) {
0158   if (!processed)
0159     return;
0160   TCanvas tc(commonName.c_str(), commonName.c_str());
0161   plot(tc);
0162   tc.Print((name + commonName + ext).c_str());
0163 }
0164
0165 void BTagDifferentialPlot::process(DQMStore::IBooker& ibook) {
0166   setVariableName();  // also sets noProcessing if not OK
0167   if (noProcessing)
0168     return;
0169   bookHisto(ibook);
0170   fillHisto();
0171   processed = true;
0172 }
0173
0174 void BTagDifferentialPlot::setVariableName() {
0175   if (constVar == constETA) {
0176     constVariableName = "eta";
0177     diffVariableName = "pt";
0178     constVariableValue =
0179         make_pair(theBinPlotters[0]->etaPtBin().getEtaMin(), theBinPlotters[0]->etaPtBin().getEtaMax());
0180   }
0181   if (constVar == constPT) {
0182     constVariableName = "pt";
0183     diffVariableName = "eta";
0184     constVariableValue = make_pair(theBinPlotters[0]->etaPtBin().getPtMin(), theBinPlotters[0]->etaPtBin().getPtMax());
0185   }
0186 }
0187
0188 void BTagDifferentialPlot::bookHisto(DQMStore::IBooker& ibook) {
0189   // vector with ranges
0190   vector<float> variableRanges;
0191
0192   for (vector<std::shared_ptr<JetTagPlotter>>::const_iterator iP = theBinPlotters.begin(); iP != theBinPlotters.end();
0193        ++iP) {
0194     const EtaPtBin& currentBin = (*iP)->etaPtBin();
0195     if (diffVariableName == "eta") {
0196       // only active bins in the variable on x-axis
0197       if (currentBin.getEtaActive()) {
0198         variableRanges.push_back(currentBin.getEtaMin());
0199         // also max if last one
0200         if (iP == --theBinPlotters.end())
0201           variableRanges.push_back(currentBin.getEtaMax());
0202       }
0203     }
0204     if (diffVariableName == "pt") {
0205       // only active bins in the variable on x-axis
0206       if (currentBin.getPtActive()) {
0207         variableRanges.push_back(currentBin.getPtMin());
0208         // also max if last one
0209         if (iP == --theBinPlotters.end())
0210           variableRanges.push_back(currentBin.getPtMax());
0211       }
0212     }
0213   }
0214
0215   // to book histo with variable binning -> put into array
0216   int nBins = variableRanges.size() - 1;
0217   float* binArray = &variableRanges[0];
0218
0219   // part of the name common to all flavours
0220   std::stringstream stream("");
0221   stream << fixedBEfficiency << "_Const_" << constVariableName << "_" << constVariableValue.first << "-";
0222   stream << constVariableValue.second << "_"
0223          << "_Vs_" << diffVariableName;
0224   commonName += stream.str();
0225   std::remove(commonName.begin(), commonName.end(), ' ');
0226   std::replace(commonName.begin(), commonName.end(), '.', 'v');
0227
0228   std::string label(commonName);
0229   HistoProviderDQM prov("Btag", label, ibook);
0230
0231   theDifferentialHistoB_b = (prov.book1D("B_" + commonName, "B_" + commonName, nBins, binArray));
0232   theDifferentialHistoB_c = (prov.book1D("C_" + commonName, "C_" + commonName, nBins, binArray));
0233   theDifferentialHistoB_dusg = (prov.book1D("DUSG_" + commonName, "DUSG_" + commonName, nBins, binArray));
0234   theDifferentialHistoB_ni = (prov.book1D("NI_" + commonName, "NI_" + commonName, nBins, binArray));
0235   theDifferentialHistoB_pu = (prov.book1D("PU_" + commonName, "PU_" + commonName, nBins, binArray));
0236
0237   if (mcPlots_ > 2) {
0238     theDifferentialHistoB_d = (prov.book1D("D_" + commonName, "D_" + commonName, nBins, binArray));
0239     theDifferentialHistoB_u = (prov.book1D("U_" + commonName, "U_" + commonName, nBins, binArray));
0240     theDifferentialHistoB_s = (prov.book1D("S_" + commonName, "S_" + commonName, nBins, binArray));
0241     theDifferentialHistoB_g = (prov.book1D("G_" + commonName, "G_" + commonName, nBins, binArray));
0242     theDifferentialHistoB_dus = (prov.book1D("DUS_" + commonName, "DUS_" + commonName, nBins, binArray));
0243   }
0244 }
0245
0246 void BTagDifferentialPlot::fillHisto() {
0247   // loop over bins and find corresponding misid. in the MisIdVs..... histo
0248   for (vector<std::shared_ptr<JetTagPlotter>>::const_iterator iP = theBinPlotters.begin(); iP != theBinPlotters.end();
0249        ++iP) {
0250     const EtaPtBin& currentBin = (*iP)->etaPtBin();
0251     EffPurFromHistos& currentEffPurFromHistos = (*iP)->getEffPurFromHistos();
0252     //
0253     bool isActive = true;
0254     double valueXAxis = -999.99;
0255     // find right bin based on middle of the interval
0256     if (diffVariableName == "eta") {
0257       isActive = currentBin.getEtaActive();
0258       valueXAxis = 0.5 * (currentBin.getEtaMin() + currentBin.getEtaMax());
0259     } else if (diffVariableName == "pt") {
0260       isActive = currentBin.getPtActive();
0261       valueXAxis = 0.5 * (currentBin.getPtMin() + currentBin.getPtMax());
0262     } else {
0263       throw cms::Exception("Configuration")
0264           << "====>>>> BTagDifferentialPlot::fillHisto() : illegal diffVariableName = " << diffVariableName << endl;
0265     }
0266
0267     // for the moment: ignore inactive bins
0268     //(maybe later: if a Bin is inactive -> set value to fill well below left edge of histogram to have it in the underflow)
0269
0270     if (!isActive)
0271       continue;
0272
0273     // to have less lines of code ....
0274     vector<pair<TH1F*, TH1F*>> effPurDifferentialPairs;
0275
0276     // all flavours(b is a good cross check! must be constant and = fixed b-efficiency)
0277     // get histo; find the bin of the fixed b-efficiency in the histo and get misid; fill
0278
0279     effPurDifferentialPairs.push_back(
0280         make_pair(currentEffPurFromHistos.getEffFlavVsBEff_b(), theDifferentialHistoB_b->getTH1F()));
0281     effPurDifferentialPairs.push_back(
0282         make_pair(currentEffPurFromHistos.getEffFlavVsBEff_c(), theDifferentialHistoB_c->getTH1F()));
0283     effPurDifferentialPairs.push_back(
0284         make_pair(currentEffPurFromHistos.getEffFlavVsBEff_dusg(), theDifferentialHistoB_dusg->getTH1F()));
0285     effPurDifferentialPairs.push_back(
0286         make_pair(currentEffPurFromHistos.getEffFlavVsBEff_ni(), theDifferentialHistoB_ni->getTH1F()));
0287     effPurDifferentialPairs.push_back(
0288         make_pair(currentEffPurFromHistos.getEffFlavVsBEff_pu(), theDifferentialHistoB_pu->getTH1F()));
0289     if (mcPlots_ > 2) {
0290       effPurDifferentialPairs.push_back(
0291           make_pair(currentEffPurFromHistos.getEffFlavVsBEff_d(), theDifferentialHistoB_d->getTH1F()));
0292       effPurDifferentialPairs.push_back(
0293           make_pair(currentEffPurFromHistos.getEffFlavVsBEff_u(), theDifferentialHistoB_u->getTH1F()));
0294       effPurDifferentialPairs.push_back(
0295           make_pair(currentEffPurFromHistos.getEffFlavVsBEff_s(), theDifferentialHistoB_s->getTH1F()));
0296       effPurDifferentialPairs.push_back(
0297           make_pair(currentEffPurFromHistos.getEffFlavVsBEff_g(), theDifferentialHistoB_g->getTH1F()));
0298       effPurDifferentialPairs.push_back(
0299           make_pair(currentEffPurFromHistos.getEffFlavVsBEff_dus(), theDifferentialHistoB_dus->getTH1F()));
0300     }
0301
0302     for (vector<pair<TH1F*, TH1F*>>::const_iterator itP = effPurDifferentialPairs.begin();
0303          itP != effPurDifferentialPairs.end();
0304          ++itP) {
0305       TH1F* effPurHist = itP->first;
0306       TH1F* diffHist = itP->second;
0307       pair<double, double> mistag = getMistag(fixedBEfficiency, effPurHist);
0308       int iBinSet = diffHist->FindBin(valueXAxis);
0309       diffHist->SetBinContent(iBinSet, mistag.first);
0310       diffHist->SetBinError(iBinSet, mistag.second);
0311     }
0312   }
0313 }
0314
0315 pair<double, double> BTagDifferentialPlot::getMistag(const double& fixedBEfficiency, TH1F* effPurHist) {
0316   int iBinGet = effPurHist->FindBin(fixedBEfficiency);
0317   double effForBEff = effPurHist->GetBinContent(iBinGet);
0318   double effForBEffErr = effPurHist->GetBinError(iBinGet);
0319
0320   if (effForBEff == 0. && effForBEffErr == 0.) {
0321     // The bin was empty. Could be that it was not filled, as the scan-plot
0322     //  did not have an entry at the requested value, or that the curve
0323     // is above or below.
0324     // Fit a plynomial, and evaluate the mistag at the requested value.
0325     int fitStatus;
0326     //need our own copy for thread safety
0327     TF1 myfunc("myfunc", "pol4");
0328     try {
0329       fitStatus = effPurHist->Fit(&myfunc, "q");
0330     } catch (cms::Exception& iException) {
0331       return pair<double, double>(effForBEff, effForBEffErr);
0332     }
0333     if (fitStatus != 0) {
0334       edm::LogWarning("BTagDifferentialPlot") << "Fit failed to hisogram " << effPurHist->GetTitle()
0335                                               << ", perhaps because too few entries = " << effPurHist->GetEntries()
0336                                               << ". This bin will be missing in plots at fixed b efficiency.";
0337       //    } else {
0338       //      edm::LogInfo("BTagDifferentialPlot")<<"Fit OK to hisogram " << effPurHist->GetTitle() << " entries = " << effPurHist->GetEntries();
0339       return pair<double, double>(effForBEff, effForBEffErr);
0340     }
0341     effForBEff = myfunc.Eval(fixedBEfficiency);
0342     effPurHist->RecursiveRemove(&myfunc);
0343     //Error: first non-empty bin on the right and take its error
0344     for (int i = iBinGet + 1; i < effPurHist->GetNbinsX(); ++i) {
0345       if (effPurHist->GetBinContent(i) != 0) {
0346         effForBEffErr = effPurHist->GetBinError(i);
0347         break;
0348       }
0349     }
0350   }
0351
0352   return pair<double, double>(effForBEff, effForBEffErr);
0353 }