CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EffPurFromHistos2D.cc
Go to the documentation of this file.
4 #include "TStyle.h"
5 #include "TCanvas.h"
6 
7 #include <iostream>
8 #include <math.h>
9 
13 
14 using namespace std;
15 using namespace RecoBTag;
16 
17 EffPurFromHistos2D::EffPurFromHistos2D ( const std::string & ext, TH2F * h_d, TH2F * h_u,
18  TH2F * h_s, TH2F * h_c, TH2F * h_b, TH2F * h_g, TH2F * h_ni,
19  TH2F * h_dus, TH2F * h_dusg, TH2F * h_pu,
20  const std::string& label, const unsigned int& mc,
21  int nBinX, double startOX, double endOX) :
22  fromDiscriminatorDistr(false),
23  mcPlots_(mc), doCTagPlots_(false), label_(label),
24  histoExtension(ext), effVersusDiscr_d(h_d), effVersusDiscr_u(h_u),
25  effVersusDiscr_s(h_s), effVersusDiscr_c(h_c), effVersusDiscr_b(h_b),
26  effVersusDiscr_g(h_g), effVersusDiscr_ni(h_ni), effVersusDiscr_dus(h_dus),
27  effVersusDiscr_dusg(h_dusg), effVersusDiscr_pu(h_pu),
28  nBinOutputX(nBinX), startOutputX(startOX), endOutputX(endOX)
29 {
30  // consistency check
31  check();
32 }
33 
35  const std::string& label, const unsigned int& mc,
36  DQMStore::IBooker & ibook,
37  int nBinX, double startOX, double endOX) :
38  fromDiscriminatorDistr(true),
39  mcPlots_(mc), doCTagPlots_(false), label_(label),
40  nBinOutputX(nBinX), startOutputX(startOX), endOutputX(endOX)
41 {
42  histoExtension = "_"+dDiscriminatorFC->baseNameTitle();
43 
45  "totalEntries" + histoExtension, "Total Entries: " + dDiscriminatorFC->baseNameDescription(),
46  dDiscriminatorFC->nBinsX(), dDiscriminatorFC->lowerBoundX(), dDiscriminatorFC->upperBoundX(), dDiscriminatorFC->nBinsY(),
47  dDiscriminatorFC->lowerBoundY(), dDiscriminatorFC->upperBoundY(),
48  false, label, mcPlots_, false, ibook));
49  // conditional discriminator cut for efficiency histos
50 
52  "effVsDiscrCut" + histoExtension, "Eff. vs Disc. Cut: " + dDiscriminatorFC->baseNameDescription(),
53  dDiscriminatorFC->nBinsX(), dDiscriminatorFC->lowerBoundX(), dDiscriminatorFC->upperBoundX(), dDiscriminatorFC->nBinsY(),
54  dDiscriminatorFC->lowerBoundY(), dDiscriminatorFC->upperBoundY(),
55  false, label, mcPlots_, false, ibook));
56  discrCutEfficScan->SetMinimum(1E-4);
57 
58  if (mcPlots_){
59 
60  if(mcPlots_>2){
61  effVersusDiscr_d = discrCutEfficScan->histo_d ();
62  effVersusDiscr_u = discrCutEfficScan->histo_u ();
63  effVersusDiscr_s = discrCutEfficScan->histo_s ();
64  effVersusDiscr_g = discrCutEfficScan->histo_g ();
65  effVersusDiscr_dus = discrCutEfficScan->histo_dus ();
66  }
67  else{
68  effVersusDiscr_d = 0;
69  effVersusDiscr_u = 0;
70  effVersusDiscr_s = 0;
71  effVersusDiscr_g = 0;
73  }
74  effVersusDiscr_c = discrCutEfficScan->histo_c ();
75  effVersusDiscr_b = discrCutEfficScan->histo_b ();
76  effVersusDiscr_ni = discrCutEfficScan->histo_ni ();
77  effVersusDiscr_dusg = discrCutEfficScan->histo_dusg();
79 
80 
81  if(mcPlots_>2){
82  effVersusDiscr_d->SetXTitle ( "Discriminant" );
83  effVersusDiscr_d->GetXaxis()->SetTitleOffset ( 0.75 );
84  effVersusDiscr_u->SetXTitle ( "Discriminant" );
85  effVersusDiscr_u->GetXaxis()->SetTitleOffset ( 0.75 );
86  effVersusDiscr_s->SetXTitle ( "Discriminant" );
87  effVersusDiscr_s->GetXaxis()->SetTitleOffset ( 0.75 );
88  effVersusDiscr_g->SetXTitle ( "Discriminant" );
89  effVersusDiscr_g->GetXaxis()->SetTitleOffset ( 0.75 );
90  effVersusDiscr_dus->SetXTitle ( "Discriminant" );
91  effVersusDiscr_dus->GetXaxis()->SetTitleOffset ( 0.75 );
92  }
93  effVersusDiscr_c->SetXTitle ( "Discriminant" );
94  effVersusDiscr_c->GetXaxis()->SetTitleOffset ( 0.75 );
95  effVersusDiscr_b->SetXTitle ( "Discriminant" );
96  effVersusDiscr_b->GetXaxis()->SetTitleOffset ( 0.75 );
97  effVersusDiscr_ni->SetXTitle ( "Discriminant" );
98  effVersusDiscr_ni->GetXaxis()->SetTitleOffset ( 0.75 );
99  effVersusDiscr_dusg->SetXTitle ( "Discriminant" );
100  effVersusDiscr_dusg->GetXaxis()->SetTitleOffset ( 0.75 );
101  effVersusDiscr_pu->SetXTitle ( "Discriminant" );
102  effVersusDiscr_pu->GetXaxis()->SetTitleOffset ( 0.75 );
103  }
104  else{
105  effVersusDiscr_d = 0;
106  effVersusDiscr_u = 0;
107  effVersusDiscr_s = 0;
108  effVersusDiscr_c = 0;
109  effVersusDiscr_b = 0;
110  effVersusDiscr_g = 0;
111  effVersusDiscr_ni = 0;
112  effVersusDiscr_dus = 0;
114  effVersusDiscr_pu = 0;
115  }
116 
117  // discr. for computation
118  vector<TH2F*> discrCfHistos = dDiscriminatorFC->getHistoVector();
119  // discr no cut
120  vector<TH2F*> discrNoCutHistos = discrNoCutEffic->getHistoVector();
121  // discr no cut
122  vector<TH2F*> discrCutHistos = discrCutEfficScan->getHistoVector();
123 
124  const int& dimHistos = discrCfHistos.size(); // they all have the same size
125 
126  // DISCR-CUT LOOP:
127  // fill the histos for eff-pur computations by scanning the discriminatorFC histogram
128 
129  // better to loop over bins -> discrCut no longer needed
130  const int& nBinsX = dDiscriminatorFC->nBinsX();
131  const int& nBinsY = dDiscriminatorFC->nBinsY();
132 
133  // loop over flavours
134  for ( int iFlav = 0; iFlav < dimHistos; iFlav++ ) {
135  if (discrCfHistos[iFlav] == 0) continue;
136  discrNoCutHistos[iFlav]->SetXTitle ( "Discriminant A" );
137  discrNoCutHistos[iFlav]->GetXaxis()->SetTitleOffset ( 0.75 );
138  discrNoCutHistos[iFlav]->SetYTitle ( "Discriminant B" );
139  discrNoCutHistos[iFlav]->GetYaxis()->SetTitleOffset ( 0.75 );
140 
141  // In Root histos, bin counting starts at 1 to nBins.
142  // bin 0 is the underflow, and nBins+1 is the overflow.
143  const double& nJetsFlav = discrCfHistos[iFlav]->GetEntries ();
144 
145  for ( int iDiscrX = nBinsX; iDiscrX > 0 ; --iDiscrX ) {
146  for ( int iDiscrY = nBinsY; iDiscrY > 0 ; --iDiscrY ) {
147  // fill all jets into NoCut histo
148  discrNoCutHistos[iFlav]->SetBinContent ( iDiscrX, iDiscrY, nJetsFlav );
149  discrNoCutHistos[iFlav]->SetBinError ( iDiscrX, iDiscrY, sqrt(nJetsFlav) );
150  const double& sum = nJetsFlav - discrCfHistos[iFlav]->Integral( 0, iDiscrX-1, 0, iDiscrY-1 );
151  discrCutHistos[iFlav]->SetBinContent ( iDiscrX, iDiscrY, sum );
152  discrCutHistos[iFlav]->SetBinError ( iDiscrX, iDiscrY, sqrt(sum) );
153  }
154  }
155  }
156 
157  // divide to get efficiency vs. discriminator cut from absolute numbers
158  discrCutEfficScan->divide ( *discrNoCutEffic ); // does: histos including discriminator cut / flat histo
159  discrCutEfficScan->setEfficiencyFlag();
160 }
161 
163 
165 {
166  if ( fromDiscriminatorDistr) {
167  discrNoCutEffic->epsPlot(name);
168  discrCutEfficScan->epsPlot(name);
169  }
170  plot(name, ".eps");
171 }
172 
174 {
175  plot(name, ".ps");
176 }
177 
179 {
180  std::string hX = "";
181  std::string Title = "";
182  if(!doCTagPlots_){
183  hX = "FlavEffVsBEff";
184  Title = "b";
185  }
186  else{
187  hX = "FlavEffVsCEff";
188  Title = "c";
189  }
190  TCanvas tc ((hX +histoExtension).c_str() ,
191  ("Flavour misidentification vs. " + Title + "-tagging efficiency " + histoExtension).c_str());
192  plot(&tc);
193  tc.Print((name + hX + histoExtension + ext).c_str());
194 }
195 
196 void EffPurFromHistos2D::plot (TPad * plotCanvas /* = 0 */) {
197 
198  setTDRStyle()->cd();
199 
200  if (plotCanvas)
201  plotCanvas->cd();
202 
203  gPad->UseCurrentStyle();
204  gPad->SetFillColor ( 0 );
205  gPad->SetLogy ( 1 );
206  gPad->SetGridx ( 1 );
207  gPad->SetGridy ( 1 );
208 }
209 
210 
212  // number of bins
213  int nBinsX_d = 0;
214  int nBinsX_u = 0;
215  int nBinsX_s = 0;
216  int nBinsX_g = 0;
217  int nBinsX_dus = 0;
218  int nBinsY_d = 0;
219  int nBinsY_u = 0;
220  int nBinsY_s = 0;
221  int nBinsY_g = 0;
222  int nBinsY_dus = 0;
223  if(mcPlots_>2){
224  nBinsX_d = effVersusDiscr_d -> GetNbinsX();
225  nBinsX_u = effVersusDiscr_u -> GetNbinsX();
226  nBinsX_s = effVersusDiscr_s -> GetNbinsX();
227  nBinsX_g = effVersusDiscr_g -> GetNbinsX();
228  nBinsX_dus = effVersusDiscr_dus -> GetNbinsX();
229  nBinsY_d = effVersusDiscr_d -> GetNbinsY();
230  nBinsY_u = effVersusDiscr_u -> GetNbinsY();
231  nBinsY_s = effVersusDiscr_s -> GetNbinsY();
232  nBinsY_g = effVersusDiscr_g -> GetNbinsY();
233  nBinsY_dus = effVersusDiscr_dus -> GetNbinsY();
234  }
235  const int& nBinsX_c = effVersusDiscr_c -> GetNbinsX();
236  const int& nBinsX_b = effVersusDiscr_b -> GetNbinsX();
237  const int& nBinsX_ni = effVersusDiscr_ni -> GetNbinsX();
238  const int& nBinsX_dusg = effVersusDiscr_dusg -> GetNbinsX();
239  const int& nBinsX_pu = effVersusDiscr_pu -> GetNbinsX();
240  const int& nBinsY_c = effVersusDiscr_c -> GetNbinsY();
241  const int& nBinsY_b = effVersusDiscr_b -> GetNbinsY();
242  const int& nBinsY_ni = effVersusDiscr_ni -> GetNbinsY();
243  const int& nBinsY_dusg = effVersusDiscr_dusg -> GetNbinsY();
244  const int& nBinsY_pu = effVersusDiscr_pu -> GetNbinsY();
245 
246  const bool& lNBinsX =
247  ( (nBinsX_d == nBinsX_u &&
248  nBinsX_d == nBinsX_s &&
249  nBinsX_d == nBinsX_c &&
250  nBinsX_d == nBinsX_b &&
251  nBinsX_d == nBinsX_g &&
252  nBinsX_d == nBinsX_ni &&
253  nBinsX_d == nBinsX_dus &&
254  nBinsX_d == nBinsX_dusg)||
255  (nBinsX_c == nBinsX_b &&
256  nBinsX_c == nBinsX_dusg &&
257  nBinsX_c == nBinsX_ni &&
258  nBinsX_c == nBinsX_pu) );
259 
260  const bool& lNBinsY =
261  ( (nBinsY_d == nBinsY_u &&
262  nBinsY_d == nBinsY_s &&
263  nBinsY_d == nBinsY_c &&
264  nBinsY_d == nBinsY_b &&
265  nBinsY_d == nBinsY_g &&
266  nBinsY_d == nBinsY_ni &&
267  nBinsY_d == nBinsY_dus &&
268  nBinsY_d == nBinsY_dusg)||
269  (nBinsY_c == nBinsY_b &&
270  nBinsY_c == nBinsY_dusg &&
271  nBinsY_c == nBinsY_ni &&
272  nBinsY_c == nBinsY_pu) );
273 
274  if ( !lNBinsX || !lNBinsY ) {
275  throw cms::Exception("Configuration")
276  << "Input histograms do not all have the same number of bins!\n";
277  }
278 
279  // start
280  float sBin_d = 0;
281  float sBin_u = 0;
282  float sBin_s = 0;
283  float sBin_g = 0;
284  float sBin_dus = 0;
285  if(mcPlots_>2){
286  sBin_d = effVersusDiscr_d -> GetBinCenter(1);
287  sBin_u = effVersusDiscr_u -> GetBinCenter(1);
288  sBin_s = effVersusDiscr_s -> GetBinCenter(1);
289  sBin_g = effVersusDiscr_g -> GetBinCenter(1);
290  sBin_dus = effVersusDiscr_dus -> GetBinCenter(1);
291  }
292  const float& sBin_c = effVersusDiscr_c -> GetBinCenter(1);
293  const float& sBin_b = effVersusDiscr_b -> GetBinCenter(1);
294  const float& sBin_ni = effVersusDiscr_ni -> GetBinCenter(1);
295  const float& sBin_dusg = effVersusDiscr_dusg -> GetBinCenter(1);
296  const float& sBin_pu = effVersusDiscr_pu -> GetBinCenter(1);
297 
298  const bool& lSBin =
299  ( (sBin_d == sBin_u &&
300  sBin_d == sBin_s &&
301  sBin_d == sBin_c &&
302  sBin_d == sBin_b &&
303  sBin_d == sBin_g &&
304  sBin_d == sBin_ni &&
305  sBin_d == sBin_dus &&
306  sBin_d == sBin_dusg)||
307  (sBin_c == sBin_b &&
308  sBin_c == sBin_dusg &&
309  sBin_c == sBin_ni &&
310  sBin_c == sBin_pu) );
311 
312  if ( !lSBin ) {
313  throw cms::Exception("Configuration")
314  << "EffPurFromHistos::check() : Input histograms do not all have the same start bin!\n";
315  }
316 
317  // end
318  float eBin_d = 0;
319  float eBin_u = 0;
320  float eBin_s = 0;
321  float eBin_g = 0;
322  float eBin_dus = 0;
323  const int& binEnd = effVersusDiscr_b -> GetBin(nBinsX_b - 1, nBinsY_b - 1);
324  if(mcPlots_>2){
325  eBin_d = effVersusDiscr_d -> GetBinCenter( binEnd );
326  eBin_u = effVersusDiscr_u -> GetBinCenter( binEnd );
327  eBin_s = effVersusDiscr_s -> GetBinCenter( binEnd );
328  eBin_g = effVersusDiscr_g -> GetBinCenter( binEnd );
329  eBin_dus = effVersusDiscr_dus -> GetBinCenter( binEnd );
330  }
331  const float& eBin_c = effVersusDiscr_c -> GetBinCenter( binEnd );
332  const float& eBin_b = effVersusDiscr_b -> GetBinCenter( binEnd );
333  const float& eBin_ni = effVersusDiscr_ni -> GetBinCenter( binEnd );
334  const float& eBin_dusg = effVersusDiscr_dusg -> GetBinCenter( binEnd );
335  const float& eBin_pu = effVersusDiscr_pu -> GetBinCenter( binEnd );
336 
337  const bool& lEBin =
338  ( (eBin_d == eBin_u &&
339  eBin_d == eBin_s &&
340  eBin_d == eBin_c &&
341  eBin_d == eBin_b &&
342  eBin_d == eBin_g &&
343  eBin_d == eBin_ni &&
344  eBin_d == eBin_dus &&
345  eBin_d == eBin_dusg)||
346  (eBin_c == eBin_b &&
347  eBin_c == eBin_dusg &&
348  eBin_c == eBin_ni &&
349  eBin_c == eBin_pu) );
350 
351  if ( !lEBin ) {
352  throw cms::Exception("Configuration")
353  << "EffPurFromHistos::check() : Input histograms do not all have the same end bin!\n";
354  }
355 }
356 
357 void EffPurFromHistos2D::compute (DQMStore::IBooker & ibook, vector<double> fixedEff)
358 {
359 
360  if (!mcPlots_ || fixedEff.size()<1) {
361  return;
362  }
363 
364  // to have shorter names ......
365  const std::string & hE = histoExtension;
366  std::string hX = "DUSG_vs_B_eff_at_fixedCeff_";
367  if(!doCTagPlots_) hX = "DUSG_vs_C_eff_at_fixedBeff_";
368 
369  // create histograms from base name and extension as given from user
370  HistoProviderDQM prov("Btag",label_,ibook);
371 
372  for (unsigned int ieff=0; ieff<fixedEff.size(); ieff++){
373  std::string fixedEfficiency = std::to_string(fixedEff[ieff]);
374  fixedEfficiency.replace(1,1,"_");
375  X_vs_Y_eff_at_fixedZeff.push_back( (prov.book1D ( hX + fixedEfficiency + hE , hX + fixedEfficiency + hE , nBinOutputX , startOutputX , endOutputX )) );
376  X_vs_Y_eff_at_fixedZeff[ieff]->setEfficiencyFlag();
377 
378  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetXTitle ( "Light mistag" );
379  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetYTitle ( "B mistag" );
380  if(!doCTagPlots_) X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetYTitle ( "C mistag" );
381  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->GetXaxis()->SetTitleOffset ( 0.75 );
382  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->GetYaxis()->SetTitleOffset ( 0.75 );
383  }
384  // loop over eff. vs. discriminator cut b-histo and look in which bin the closest entry is;
385  // use fact that eff decreases monotonously
386 
387  // any of the histos to be created can be taken here:
388  MonitorElement * EffFlavVsXEff = X_vs_Y_eff_at_fixedZeff[0];
389 
390  const int& nBinX = EffFlavVsXEff->getTH1F()->GetNbinsX();
391 
392  for ( int iBinX = 1; iBinX <= nBinX; iBinX++ ) { // loop over the bins on the x-axis of the histograms to be filled
393  const float& effBinWidthX = EffFlavVsXEff->getTH1F()->GetBinWidth ( iBinX );
394  const float& effMidX = EffFlavVsXEff->getTH1F()->GetBinCenter ( iBinX ); // middle of efficiency bin
395  const float& effLeftX = effMidX - 0.5*effBinWidthX; // left edge of bin
396  const float& effRightX = effMidX + 0.5*effBinWidthX; // right edge of bin
397 
398  vector<int> binClosest;
399  if(doCTagPlots_) binClosest = findBinClosestYValueAtFixedZ ( effVersusDiscr_dusg, effMidX , effLeftX , effRightX , effVersusDiscr_c, fixedEff );
400  else binClosest = findBinClosestYValueAtFixedZ ( effVersusDiscr_dusg, effMidX , effLeftX , effRightX , effVersusDiscr_b, fixedEff );
401 
402  for (unsigned int ieff = 0; ieff < binClosest.size(); ieff++){
403  const bool& binFound = ( binClosest[ieff] > 0 ) ;
404  //
405  if ( binFound ) {
406  // fill the histos
407  if(doCTagPlots_) {
408  X_vs_Y_eff_at_fixedZeff[ieff]->Fill(effMidX, effVersusDiscr_b->GetBinContent ( binClosest[ieff] ));
409  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetBinError(iBinX, effVersusDiscr_b->GetBinError ( binClosest[ieff] ));
410  }
411  else {
412  X_vs_Y_eff_at_fixedZeff[ieff]->Fill(effMidX, effVersusDiscr_c->GetBinContent ( binClosest[ieff] ));
413  X_vs_Y_eff_at_fixedZeff[ieff]->getTH1F()->SetBinError(iBinX, effVersusDiscr_c->GetBinError ( binClosest[ieff] ));
414  }
415  }
416  }
417  }
418 }
419 
420 
421 
422 #include <typeinfo>
double lowerBoundX() const
double lowerBoundY() const
std::string baseNameTitle() const
std::string histoExtension
virtual MonitorElement * book1D(const std::string &name, const std::string &title, const int &nchX, const double &lowX, const double &highX)
std::unique_ptr< FlavourHistograms2D< double, double > > discrCutEfficScan
T sqrt(T t)
Definition: SSEVec.h:18
std::vector< MonitorElement * > X_vs_Y_eff_at_fixedZeff
double upperBoundX() const
void epsPlot(const std::string &name)
std::vector< TH2F * > getHistoVector() const
EffPurFromHistos2D(const std::string &ext, TH2F *h_d, TH2F *h_u, TH2F *h_s, TH2F *h_c, TH2F *h_b, TH2F *h_g, TH2F *h_ni, TH2F *h_dus, TH2F *h_dusg, TH2F *h_pu, const std::string &label, const unsigned int &mc, int nBinX=100, double startOX=0.05, double endOX=1.05)
std::string baseNameDescription() const
TH1F * getTH1F(void) const
std::unique_ptr< FlavourHistograms2D< double, double > > discrNoCutEffic
void psPlot(const std::string &name)
void setTDRStyle()
Definition: plotscripts.py:87
std::vector< int > findBinClosestYValueAtFixedZ(const TH2F *, const float &yVal, const float &yLow, const float &yHigh, const TH2F *, const std::vector< double > &zVal)
Definition: Tools.cc:286
volatile std::atomic< bool > shutdown_flag false
void plot(TPad *theCanvas=0)
double upperBoundY() const
void compute(DQMStore::IBooker &ibook, std::vector< double > fixedEff)