RecoB/interface/FlavourHistorgrams.h

0001
0002 #ifndef FlavourHistograms_H
0003 #define FlavourHistograms_H
0004
0005 #include "DQMServices/Core/interface/DQMStore.h"
0006 #include "FWCore/ServiceRegistry/interface/Service.h"
0007
0008 // #include "BTagPlotPrintC.h"
0009
0010 #include "TH1F.h"
0011 #include "TCanvas.h"
0012 #include "TROOT.h"
0013 #include "TSystem.h"
0014 #include "TStyle.h"
0015 #include "TEfficiency.h"
0016
0017 #include "DQMOffline/RecoB/interface/Tools.h"
0018 #include "DQMOffline/RecoB/interface/HistoProviderDQM.h"
0019 #include <iostream>
0020 #include <string>
0021
0022 //
0023 // class to describe Histo
0024 //
0025 template <class T>
0026 class FlavourHistograms {
0027 public:
0028   typedef dqm::legacy::DQMStore DQMStore;
0029   typedef dqm::legacy::MonitorElement MonitorElement;
0030
0031   FlavourHistograms(const std::string& baseNameTitle_,
0032                     const std::string& baseNameDescription_,
0033                     const int& nBins_,
0034                     const double& lowerBound_,
0035                     const double& upperBound_,
0036                     const std::string& plotFirst_,
0037                     const std::string& folder,
0038                     const unsigned int& mc,
0039                     DQMStore::IGetter& iget);
0040
0041   FlavourHistograms(const std::string& baseNameTitle_,
0042                     const std::string& baseNameDescription_,
0043                     const int& nBins_,
0044                     const double& lowerBound_,
0045                     const double& upperBound_,
0046                     const bool& statistics_,
0047                     const bool& plotLog_,
0048                     const bool& plotNormalized_,
0049                     const std::string& plotFirst_,
0050                     const std::string& folder,
0051                     const unsigned int& mc,
0052                     DQMStore::IBooker& ibook);
0053
0054   virtual ~FlavourHistograms();
0055
0056   // fill entry
0057   // For single variables and arrays (for arrays only a single index can be filled)
0058   void fill(const int& flavour, const T& variable) const;
0059   void fill(const int& flavour, const T& variable, const T& w) const;
0060
0061   // For single variables and arrays
0062   void fill(const int& flavour, const T* variable) const;
0063
0064   void settitle(const char* title);
0065
0066   void plot(TPad* theCanvas = nullptr);
0067
0068   void epsPlot(const std::string& name);
0069
0070   void divide(const FlavourHistograms<T>& bHD);
0071   void setEfficiencyFlag();
0072
0073   inline void SetMaximum(const double& max) { theMax = max; }
0074   inline void SetMinimum(const double& min) { theMin = min; }
0075
0076   // trivial access functions
0077   inline std::string baseNameTitle() const { return theBaseNameTitle; }
0078   inline std::string baseNameDescription() const { return theBaseNameDescription; }
0079   inline int nBins() const { return theNBins; }
0080   inline double lowerBound() const { return theLowerBound; }
0081   inline double upperBound() const { return theUpperBound; }
0082   inline bool statistics() const { return theStatistics; }
0083   inline bool plotLog() const { return thePlotLog; }
0084   inline bool plotNormalized() const { return thePlotNormalized; }
0085   inline std::string plotFirst() const { return thePlotFirst; }
0086   inline int* arrayDimension() const { return theArrayDimension; }
0087   inline int maxDimension() const { return theMaxDimension; }
0088   inline int indexToPlot() const { return theIndexToPlot; }
0089
0090   // access to the histos
0091   inline TH1F* histo_all() const { return theHisto_all->getTH1F(); }
0092   inline TH1F* histo_d() const { return theHisto_d->getTH1F(); }
0093   inline TH1F* histo_u() const { return theHisto_u->getTH1F(); }
0094   inline TH1F* histo_s() const { return theHisto_s->getTH1F(); }
0095   inline TH1F* histo_c() const { return theHisto_c->getTH1F(); }
0096   inline TH1F* histo_b() const { return theHisto_b->getTH1F(); }
0097   inline TH1F* histo_g() const { return theHisto_g->getTH1F(); }
0098   inline TH1F* histo_ni() const { return theHisto_ni->getTH1F(); }
0099   inline TH1F* histo_dus() const { return theHisto_dus->getTH1F(); }
0100   inline TH1F* histo_dusg() const { return theHisto_dusg->getTH1F(); }
0101   inline TH1F* histo_pu() const { return theHisto_pu->getTH1F(); }
0102
0103   std::vector<TH1F*> getHistoVector() const;
0104
0105 protected:
0106   void fillVariable(const int& flavour, const T& var, const T& w) const;
0107   double ClopperPearsonUnc(double num, double den);
0108   void ComputeEfficiency(TH1F* num, TH1F* den, int bin);
0109
0110   //
0111   // the data members
0112   //
0113
0114   //   T *   theVariable   ;
0115
0116   // for arrays
0117   int* theArrayDimension;
0118   int theMaxDimension;
0119   int theIndexToPlot;  // in case that not the complete array has to be plotted
0120
0121   std::string theBaseNameTitle;
0122   std::string theBaseNameDescription;
0123   int theNBins;
0124   double theLowerBound;
0125   double theUpperBound;
0126   bool theStatistics;
0127   bool thePlotLog;
0128   bool thePlotNormalized;
0129   std::string thePlotFirst;  // one character; gives flavour to plot first: l (light) , c , b
0130   double theMin, theMax;
0131
0132   // the histos
0133   MonitorElement* theHisto_all;
0134   MonitorElement* theHisto_d;
0135   MonitorElement* theHisto_u;
0136   MonitorElement* theHisto_s;
0137   MonitorElement* theHisto_c;
0138   MonitorElement* theHisto_b;
0139   MonitorElement* theHisto_g;
0140   MonitorElement* theHisto_ni;
0141   MonitorElement* theHisto_dus;
0142   MonitorElement* theHisto_dusg;
0143   MonitorElement* theHisto_pu;
0144
0145   //  DQMStore * dqmStore_;
0146
0147   // the canvas to plot
0148   TCanvas* theCanvas;
0149
0150 private:
0151   FlavourHistograms() {}
0152
0153   unsigned int mcPlots_;
0154 };
0155
0156 template <class T>
0157 FlavourHistograms<T>::FlavourHistograms(const std::string& baseNameTitle_,
0158                                         const std::string& baseNameDescription_,
0159                                         const int& nBins_,
0160                                         const double& lowerBound_,
0161                                         const double& upperBound_,
0162                                         const std::string& plotFirst_,
0163                                         const std::string& folder,
0164                                         const unsigned int& mc,
0165                                         DQMStore::IGetter& iget)
0166     : theMaxDimension(-1),
0167       theIndexToPlot(-1),
0168       theBaseNameTitle(baseNameTitle_),
0169       theBaseNameDescription(baseNameDescription_),
0170       theNBins(nBins_),
0171       theLowerBound(lowerBound_),
0172       theUpperBound(upperBound_),
0173       theStatistics(false),
0174       thePlotLog(false),
0175       thePlotNormalized(false),
0176       thePlotFirst(plotFirst_),
0177       theMin(-1.),
0178       theMax(-1.),
0179       mcPlots_(mc) {
0180   // defaults for array dimensions
0181   theArrayDimension = nullptr;
0182
0183   // check plot order string
0184   if (thePlotFirst == "l" || thePlotFirst == "c" || thePlotFirst == "b") {
0185     // OK
0186   } else {
0187     // not correct: print warning and set default (l)
0188     std::cout << "FlavourHistograms::FlavourHistograms : thePlotFirst was not correct : " << thePlotFirst << std::endl;
0189     std::cout << "FlavourHistograms::FlavourHistograms : Set it to default value (l)! " << std::endl;
0190     thePlotFirst = "l";
0191   }
0192
0193   if (mcPlots_ % 2 == 0)
0194     theHisto_all = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "ALL");
0195   else
0196     theHisto_all = nullptr;
0197   if (mcPlots_) {
0198     if (mcPlots_ > 2) {
0199       theHisto_d = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "D");
0200       theHisto_u = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "U");
0201       theHisto_s = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "S");
0202       theHisto_g = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "G");
0203       theHisto_dus = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "DUS");
0204     } else {
0205       theHisto_d = nullptr;
0206       theHisto_u = nullptr;
0207       theHisto_s = nullptr;
0208       theHisto_g = nullptr;
0209       theHisto_dus = nullptr;
0210     }
0211     theHisto_c = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "C");
0212     theHisto_b = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "B");
0213     theHisto_ni = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "NI");
0214     theHisto_dusg = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "DUSG");
0215     theHisto_pu = iget.get("Btag/" + folder + "/" + theBaseNameTitle + "PU");
0216   } else {
0217     theHisto_d = nullptr;
0218     theHisto_u = nullptr;
0219     theHisto_s = nullptr;
0220     theHisto_c = nullptr;
0221     theHisto_b = nullptr;
0222     theHisto_g = nullptr;
0223     theHisto_ni = nullptr;
0224     theHisto_dus = nullptr;
0225     theHisto_dusg = nullptr;
0226     theHisto_pu = nullptr;
0227   }
0228   // defaults for other data members
0229   theCanvas = nullptr;
0230 }
0231
0232 template <class T>
0233 FlavourHistograms<T>::FlavourHistograms(const std::string& baseNameTitle_,
0234                                         const std::string& baseNameDescription_,
0235                                         const int& nBins_,
0236                                         const double& lowerBound_,
0237                                         const double& upperBound_,
0238                                         const bool& statistics_,
0239                                         const bool& plotLog_,
0240                                         const bool& plotNormalized_,
0241                                         const std::string& plotFirst_,
0242                                         const std::string& folder,
0243                                         const unsigned int& mc,
0244                                         DQMStore::IBooker& ibook)
0245     : theMaxDimension(-1),
0246       theIndexToPlot(-1),
0247       theBaseNameTitle(baseNameTitle_),
0248       theBaseNameDescription(baseNameDescription_),
0249       theNBins(nBins_),
0250       theLowerBound(lowerBound_),
0251       theUpperBound(upperBound_),
0252       theStatistics(statistics_),
0253       thePlotLog(plotLog_),
0254       thePlotNormalized(plotNormalized_),
0255       thePlotFirst(plotFirst_),
0256       theMin(-1.),
0257       theMax(-1.),
0258       mcPlots_(mc) {
0259   // defaults for array dimensions
0260   theArrayDimension = nullptr;
0261
0262   // check plot order string
0263   if (thePlotFirst == "l" || thePlotFirst == "c" || thePlotFirst == "b") {
0264     // OK
0265   } else {
0266     // not correct: print warning and set default (l)
0267     std::cout << "FlavourHistograms::FlavourHistograms : thePlotFirst was not correct : " << thePlotFirst << std::endl;
0268     std::cout << "FlavourHistograms::FlavourHistograms : Set it to default value (l)! " << std::endl;
0269     thePlotFirst = "l";
0270   }
0271
0272   // book histos
0273   HistoProviderDQM prov("Btag", folder, ibook);
0274   if (mcPlots_ % 2 == 0)
0275     theHisto_all = (prov.book1D(
0276         theBaseNameTitle + "ALL", theBaseNameDescription + " all jets", theNBins, theLowerBound, theUpperBound));
0277   else
0278     theHisto_all = nullptr;
0279   if (mcPlots_) {
0280     if (mcPlots_ > 2) {
0281       theHisto_d = (prov.book1D(
0282           theBaseNameTitle + "D", theBaseNameDescription + " d-jets", theNBins, theLowerBound, theUpperBound));
0283       theHisto_u = (prov.book1D(
0284           theBaseNameTitle + "U", theBaseNameDescription + " u-jets", theNBins, theLowerBound, theUpperBound));
0285       theHisto_s = (prov.book1D(
0286           theBaseNameTitle + "S", theBaseNameDescription + " s-jets", theNBins, theLowerBound, theUpperBound));
0287       theHisto_g = (prov.book1D(
0288           theBaseNameTitle + "G", theBaseNameDescription + " g-jets", theNBins, theLowerBound, theUpperBound));
0289       theHisto_dus = (prov.book1D(
0290           theBaseNameTitle + "DUS", theBaseNameDescription + " dus-jets", theNBins, theLowerBound, theUpperBound));
0291     } else {
0292       theHisto_d = nullptr;
0293       theHisto_u = nullptr;
0294       theHisto_s = nullptr;
0295       theHisto_g = nullptr;
0296       theHisto_dus = nullptr;
0297     }
0298     theHisto_c = (prov.book1D(
0299         theBaseNameTitle + "C", theBaseNameDescription + " c-jets", theNBins, theLowerBound, theUpperBound));
0300     theHisto_b = (prov.book1D(
0301         theBaseNameTitle + "B", theBaseNameDescription + " b-jets", theNBins, theLowerBound, theUpperBound));
0302     theHisto_ni = (prov.book1D(
0303         theBaseNameTitle + "NI", theBaseNameDescription + " ni-jets", theNBins, theLowerBound, theUpperBound));
0304     theHisto_dusg = (prov.book1D(
0305         theBaseNameTitle + "DUSG", theBaseNameDescription + " dusg-jets", theNBins, theLowerBound, theUpperBound));
0306     theHisto_pu = (prov.book1D(
0307         theBaseNameTitle + "PU", theBaseNameDescription + " pu-jets", theNBins, theLowerBound, theUpperBound));
0308   } else {
0309     theHisto_d = nullptr;
0310     theHisto_u = nullptr;
0311     theHisto_s = nullptr;
0312     theHisto_c = nullptr;
0313     theHisto_b = nullptr;
0314     theHisto_g = nullptr;
0315     theHisto_ni = nullptr;
0316     theHisto_dus = nullptr;
0317     theHisto_dusg = nullptr;
0318     theHisto_pu = nullptr;
0319   }
0320
0321   // statistics if requested
0322   if (theStatistics) {
0323     if (theHisto_all)
0324       theHisto_all->enableSumw2();
0325     if (mcPlots_) {
0326       if (mcPlots_ > 2) {
0327         theHisto_d->enableSumw2();
0328         theHisto_u->enableSumw2();
0329         theHisto_s->enableSumw2();
0330         theHisto_g->enableSumw2();
0331         theHisto_dus->enableSumw2();
0332       }
0333       theHisto_c->enableSumw2();
0334       theHisto_b->enableSumw2();
0335       theHisto_ni->enableSumw2();
0336       theHisto_dusg->enableSumw2();
0337       theHisto_pu->enableSumw2();
0338     }
0339   }
0340   // defaults for other data members
0341   theCanvas = nullptr;
0342 }
0343
0344 template <class T>
0345 FlavourHistograms<T>::~FlavourHistograms() {
0346   // delete the canvas*/
0347   delete theCanvas;
0348 }
0349
0350 // fill entry
0351 template <class T>
0352 void FlavourHistograms<T>::fill(const int& flavour, const T& variable) const {
0353   // For single variables and arrays (for arrays only a single index can be filled)
0354   fillVariable(flavour, variable, 1.);
0355 }
0356
0357 template <class T>
0358 void FlavourHistograms<T>::fill(const int& flavour, const T& variable, const T& w) const {
0359   // For single variables and arrays (for arrays only a single index can be filled)
0360   fillVariable(flavour, variable, w);
0361 }
0362
0363 template <class T>
0364 void FlavourHistograms<T>::fill(const int& flavour, const T* variable) const {
0365   if (theArrayDimension == nullptr) {
0366     // single variable
0367     fillVariable(flavour, *variable, 1.);
0368   } else {
0369     // array
0370     int iMax = (*theArrayDimension > theMaxDimension) ? theMaxDimension : *theArrayDimension;
0371     //
0372     for (int i = 0; i != iMax; ++i) {
0373       // check if only one index to be plotted (<0: switched off -> plot all)
0374       if ((theIndexToPlot < 0) || (i == theIndexToPlot)) {
0375         fillVariable(flavour, *(variable + i), 1.);
0376       }
0377     }
0378
0379     // if single index to be filled but not enough entries: fill 0.0 (convention!)
0380     if (theIndexToPlot >= iMax) {
0381       // cout << "==>> The index to be filled is too big -> fill 0.0 : " << theBaseNameTitle << " : " << theIndexToPlot << " >= " << iMax << endl ;
0382       const T& theZero = static_cast<T>(0.0);
0383       fillVariable(flavour, theZero, 1.);
0384     }
0385   }
0386 }
0387
0388 template <class T>
0389 void FlavourHistograms<T>::settitle(const char* title) {
0390   if (theHisto_all)
0391     theHisto_all->setAxisTitle(title);
0392   if (mcPlots_) {
0393     if (mcPlots_ > 2) {
0394       theHisto_d->setAxisTitle(title);
0395       theHisto_u->setAxisTitle(title);
0396       theHisto_s->setAxisTitle(title);
0397       theHisto_g->setAxisTitle(title);
0398       theHisto_dus->setAxisTitle(title);
0399     }
0400     theHisto_c->setAxisTitle(title);
0401     theHisto_b->setAxisTitle(title);
0402     theHisto_ni->setAxisTitle(title);
0403     theHisto_dusg->setAxisTitle(title);
0404     theHisto_pu->setAxisTitle(title);
0405   }
0406 }
0407
0408 template <class T>
0409 void FlavourHistograms<T>::plot(TPad* theCanvas /* = 0 */) {
0410   //fixme:
0411   bool btppNI = false;
0412   bool btppColour = true;
0413
0414   if (theCanvas)
0415     theCanvas->cd();
0416
0417   RecoBTag::setTDRStyle()->cd();
0418   gPad->UseCurrentStyle();
0419   //   if ( !btppTitle ) gStyle->SetOptTitle ( 0 ) ;
0420   //
0421   //   // here: plot histograms in a canvas
0422   //   theCanvas = new TCanvas ( "C" + theBaseNameTitle , "C" + theBaseNameDescription ,
0423   //                 btppXCanvas , btppYCanvas ) ;
0424   //   theCanvas->SetFillColor ( 0 ) ;
0425   //   theCanvas->cd  ( 1 ) ;
0426   gPad->SetLogy(0);
0427   if (thePlotLog)
0428     gPad->SetLogy(1);
0429   gPad->SetGridx(0);
0430   gPad->SetGridy(0);
0431   gPad->SetTitle(nullptr);
0432
0433   MonitorElement* histo[4];
0434   int col[4], lineStyle[4], markerStyle[4];
0435   int lineWidth = 1;
0436
0437   const double markerSize = gPad->GetWh() * gPad->GetHNDC() / 500.;
0438
0439   // default (l)
0440   histo[0] = theHisto_dusg;
0441   //CW histo_1 = theHisto_dus ;
0442   histo[1] = theHisto_b;
0443   histo[2] = theHisto_c;
0444   histo[3] = nullptr;
0445
0446   double max = theMax;
0447   if (theMax <= 0.) {
0448     max = theHisto_dusg->getTH1F()->GetMaximum();
0449     if (theHisto_b->getTH1F()->GetMaximum() > max)
0450       max = theHisto_b->getTH1F()->GetMaximum();
0451     if (theHisto_c->getTH1F()->GetMaximum() > max)
0452       max = theHisto_c->getTH1F()->GetMaximum();
0453   }
0454
0455   if (btppNI) {
0456     histo[3] = theHisto_ni;
0457     if (theHisto_ni->getTH1F()->GetMaximum() > max)
0458       max = theHisto_ni->getTH1F()->GetMaximum();
0459   }
0460
0461   if (btppColour) {  // print colours
0462     col[0] = 4;
0463     col[1] = 2;
0464     col[2] = 6;
0465     col[3] = 3;
0466     lineStyle[0] = 1;
0467     lineStyle[1] = 1;
0468     lineStyle[2] = 1;
0469     lineStyle[3] = 1;
0470     markerStyle[0] = 20;
0471     markerStyle[1] = 21;
0472     markerStyle[2] = 22;
0473     markerStyle[3] = 23;
0474     lineWidth = 1;
0475   } else {  // different marker/line styles
0476     col[1] = 1;
0477     col[2] = 1;
0478     col[3] = 1;
0479     col[0] = 1;
0480     lineStyle[0] = 2;
0481     lineStyle[1] = 1;
0482     lineStyle[2] = 3;
0483     lineStyle[3] = 4;
0484     markerStyle[0] = 20;
0485     markerStyle[1] = 21;
0486     markerStyle[2] = 22;
0487     markerStyle[3] = 23;
0488   }
0489
0490   // if changing order (NI stays always last)
0491
0492   // c to plot first
0493   if (thePlotFirst == "c") {
0494     histo[0] = theHisto_c;
0495     if (btppColour)
0496       col[0] = 6;
0497     if (!btppColour)
0498       lineStyle[0] = 3;
0499     histo[2] = theHisto_dusg;
0500     if (btppColour)
0501       col[2] = 4;
0502     if (!btppColour)
0503       lineStyle[2] = 2;
0504   }
0505
0506   // b to plot first
0507   if (thePlotFirst == "b") {
0508     histo[0] = theHisto_b;
0509     if (btppColour)
0510       col[0] = 2;
0511     if (!btppColour)
0512       lineStyle[0] = 1;
0513     histo[1] = theHisto_dusg;
0514     if (btppColour)
0515       col[1] = 4;
0516     if (!btppColour)
0517       lineStyle[1] = 2;
0518   }
0519
0520   histo[0]->setAxisTitle(theBaseNameDescription);
0521   histo[0]->getTH1F()->GetYaxis()->SetTitle("Arbitrary Units");
0522   histo[0]->getTH1F()->GetYaxis()->SetTitleOffset(1.25);
0523
0524   for (int i = 0; i != 4; ++i) {
0525     if (histo[i] == nullptr)
0526       continue;
0527     histo[i]->getTH1F()->SetStats(false);
0528     histo[i]->getTH1F()->SetLineStyle(lineStyle[i]);
0529     histo[i]->getTH1F()->SetLineWidth(lineWidth);
0530     histo[i]->getTH1F()->SetLineColor(col[i]);
0531     histo[i]->getTH1F()->SetMarkerStyle(markerStyle[i]);
0532     histo[i]->getTH1F()->SetMarkerColor(col[i]);
0533     histo[i]->getTH1F()->SetMarkerSize(markerSize);
0534   }
0535
0536   if (thePlotNormalized) {
0537     if (histo[0]->getTH1F()->GetEntries() != 0) {
0538       histo[0]->getTH1F()->DrawNormalized();
0539     } else {
0540       histo[0]->getTH1F()->SetMaximum(1.0);
0541       histo[0]->getTH1F()->Draw();
0542     }
0543     if (histo[1]->getTH1F()->GetEntries() != 0)
0544       histo[1]->getTH1F()->DrawNormalized("Same");
0545     if (histo[2]->getTH1F()->GetEntries() != 0)
0546       histo[2]->getTH1F()->DrawNormalized("Same");
0547     if ((histo[3] != nullptr) && (histo[3]->getTH1F()->GetEntries() != 0))
0548       histo[3]->getTH1F()->DrawNormalized("Same");
0549   } else {
0550     histo[0]->getTH1F()->SetMaximum(max * 1.05);
0551     if (theMin != -1.)
0552       histo[0]->getTH1F()->SetMinimum(theMin);
0553     histo[0]->getTH1F()->Draw();
0554     histo[1]->getTH1F()->Draw("Same");
0555     histo[2]->getTH1F()->Draw("Same");
0556     if (histo[3] != nullptr)
0557       histo[3]->getTH1F()->Draw("Same");
0558   }
0559 }
0560
0561 template <class T>
0562 void FlavourHistograms<T>::epsPlot(const std::string& name) {
0563   TCanvas tc(theBaseNameTitle.c_str(), theBaseNameDescription.c_str());
0564
0565   plot(&tc);
0566   tc.Print((name + theBaseNameTitle + ".eps").c_str());
0567 }
0568
0569 template <class T>
0570 double FlavourHistograms<T>::ClopperPearsonUnc(double num, double den) {
0571   double effVal = num / den;
0572   double errLo = TEfficiency::ClopperPearson(static_cast<int>(den), static_cast<int>(num), 0.683, false);
0573   double errUp = TEfficiency::ClopperPearson(static_cast<int>(den), static_cast<int>(num), 0.683, true);
0574   return std::max(effVal - errLo, errUp - effVal);
0575 }
0576
0577 template <class T>
0578 void FlavourHistograms<T>::ComputeEfficiency(TH1F* num, TH1F* den, int bin) {
0579   double effVal = 1.;
0580   double errVal = 0.;
0581   double numVal = num->GetBinContent(bin);
0582   double denVal = den->GetBinContent(bin);
0583   if (denVal > 0 && numVal <= denVal) {
0584     effVal = numVal / denVal;
0585     errVal = ClopperPearsonUnc(numVal, denVal);
0586   }
0587   num->SetBinContent(bin, effVal);
0588   num->SetBinError(bin, errVal);
0589 }
0590
0591 template <class T>
0592 void FlavourHistograms<T>::divide(const FlavourHistograms<T>& bHD) {
0593   for (int bin = 0; bin < theNBins + 2; bin++) {
0594     if (theHisto_all)
0595       ComputeEfficiency(theHisto_all->getTH1F(), bHD.histo_all(), bin);
0596     if (mcPlots_) {
0597       if (mcPlots_ > 2) {
0598         ComputeEfficiency(theHisto_d->getTH1F(), bHD.histo_d(), bin);
0599         ComputeEfficiency(theHisto_u->getTH1F(), bHD.histo_u(), bin);
0600         ComputeEfficiency(theHisto_s->getTH1F(), bHD.histo_s(), bin);
0601         ComputeEfficiency(theHisto_g->getTH1F(), bHD.histo_g(), bin);
0602         ComputeEfficiency(theHisto_dus->getTH1F(), bHD.histo_dus(), bin);
0603       }
0604       ComputeEfficiency(theHisto_c->getTH1F(), bHD.histo_c(), bin);
0605       ComputeEfficiency(theHisto_b->getTH1F(), bHD.histo_b(), bin);
0606       ComputeEfficiency(theHisto_ni->getTH1F(), bHD.histo_ni(), bin);
0607       ComputeEfficiency(theHisto_dusg->getTH1F(), bHD.histo_dusg(), bin);
0608       ComputeEfficiency(theHisto_pu->getTH1F(), bHD.histo_pu(), bin);
0609     }
0610   }
0611 }
0612
0613 template <class T>
0614 void FlavourHistograms<T>::setEfficiencyFlag() {
0615   if (theHisto_all)
0616     theHisto_all->setEfficiencyFlag();
0617   if (mcPlots_) {
0618     if (mcPlots_ > 2) {
0619       theHisto_d->setEfficiencyFlag();
0620       theHisto_u->setEfficiencyFlag();
0621       theHisto_s->setEfficiencyFlag();
0622       theHisto_g->setEfficiencyFlag();
0623       theHisto_dus->setEfficiencyFlag();
0624     }
0625     theHisto_c->setEfficiencyFlag();
0626     theHisto_b->setEfficiencyFlag();
0627     theHisto_ni->setEfficiencyFlag();
0628     theHisto_dusg->setEfficiencyFlag();
0629     theHisto_pu->setEfficiencyFlag();
0630   }
0631 }
0632
0633 template <class T>
0634 void FlavourHistograms<T>::fillVariable(const int& flavour, const T& var, const T& w) const {
0635   // all, except for the Jet Multiplicity which is not filled for each jets but for each events
0636   if ((theBaseNameDescription != "Jet Multiplicity" || flavour == -1) && theHisto_all)
0637     theHisto_all->Fill(var, w);
0638
0639   // flavour specific
0640   if (!mcPlots_ || (theBaseNameDescription == "Jet Multiplicity" && flavour == -1))
0641     return;
0642
0643   switch (std::abs(flavour)) {
0644     case 1:
0645       if (mcPlots_ > 2) {
0646         theHisto_d->Fill(var, w);
0647         if (theBaseNameDescription != "Jet Multiplicity")
0648           theHisto_dus->Fill(var, w);
0649       }
0650       if (theBaseNameDescription != "Jet Multiplicity")
0651         theHisto_dusg->Fill(var, w);
0652       return;
0653     case 2:
0654       if (mcPlots_ > 2) {
0655         theHisto_u->Fill(var, w);
0656         if (theBaseNameDescription != "Jet Multiplicity")
0657           theHisto_dus->Fill(var, w);
0658       }
0659       if (theBaseNameDescription != "Jet Multiplicity")
0660         theHisto_dusg->Fill(var, w);
0661       return;
0662     case 3:
0663       if (mcPlots_ > 2) {
0664         theHisto_s->Fill(var, w);
0665         if (theBaseNameDescription != "Jet Multiplicity")
0666           theHisto_dus->Fill(var, w);
0667       }
0668       if (theBaseNameDescription != "Jet Multiplicity")
0669         theHisto_dusg->Fill(var, w);
0670       return;
0671     case 4:
0672       theHisto_c->Fill(var, w);
0673       return;
0674     case 5:
0675       theHisto_b->Fill(var, w);
0676       return;
0677     case 21:
0678       if (mcPlots_ > 2)
0679         theHisto_g->Fill(var, w);
0680       if (theBaseNameDescription != "Jet Multiplicity")
0681         theHisto_dusg->Fill(var, w);
0682       return;
0683     case 123:
0684       if (mcPlots_ > 2 && theBaseNameDescription == "Jet Multiplicity")
0685         theHisto_dus->Fill(var, w);
0686       return;
0687     case 12321:
0688       if (theBaseNameDescription == "Jet Multiplicity")
0689         theHisto_dusg->Fill(var, w);
0690       return;
0691     case 20:
0692       theHisto_pu->Fill(var, w);
0693       return;
0694     default:
0695       theHisto_ni->Fill(var, w);
0696       return;
0697   }
0698 }
0699
0700 template <class T>
0701 std::vector<TH1F*> FlavourHistograms<T>::getHistoVector() const {
0702   std::vector<TH1F*> histoVector;
0703   if (theHisto_all)
0704     histoVector.push_back(theHisto_all->getTH1F());
0705   if (mcPlots_) {
0706     if (mcPlots_ > 2) {
0707       histoVector.push_back(theHisto_d->getTH1F());
0708       histoVector.push_back(theHisto_u->getTH1F());
0709       histoVector.push_back(theHisto_s->getTH1F());
0710       histoVector.push_back(theHisto_g->getTH1F());
0711       histoVector.push_back(theHisto_dus->getTH1F());
0712     }
0713     histoVector.push_back(theHisto_c->getTH1F());
0714     histoVector.push_back(theHisto_b->getTH1F());
0715     histoVector.push_back(theHisto_ni->getTH1F());
0716     histoVector.push_back(theHisto_dusg->getTH1F());
0717     histoVector.push_back(theHisto_pu->getTH1F());
0718   }
0719   return histoVector;
0720 }
0721 #endif