Logo ROOT  
Reference Guide
TEfficiency.cxx
Go to the documentation of this file.
1#ifndef ROOT_TEfficiency_cxx
2#define ROOT_TEfficiency_cxx
3
4//standard header
5#include <vector>
6#include <string>
7#include <cmath>
8#include <stdlib.h>
9#include <cassert>
10
11//ROOT headers
14#include "TDirectory.h"
15#include "TF1.h"
16#include "TGraphAsymmErrors.h"
17#include "TH1.h"
18#include "TH2.h"
19#include "TH3.h"
20#include "TList.h"
21#include "TMath.h"
22#include "TROOT.h"
23#include "TStyle.h"
24#include "TVirtualPad.h"
25#include "TError.h"
28
29//custom headers
30#include "TEfficiency.h"
31
32// file with extra class for FC method
33#include "TEfficiencyHelper.h"
34
35//default values
38const Double_t kDefConfLevel = 0.682689492137; // 1 sigma
41
43
44////////////////////////////////////////////////////////////////////////////////
45/** \class TEfficiency
46 \ingroup Hist
47 \brief Class to handle efficiency histograms
48
49## I. Overview
50This class handles the calculation of efficiencies and their uncertainties. It
51provides several statistical methods for calculating frequentist and Bayesian
52confidence intervals as well as a function for combining several efficiencies.
53
54Efficiencies have a lot of applications and meanings but in principle, they can
55be described by the fraction of good/passed events k out of sample containing
56N events. One is usually interested in the dependency of the efficiency on other
57(binned) variables. The number of passed and total events is therefore stored
58internally in two histograms (TEfficiency::fTotalHistogram and TEfficiency::fPassedHistogram).
59Then the efficiency, as well as its upper and lower error, can be calculated for each bin
60individually.
61
62As the efficiency can be regarded as a parameter of a binomial distribution, the
63number of passed and total events must always be integer numbers. Therefore a
64filling with weights is not possible. However, you can assign a global weight to each
65TEfficiency object (TEfficiency::SetWeight).
66It is necessary to create one TEfficiency object
67for each weight if you investigate a process involving different weights. This
68procedure needs more effort but enables you to re-use the filled object in cases
69where you want to change one or more weights. This would not be possible if all
70events with different weights were filled in the same histogram.
71
72## II. Creating a TEfficiency object
73If you start a new analysis, it is highly recommended to use the TEfficiency class
74from the beginning. You can then use one of the constructors for fixed or
75variable bin size and your desired dimension. These constructors append the
76created TEfficiency object to the current directory. So it will be written
77automatically to a file during the next TFile::Write command.
78
79Example: create a two-dimensional TEfficiency object with
80- name = "eff"
81- title = "my efficiency"
82- axis titles: x, y and LaTeX-formatted epsilon as a label for Z axis
83- 10 bins with constant bin width (= 1) along X axis starting at 0 (lower edge
84 from the first bin) up to 10 (upper edge of last bin)
85- 20 bins with constant bin width (= 0.5) along Y axis starting at -5 (lower
86 edge from the first bin) up to 5 (upper edge of last bin)
87
88 TEfficiency* pEff = new TEfficiency("eff","my efficiency;x;y;#epsilon",10,0,10,20,-5,5);
89
90If you already have two histograms filled with the number of passed and total
91events, you will use the constructor TEfficiency(const TH1& passed,const TH1& total)
92to construct the TEfficiency object. The histograms "passed" and "total" have
93to fulfill the conditions mentioned in TEfficiency::CheckConsistency, otherwise the construction will fail.
94As the histograms already exist, the new TEfficiency is by default **not** attached
95to the current directory to avoid duplication of data. If you want to store the
96new object anyway, you can either write it directly by calling TObject::Write or attach it to a directory using TEfficiency::SetDirectory.
97This also applies to TEfficiency objects created by the copy constructor TEfficiency::TEfficiency(const TEfficiency& rEff).
98
99
100### Example 1
101
102~~~~~~~~~~~~~~~{.cpp}
103TEfficiency* pEff = 0;
104TFile* pFile = new TFile("myfile.root","recreate");
105
106//h_pass and h_total are valid and consistent histograms
107if(TEfficiency::CheckConsistency(h_pass,h_total))
108{
109 pEff = new TEfficiency(h_pass,h_total);
110 // this will write the TEfficiency object to "myfile.root"
111 // AND pEff will be attached to the current directory
112 pEff->Write();
113}
114~~~~~~~~~~~~~~~
115
116### Example 2
117
118~~~~~~~~~~~~~~~{.cpp}
119TEfficiency* pEff = 0;
120TFile* pFile = new TFile("myfile.root","recreate");
121
122//h_pass and h_total are valid and consistent histograms
123if(TEfficiency::CheckConsistency(h_pass,h_total))
124{
125 pEff = new TEfficiency(h_pass,h_total);
126 //this will attach the TEfficiency object to the current directory
127 pEff->SetDirectory(gDirectory);
128 //now all objects in gDirectory will be written to "myfile.root"
129 pFile->Write();
130}
131~~~~~~~~~~~~~~~
132
133In case you already have two filled histograms and you only want to
134plot them as a graph, you should rather use TGraphAsymmErrors::TGraphAsymmErrors(const TH1* pass,const TH1* total,Option_t* opt)
135to create a graph object.
136
137## III. Filling with events
138You can fill the TEfficiency object by calling the TEfficiency::Fill(Bool_t bPassed,Double_t x,Double_t y,Double_t z) method.
139The "bPassed" boolean flag indicates whether the current event is good
140(both histograms are filled) or not (only TEfficiency::fTotalHistogram is filled).
141The x, y and z variables determine the bin which is filled. For lower dimensions, the z- or even the y-value may be omitted.
142
143Begin_Macro(source)
144{
145 //canvas only needed for this documentation
146 TCanvas* c1 = new TCanvas("example","",600,400);
147 c1->SetFillStyle(1001);
148 c1->SetFillColor(kWhite);
149
150 //create one-dimensional TEfficiency object with fixed bin size
151 TEfficiency* pEff = new TEfficiency("eff","my efficiency;x;#epsilon",20,0,10);
152 TRandom3 rand3;
153
154 bool bPassed;
155 double x;
156 for(int i=0; i<10000; ++i)
157 {
158 //simulate events with variable under investigation
159 x = rand3.Uniform(10);
160 //check selection: bPassed = DoesEventPassSelection(x)
161 bPassed = rand3.Rndm() < TMath::Gaus(x,5,4);
162 pEff->Fill(bPassed,x);
163 }
164
165 pEff->Draw("AP");
166
167 //only for this documentation
168 return c1;
169}
170End_Macro
171
172You can also set the number of passed or total events for a bin directly by
173using the TEfficiency::SetPassedEvents or TEfficiency::SetTotalEvents method.
174
175## IV. Statistic options
176The calculation of the estimated efficiency depends on the chosen statistic
177option. Let k denotes the number of passed events and N the number of total
178events.
179
180###Frequentist methods
181The expectation value of the number of passed events is given by the true
182efficiency times the total number of events. One can estimate the efficiency
183by replacing the expected number of passed events by the observed number of
184passed events.
185
186\f[
187 k = \epsilon \times N \Rightarrow \hat{\varepsilon} = \frac{k}{N}
188\f]
189
190### Bayesian methods
191In Bayesian statistics a likelihood-function (how probable is it to get the
192observed data assuming a true efficiency) and a prior probability (what is the
193probability that a certain true efficiency is actually realised) are used to
194determine a posterior probability by using Bayes theorem. At the moment, only
195beta distributions (have 2 free parameters) are supported as prior
196probabilities.
197
198\f{eqnarray*}{
199 P(\epsilon | k ; N) &=& \frac{1}{norm} \times P(k | \epsilon ; N) \times Prior(\epsilon) \\
200 P(k | \epsilon ; N) &=& Binomial(N,k) \times \epsilon^{k} \times (1 - \epsilon)^{N - k} ...\ binomial\ distribution \\
201 Prior(\epsilon) &=& \frac{1}{B(\alpha,\beta)} \times \epsilon ^{\alpha - 1} \times (1 - \epsilon)^{\beta - 1} \equiv Beta(\epsilon; \alpha,\beta) \\
202 \Rightarrow P(\epsilon | k ; N) &=& \frac{1}{norm'} \times \epsilon^{k + \alpha - 1} \times (1 - \epsilon)^{N - k + \beta - 1} \equiv Beta(\epsilon; k + \alpha, N - k + \beta)
203\f}
204
205By default the expectation value of this posterior distribution is used as an estimator for the efficiency:
206
207\f[
208 \hat{\varepsilon} = \frac{k + \alpha}{N + \alpha + \beta}
209\f]
210
211Optionally the mode can also be used as a value for the estimated efficiency. This can be done by calling
212SetBit(kPosteriorMode) or TEfficiency::SetPosteriorMode. In this case, the estimated efficiency is:
213
214\f[
215 \hat{\varepsilon} = \frac{k + \alpha -1}{N + \alpha + \beta - 2}
216\f]
217
218In the case of a uniform prior distribution, B(x,1,1), the posterior mode is k/n, equivalent to the frequentist
219estimate (the maximum likelihood value).
220
221The statistic options also specify which confidence interval is used for calculating
222the uncertainties of the efficiency. The following properties define the error
223calculation:
224- **fConfLevel:** desired confidence level: 0 < fConfLevel < 1 (TEfficiency::GetConfidenceLevel / TEfficiency::SetConfidenceLevel)
225- **fStatisticOption** defines which method is used to calculate the boundaries of the confidence interval (TEfficiency::SetStatisticOption)
226- **fBeta_alpha, fBeta_beta:** parameters for the prior distribution which is only used in the bayesian case (TEfficiency::GetBetaAlpha / TEfficiency::GetBetaBeta / TEfficiency::SetBetaAlpha / TEfficiency::SetBetaBeta)
227- **kIsBayesian:** flag whether bayesian statistics are used or not (TEfficiency::UsesBayesianStat)
228- **kShortestInterval:** flag whether shortest interval (instead of central one) are used in case of Bayesian statistics (TEfficiency::UsesShortestInterval). Normally shortest interval should be used in combination with the mode (see TEfficiency::UsesPosteriorMode)
229- **fWeight:** global weight for this TEfficiency object which is used during combining or merging with other TEfficiency objects(TEfficiency::GetWeight / TEfficiency::SetWeight)
230
231In the following table, the implemented confidence intervals are listed
232with their corresponding statistic option. For more details on the calculation,
233please have a look at the mentioned functions.
234
235
236| name | statistic option | function | kIsBayesian | parameters |
237|------------------|------------------|---------------------|-------------|------------|
238| Clopper-Pearson | kFCP | TEfficiency::ClopperPearson |false |total events, passed events, confidence level |
239| normal approximation | kFNormal | TEfficiency::Normal | false | total events, passed events, confidence level |
240| Wilson | kFWilson | TEfficiency::Wilson | false | total events, passed events, confidence level |
241| Agresti-Coull | kFAC | TEfficiency::AgrestiCoull | false | total events, passed events. confidence level |
242| Feldman-Cousins | kFFC | TEfficiency::FeldmanCousins | false | total events, passed events, confidence level |
243| Mid-P Lancaster | kMidP | TEfficiency::MidPInterval | false | total events, passed events, confidence level |
244| Jeffrey | kBJeffrey | TEfficiency::Bayesian | true | total events, passed events, confidence level, fBeta_alpha = 0.5, fBeta_beta = 0.5 |
245| Uniform prior | kBUniform |TEfficiency::Bayesian | true |total events, passed events, confidence level, fBeta_alpha = 1, fBeta_beta = 1 |
246| custom prior | kBBayesian |TEfficiency::Bayesian | true |total events, passed events, confidence level, fBeta_alpha, fBeta_beta |
247
248The following example demonstrates the effect of different statistic options and
249confidence levels.
250
251Begin_Macro(source)
252{
253 //canvas only needed for the documentation
254 TCanvas* c1 = new TCanvas("c1","",600,400);
255 c1->Divide(2);
256 c1->SetFillStyle(1001);
257 c1->SetFillColor(kWhite);
258
259 //create one-dimensional TEfficiency object with fixed bin size
260 TEfficiency* pEff = new TEfficiency("eff","different confidence levels;x;#epsilon",20,0,10);
261 TRandom3 rand3;
262
263 bool bPassed;
264 double x;
265 for(int i=0; i<1000; ++i)
266 {
267 //simulate events with variable under investigation
268 x = rand3.Uniform(10);
269 //check selection: bPassed = DoesEventPassSelection(x)
270 bPassed = rand3.Rndm() < TMath::Gaus(x,5,4);
271 pEff->Fill(bPassed,x);
272 }
273
274 //set style attributes
275 pEff->SetFillStyle(3004);
276 pEff->SetFillColor(kRed);
277
278 //copy current TEfficiency object and set new confidence level
279 TEfficiency* pCopy = new TEfficiency(*pEff);
280 pCopy->SetConfidenceLevel(0.90);
281
282 //set style attributes
283 pCopy->SetFillStyle(3005);
284 pCopy->SetFillColor(kBlue);
285
286 c1->cd(1);
287
288 //add legend
289 TLegend* leg1 = new TLegend(0.3,0.1,0.7,0.5);
290 leg1->AddEntry(pEff,"95%","F");
291 leg1->AddEntry(pCopy,"68.3%","F");
292
293 pEff->Draw("A4");
294 pCopy->Draw("same4");
295 leg1->Draw("same");
296
297 //use same confidence level but different statistic methods
298 TEfficiency* pEff2 = new TEfficiency(*pEff);
299 TEfficiency* pCopy2 = new TEfficiency(*pEff);
300
301 pEff2->SetStatisticOption(TEfficiency::kFNormal);
302 pCopy2->SetStatisticOption(TEfficiency::kFAC);
303
304 pEff2->SetTitle("different statistic options;x;#epsilon");
305
306 //set style attributes
307 pCopy2->SetFillStyle(3005);
308 pCopy2->SetFillColor(kBlue);
309
310 c1->cd(2);
311
312 //add legend
313 TLegend* leg2 = new TLegend(0.3,0.1,0.7,0.5);
314 leg2->AddEntry(pEff2,"kFNormal","F");
315 leg2->AddEntry(pCopy2,"kFAC","F");
316
317 pEff2->Draw("a4");
318 pCopy2->Draw("same4");
319 leg2->Draw("same");
320
321 //only for this documentation
322 c1->cd(0);
323 return c1;
324}
325End_Macro
326
327The prior probability of the efficiency in Bayesian statistics can be given
328in terms of a beta distribution. The beta distribution has two positive shape
329parameters. The resulting priors for different combinations of these shape
330parameters are shown in the plot below.
331
332Begin_Macro(source)
333{
334 //canvas only needed for the documentation
335 TCanvas* c1 = new TCanvas("c1","",600,400);
336 c1->SetFillStyle(1001);
337 c1->SetFillColor(kWhite);
338
339 //create different beta distributions
340 TF1* f1 = new TF1("f1","TMath::BetaDist(x,1,1)",0,1);
341 f1->SetLineColor(kBlue);
342 TF1* f2 = new TF1("f2","TMath::BetaDist(x,0.5,0.5)",0,1);
343 f2->SetLineColor(kRed);
344 TF1* f3 = new TF1("f3","TMath::BetaDist(x,1,5)",0,1);
345 f3->SetLineColor(kGreen+3);
346 f3->SetTitle("Beta distributions as priors;#epsilon;P(#epsilon)");
347 TF1* f4 = new TF1("f4","TMath::BetaDist(x,4,3)",0,1);
348 f4->SetLineColor(kViolet);
349
350 //add legend
351 TLegend* leg = new TLegend(0.25,0.5,0.85,0.89);
352 leg->SetFillColor(kWhite);
353 leg->SetFillStyle(1001);
354 leg->AddEntry(f1,"a=1, b=1","L");
355 leg->AddEntry(f2,"a=0.5, b=0.5","L");
356 leg->AddEntry(f3,"a=1, b=5","L");
357 leg->AddEntry(f4,"a=4, b=3","L");
358
359 f3->Draw();
360 f1->Draw("same");
361 f2->Draw("Same");
362 f4->Draw("same");
363 leg->Draw("same");
364
365 //only for this documentation
366 return c1;
367}
368End_Macro
369
370
371## IV.1 Coverage probabilities for different methods
372The following pictures illustrate the actual coverage probability for the
373different values of the true efficiency and the total number of events when a
374confidence level of 95% is desired.
375
376\image html normal95.gif "Normal Approximation"
377
378
379\image html wilson95.gif "Wilson"
380
381
382\image html ac95.gif "Agresti Coull"
383
384
385\image html cp95.gif "Clopper Pearson"
386
387
388\image html uni95.gif "Bayesian with Uniform Prior"
389
390
391\image html jeffrey95.gif "Bayesian with Jeffrey Prior"
392
393The average (over all possible true efficiencies) coverage probability for
394different number of total events is shown in the next picture.
395\image html av_cov.png "Average Coverage"
396
397## V. Merging and combining TEfficiency objects
398In many applications, the efficiency should be calculated for an inhomogeneous
399sample in the sense that it contains events with different weights. In order
400to be able to determine the correct overall efficiency, it is necessary to
401use for each subsample (= all events with the same weight) a different
402TEfficiency object. After finishing your analysis you can then construct the
403overall efficiency with its uncertainty.
404
405This procedure has the advantage that you can change the weight of one
406subsample easily without rerunning the whole analysis. On the other hand, more
407effort is needed to handle several TEfficiency objects instead of one
408histogram. In the case of many different or even continuously distributed
409weights, this approach becomes cumbersome. One possibility to overcome this
410problem is the usage of binned weights.
411
412### Example
413In particle physics weights arises from the fact that you want to
414normalise your results to a certain reference value. A very common formula for
415calculating weights is
416
417\f{eqnarray*}{
418 w &=& \frac{\sigma L}{N_{gen} \epsilon_{trig}} \\
419 &-& \sigma ...\ cross\ section \\
420 &-& L ...\ luminosity \\
421 &-& N_{gen}\ ... number\ of\ generated\ events \\
422 &-& \epsilon_{trig}\ ...\ (known)\ trigger\ efficiency \\
423\f}
424
425The reason for different weights can therefore be:
426- different processes
427- other integrated luminosity
428- varying trigger efficiency
429- different sample sizes
430- ...
431- or even combination of them
432
433Depending on the actual meaning of different weights in your case, you
434should either merge or combine them to get the overall efficiency.
435
436### V.1 When should I use merging?
437If the weights are artificial and do not represent real alternative hypotheses,
438you should merge the different TEfficiency objects. That means especially for
439the Bayesian case that the prior probability should be the same for all merged
440TEfficiency objects. The merging can be done by invoking one of the following
441operations:
442- eff1.Add(eff2)
443- eff1 += eff2
444- eff1 = eff1 + eff2
445
446The result of the merging is stored in the TEfficiency object which is marked
447bold above. The contents of the internal histograms of both TEfficiency
448objects are added and a new weight is assigned. The statistic options are not
449changed.
450
451\f[
452 \frac{1}{w_{new}} = \frac{1}{w_{1}} + \frac{1}{w_{2}}
453\f]
454
455### Example:
456If you use two samples with different numbers of generated events for the same
457process and you want to normalise both to the same integrated luminosity and
458trigger efficiency, the different weights then arise just from the fact that
459you have different numbers of events. The TEfficiency objects should be merged
460because the samples do not represent true alternatives. You expect the same
461result as if you would have a big sample with all events in it.
462
463\f[
464 w_{1} = \frac{\sigma L}{\epsilon N_{1}}, w_{2} = \frac{\sigma L}{\epsilon N_{2}} \Rightarrow w_{new} = \frac{\sigma L}{\epsilon (N_{1} + N_{2})} = \frac{1}{\frac{1}{w_{1}} + \frac{1}{w_{2}}}
465\f]
466
467### V.2 When should I use combining?
468You should combine TEfficiency objects whenever the weights represent
469alternatives processes for the efficiency. As the combination of two TEfficiency
470objects is not always consistent with the representation by two internal
471histograms, the result is not stored in a TEfficiency object but a TGraphAsymmErrors
472is returned which shows the estimated combined efficiency and its uncertainty
473for each bin.
474At the moment the combination method TEfficiency::Combine only supports a combination of 1-dimensional
475efficiencies in a Bayesian approach.
476
477
478For calculating the combined efficiency and its uncertainty for each bin only Bayesian statistics
479is used. No frequentists methods are presently supported for computing the combined efficiency and
480its confidence interval.
481In the case of the Bayesian statistics, a combined posterior is constructed taking into account the
482weight of each TEfficiency object. The same prior is used for all the TEfficiency objects.
483
484\f{eqnarray*}{
485 P_{comb}(\epsilon | {w_{i}}, {k_{i}} , {N_{i}}) = \frac{1}{norm} \prod_{i}{L(k_{i} | N_{i}, \epsilon)}^{w_{i}} \Pi( \epsilon )\\
486L(k_{i} | N_{i}, \epsilon)\ is\ the\ likelihood\ function\ for\ the\ sample\ i\ (a\ Binomial\ distribution)\\
487\Pi( \epsilon)\ is\ the\ prior,\ a\ beta\ distribution\ B(\epsilon, \alpha, \beta).\\
488The\ resulting\ combined\ posterior\ is \\
489P_{comb}(\epsilon |{w_{i}}; {k_{i}}; {N_{i}}) = B(\epsilon, \sum_{i}{ w_{i} k_{i}} + \alpha, \sum_{i}{ w_{i}(n_{i}-k_{i})}+\beta) \\
490\hat{\varepsilon} = \int_{0}^{1} \epsilon \times P_{comb}(\epsilon | {k_{i}} , {N_{i}}) d\epsilon \\
491confidence\ level = 1 - \alpha \\
492\frac{\alpha}{2} = \int_{0}^{\epsilon_{low}} P_{comb}(\epsilon | {k_{i}} , {N_{i}}) d\epsilon ...\ defines\ lower\ boundary \\
4931- \frac{\alpha}{2} = \int_{0}^{\epsilon_{up}} P_{comb}(\epsilon | {k_{i}} , {N_{i}}) d\epsilon ...\ defines\ upper\ boundary
494\f}
495
496
497###Example:
498If you use cuts to select electrons which can originate from two different
499processes, you can determine the selection efficiency for each process. The
500overall selection efficiency is then the combined efficiency. The weights to be used in the
501combination should be the probability that an
502electron comes from the corresponding process.
503
504\f[
505p_{1} = \frac{\sigma_{1}}{\sigma_{1} + \sigma_{2}} = \frac{N_{1}w_{1}}{N_{1}w_{1} + N_{2}w_{2}}\\
506p_{2} = \frac{\sigma_{2}}{\sigma_{1} + \sigma_{2}} = \frac{N_{2}w_{2}}{N_{1}w_{1} + N_{2}w_{2}}
507\f]
508
509## VI. Further operations
510
511### VI.Information about the internal histograms
512The methods TEfficiency::GetPassedHistogram and TEfficiency::GetTotalHistogram
513return a constant pointer to the internal histograms. They can be used to
514obtain information about the internal histograms (e.g., the binning, number of passed / total events in a bin, mean values...).
515One can obtain a clone of the internal histograms by calling TEfficiency::GetCopyPassedHisto or TEfficiency::GetCopyTotalHisto.
516The returned histograms are completely independent from the current
517TEfficiency object. By default, they are not attached to a directory to
518avoid the duplication of data and the user is responsible for deleting them.
519
520
521~~~~~~~~~~~~~~~{.cpp}
522//open a root file which contains a TEfficiency object
523TFile* pFile = new TFile("myfile.root","update");
524
525//get TEfficiency object with name "my_eff"
526TEfficiency* pEff = (TEfficiency*)pFile->Get("my_eff");
527
528//get clone of total histogram
529TH1* clone = pEff->GetCopyTotalHisto();
530
531//change clone...
532//save changes of clone directly
533clone->Write();
534//or append it to the current directory and write the file
535//clone->SetDirectory(gDirectory);
536//pFile->Write();
537
538//delete histogram object
539delete clone;
540clone = 0;
541~~~~~~~~~~~~~~~
542
543It is also possible to set the internal total or passed histogram by using the
544methods TEfficiency::SetPassedHistogram or TEfficiency::SetTotalHistogram.
545
546In order to ensure the validity of the TEfficiency object, the consistency of the
547new histogram and the stored histogram is checked. It might be
548impossible sometimes to change the histograms in a consistent way. Therefore one can force
549the replacement by passing the "f" option. Then the user has to ensure that the
550other internal histogram is replaced as well and that the TEfficiency object is
551in a valid state.
552
553### VI.2 Fitting
554The efficiency can be fitted using the TEfficiency::Fit function which internally uses
555the TBinomialEfficiencyFitter::Fit method.
556As this method is using a maximum-likelihood-fit, it is necessary to initialise
557the given fit function with reasonable start values.
558The resulting fit function is attached to the list of associated functions and
559will be drawn automatically during the next TEfficiency::Draw command.
560The list of associated function can be modified by using the pointer returned
561by TEfficiency::GetListOfFunctions.
562
563Begin_Macro(source)
564{
565 //canvas only needed for this documentation
566 TCanvas* c1 = new TCanvas("example","",600,400);
567 c1->SetFillStyle(1001);
568 c1->SetFillColor(kWhite);
569
570 //create one-dimensional TEfficiency object with fixed bin size
571 TEfficiency* pEff = new TEfficiency("eff","my efficiency;x;#epsilon",20,0,10);
572 TRandom3 rand3;
573
574 bool bPassed;
575 double x;
576 for(int i=0; i<10000; ++i)
577 {
578 //simulate events with variable under investigation
579 x = rand3.Uniform(10);
580 //check selection: bPassed = DoesEventPassSelection(x)
581 bPassed = rand3.Rndm() < TMath::Gaus(x,5,4);
582 pEff->Fill(bPassed,x);
583 }
584
585 //create a function for fitting and do the fit
586 TF1* f1 = new TF1("f1","gaus",0,10);
587 f1->SetParameters(1,5,2);
588 pEff->Fit(f1);
589
590 //create a threshold function
591 TF1* f2 = new TF1("thres","0.8",0,10);
592 f2->SetLineColor(kRed);
593 //add it to the list of functions
594 //use add first because the parameters of the last function will be displayed
595 pEff->GetListOfFunctions()->AddFirst(f2);
596
597 pEff->Draw("AP");
598
599 //only for this documentation
600 return c1;
601}
602End_Macro
603
604### VI.3 Draw a TEfficiency object
605A TEfficiency object can be drawn by calling the usual TEfficiency::Draw method.
606At the moment drawing is only supported for 1- and 2-dimensional TEfficiency objects.
607In the 1-dimensional case, you can use the same options as for the TGraphAsymmErrors::Draw
608method. For 2-dimensional TEfficiency objects, you can pass the same options as
609for a TH2::Draw object.
610
611*/
612
613////////////////////////////////////////////////////////////////////////////////
614/// Default constructor
615///
616/// Should not be used explicitly
617
619fBeta_alpha(kDefBetaAlpha),
620fBeta_beta(kDefBetaBeta),
621fBoundary(0),
622fConfLevel(kDefConfLevel),
623fDirectory(0),
624fFunctions(0),
625fPaintGraph(0),
626fPaintHisto(0),
627fPassedHistogram(0),
628fTotalHistogram(0),
629fWeight(kDefWeight)
630{
632
633 // create 2 dummy histograms
634 fPassedHistogram = new TH1F("h_passed","passed",10,0,10);
635 fTotalHistogram = new TH1F("h_total","total",10,0,10);
636}
637
638////////////////////////////////////////////////////////////////////////////////
639/// Constructor using two existing histograms as input
640///
641///Input: passed - contains the events fulfilling some criteria
642/// total - contains all investigated events
643///
644///Notes: - both histograms have to fulfill the conditions of CheckConsistency
645/// - dimension of the resulting efficiency object depends
646/// on the dimension of the given histograms
647/// - Clones of both histograms are stored internally
648/// - The function SetName(total.GetName() + "_clone") is called to set
649/// the names of the new object and the internal histograms..
650/// - The created TEfficiency object is NOT appended to a directory. It
651/// will not be written to disk during the next TFile::Write() command
652/// in order to prevent duplication of data. If you want to save this
653/// TEfficiency object anyway, you can either append it to a
654/// directory by calling SetDirectory(TDirectory*) or write it
655/// explicitly to disk by calling Write().
656
657TEfficiency::TEfficiency(const TH1& passed,const TH1& total):
658fBeta_alpha(kDefBetaAlpha),
659fBeta_beta(kDefBetaBeta),
660fConfLevel(kDefConfLevel),
661fDirectory(0),
662fFunctions(0),
663fPaintGraph(0),
664fPaintHisto(0),
665fWeight(kDefWeight)
666{
667 //check consistency of histograms
668 if(CheckConsistency(passed,total)) {
669 // do not add cloned histograms to gDirectory
670 {
671 TDirectory::TContext ctx(nullptr);
672 fTotalHistogram = (TH1*)total.Clone();
673 fPassedHistogram = (TH1*)passed.Clone();
674 }
675
676 TString newName = total.GetName();
677 newName += TString("_clone");
678 SetName(newName);
679
680 // are the histograms filled with weights?
681 if(CheckWeights(passed,total))
682 {
683 Info("TEfficiency","given histograms are filled with weights");
685 }
686 }
687 else {
688 Error("TEfficiency(const TH1&,const TH1&)","histograms are not consistent -> results are useless");
689 Warning("TEfficiency(const TH1&,const TH1&)","using two empty TH1D('h1','h1',10,0,10)");
690
691 // do not add new created histograms to gDirectory
692 TDirectory::TContext ctx(nullptr);
693 fTotalHistogram = new TH1D("h1_total","h1 (total)",10,0,10);
694 fPassedHistogram = new TH1D("h1_passed","h1 (passed)",10,0,10);
695 }
696
697 SetBit(kPosteriorMode,false);
699
701 SetDirectory(0);
702}
703
704////////////////////////////////////////////////////////////////////////////////
705/// Create 1-dimensional TEfficiency object with variable bin size.
706///
707/// Constructor creates two new and empty histograms with a given binning
708///
709/// Input:
710///
711/// - `name`: the common part of the name for both histograms (no blanks)
712/// fTotalHistogram has name: name + "_total"
713/// fPassedHistogram has name: name + "_passed"
714/// - `title`: the common part of the title for both histogram
715/// fTotalHistogram has title: title + " (total)"
716/// fPassedHistogram has title: title + " (passed)"
717/// It is possible to label the axis by passing a title with
718/// the following format: "title;xlabel;ylabel".
719/// - `nbins`: number of bins on the x-axis
720/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
721/// xbins[nbinsx] ... lower edge for overflow bin
722
723TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbins,
724 const Double_t* xbins):
725fBeta_alpha(kDefBetaAlpha),
726fBeta_beta(kDefBetaBeta),
727fConfLevel(kDefConfLevel),
728fDirectory(0),
729fFunctions(0),
730fPaintGraph(0),
731fPaintHisto(0),
732fWeight(kDefWeight)
733{
734 // do not add new created histograms to gDirectory
735 {
736 // use separate scope for TContext
737 TDirectory::TContext ctx(nullptr);
738 fTotalHistogram = new TH1D("total","total",nbins,xbins);
739 fPassedHistogram = new TH1D("passed","passed",nbins,xbins);
740 }
741
742 Build(name,title);
743}
744
745////////////////////////////////////////////////////////////////////////////////
746/// Create 1-dimensional TEfficiency object with fixed bins size.
747///
748/// Constructor creates two new and empty histograms with a fixed binning.
749///
750/// Input:
751///
752/// - `name`: the common part of the name for both histograms(no blanks)
753/// fTotalHistogram has name: name + "_total"
754/// fPassedHistogram has name: name + "_passed"
755/// - `title`: the common part of the title for both histogram
756/// fTotalHistogram has title: title + " (total)"
757/// fPassedHistogram has title: title + " (passed)"
758/// It is possible to label the axis by passing a title with
759/// the following format: "title;xlabel;ylabel".
760/// - `nbinsx`: number of bins on the x-axis
761/// - `xlow`: lower edge of first bin
762/// - `xup`: upper edge of last bin
763
764TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
765 Double_t xlow,Double_t xup):
766fBeta_alpha(kDefBetaAlpha),
767fBeta_beta(kDefBetaBeta),
768fConfLevel(kDefConfLevel),
769fDirectory(0),
770fFunctions(0),
771fPaintGraph(0),
772fPaintHisto(0),
773fWeight(kDefWeight)
774{
775 // do not add new created histograms to gDirectory
776 {
777 TDirectory::TContext ctx(nullptr);
778 fTotalHistogram = new TH1D("total","total",nbinsx,xlow,xup);
779 fPassedHistogram = new TH1D("passed","passed",nbinsx,xlow,xup);
780 }
781 Build(name,title);
782}
783
784////////////////////////////////////////////////////////////////////////////////
785/// Create 2-dimensional TEfficiency object with fixed bin size.
786///
787/// Constructor creates two new and empty histograms with a fixed binning.
788///
789/// Input:
790///
791/// - `name`: the common part of the name for both histograms(no blanks)
792/// fTotalHistogram has name: name + "_total"
793/// fPassedHistogram has name: name + "_passed"
794/// - `title`: the common part of the title for both histogram
795/// fTotalHistogram has title: title + " (total)"
796/// fPassedHistogram has title: title + " (passed)"
797/// It is possible to label the axis by passing a title with
798/// the following format: "title;xlabel;ylabel;zlabel".
799/// - `nbinsx`: number of bins on the x-axis
800/// - `xlow`: lower edge of first x-bin
801/// - `xup`: upper edge of last x-bin
802/// - `nbinsy`: number of bins on the y-axis
803/// - `ylow`: lower edge of first y-bin
804/// - `yup`: upper edge of last y-bin
805
806TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
807 Double_t xlow,Double_t xup,Int_t nbinsy,
808 Double_t ylow,Double_t yup):
809fBeta_alpha(kDefBetaAlpha),
810fBeta_beta(kDefBetaBeta),
811fConfLevel(kDefConfLevel),
812fDirectory(0),
813fFunctions(0),
814fPaintGraph(0),
815fPaintHisto(0),
816fWeight(kDefWeight)
817{
818 // do not add new created histograms to gDirectory
819 {
820 TDirectory::TContext ctx(nullptr);
821 fTotalHistogram = new TH2D("total","total",nbinsx,xlow,xup,nbinsy,ylow,yup);
822 fPassedHistogram = new TH2D("passed","passed",nbinsx,xlow,xup,nbinsy,ylow,yup);
823 }
824 Build(name,title);
825}
826
827////////////////////////////////////////////////////////////////////////////////
828/// Create 2-dimensional TEfficiency object with variable bin size.
829///
830/// Constructor creates two new and empty histograms with a given binning.
831///
832/// Input:
833///
834/// - `name`: the common part of the name for both histograms(no blanks)
835/// fTotalHistogram has name: name + "_total"
836/// fPassedHistogram has name: name + "_passed"
837/// - `title`: the common part of the title for both histogram
838/// fTotalHistogram has title: title + " (total)"
839/// fPassedHistogram has title: title + " (passed)"
840/// It is possible to label the axis by passing a title with
841/// the following format: "title;xlabel;ylabel;zlabel".
842/// - `nbinsx`: number of bins on the x-axis
843/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
844/// xbins[nbinsx] ... lower edge for overflow x-bin
845/// - `nbinsy`: number of bins on the y-axis
846/// - `ybins`: array of length (nbins + 1) with low-edges for each bin
847/// ybins[nbinsy] ... lower edge for overflow y-bin
848
849TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
850 const Double_t* xbins,Int_t nbinsy,
851 const Double_t* ybins):
852fBeta_alpha(kDefBetaAlpha),
853fBeta_beta(kDefBetaBeta),
854fConfLevel(kDefConfLevel),
855fDirectory(0),
856fFunctions(0),
857fPaintGraph(0),
858fPaintHisto(0),
859fWeight(kDefWeight)
860{
861 // do not add new created histograms to gDirectory
862 {
863 TDirectory::TContext ctx(nullptr);
864 fTotalHistogram = new TH2D("total","total",nbinsx,xbins,nbinsy,ybins);
865 fPassedHistogram = new TH2D("passed","passed",nbinsx,xbins,nbinsy,ybins);
866 }
867 Build(name,title);
868}
869
870////////////////////////////////////////////////////////////////////////////////
871/// Create 3-dimensional TEfficiency object with fixed bin size.
872///
873/// Constructor creates two new and empty histograms with a fixed binning.
874///
875/// Input:
876///
877/// - `name`: the common part of the name for both histograms(no blanks)
878/// fTotalHistogram has name: name + "_total"
879/// fPassedHistogram has name: name + "_passed"
880/// - `title`: the common part of the title for both histogram
881/// fTotalHistogram has title: title + " (total)"
882/// fPassedHistogram has title: title + " (passed)"
883/// It is possible to label the axis by passing a title with
884/// the following format: "title;xlabel;ylabel;zlabel".
885/// - `nbinsx`: number of bins on the x-axis
886/// - `xlow`: lower edge of first x-bin
887/// - `xup`: upper edge of last x-bin
888/// - `nbinsy`: number of bins on the y-axis
889/// - `ylow`: lower edge of first y-bin
890/// - `yup`: upper edge of last y-bin
891/// - `nbinsz`: number of bins on the z-axis
892/// - `zlow`: lower edge of first z-bin
893/// - `zup`: upper edge of last z-bin
894
895TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
896 Double_t xlow,Double_t xup,Int_t nbinsy,
897 Double_t ylow,Double_t yup,Int_t nbinsz,
898 Double_t zlow,Double_t zup):
899fBeta_alpha(kDefBetaAlpha),
900fBeta_beta(kDefBetaBeta),
901fConfLevel(kDefConfLevel),
902fDirectory(0),
903fFunctions(0),
904fPaintGraph(0),
905fPaintHisto(0),
906fWeight(kDefWeight)
907{
908 // do not add new created histograms to gDirectory
909 {
910 TDirectory::TContext ctx(nullptr);
911 fTotalHistogram = new TH3D("total","total",nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup);
912 fPassedHistogram = new TH3D("passed","passed",nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup);
913 }
914 Build(name,title);
915}
916
917////////////////////////////////////////////////////////////////////////////////
918/// Create 3-dimensional TEfficiency object with variable bin size.
919///
920/// Constructor creates two new and empty histograms with a given binning.
921///
922/// Input:
923///
924/// - `name`: the common part of the name for both histograms(no blanks)
925/// fTotalHistogram has name: name + "_total"
926/// fPassedHistogram has name: name + "_passed"
927/// - `title`: the common part of the title for both histogram
928/// fTotalHistogram has title: title + " (total)"
929/// fPassedHistogram has title: title + " (passed)"
930/// It is possible to label the axis by passing a title with
931/// the following format: "title;xlabel;ylabel;zlabel".
932/// - `nbinsx`: number of bins on the x-axis
933/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
934/// xbins[nbinsx] ... lower edge for overflow x-bin
935/// - `nbinsy`: number of bins on the y-axis
936/// - `ybins`: array of length (nbins + 1) with low-edges for each bin
937/// xbins[nbinsx] ... lower edge for overflow y-bin
938/// - `nbinsz`: number of bins on the z-axis
939/// - `zbins`: array of length (nbins + 1) with low-edges for each bin
940/// xbins[nbinsx] ... lower edge for overflow z-bin
941
942TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
943 const Double_t* xbins,Int_t nbinsy,
944 const Double_t* ybins,Int_t nbinsz,
945 const Double_t* zbins):
946fBeta_alpha(kDefBetaAlpha),
947fBeta_beta(kDefBetaBeta),
948fConfLevel(kDefConfLevel),
949fDirectory(0),
950fFunctions(0),
951fPaintGraph(0),
952fPaintHisto(0),
953fWeight(kDefWeight)
954{
955 // do not add new created histograms to gDirectory
956 {
957 TDirectory::TContext ctx(nullptr);
958 fTotalHistogram = new TH3D("total","total",nbinsx,xbins,nbinsy,ybins,nbinsz,zbins);
959 fPassedHistogram = new TH3D("passed","passed",nbinsx,xbins,nbinsy,ybins,nbinsz,zbins);
960 }
961 Build(name,title);
962}
963
964////////////////////////////////////////////////////////////////////////////////
965/// Copy constructor.
966///
967///The list of associated objects (e.g. fitted functions) is not copied.
968///
969///Note:
970///
971/// - SetName(rEff.GetName() + "_copy") is called to set the names of the
972/// object and the histograms.
973/// - The titles are set by calling SetTitle("[copy] " + rEff.GetTitle()).
974/// - The copied TEfficiency object is NOT appended to a directory. It
975/// will not be written to disk during the next TFile::Write() command
976/// in order to prevent duplication of data. If you want to save this
977/// TEfficiency object anyway, you can either append it to a directory
978/// by calling SetDirectory(TDirectory*) or write it explicitly to disk
979/// by calling Write().
980
982TNamed(),
983TAttLine(),
984TAttFill(),
985TAttMarker(),
986fBeta_alpha(rEff.fBeta_alpha),
987fBeta_beta(rEff.fBeta_beta),
988fBeta_bin_params(rEff.fBeta_bin_params),
989fConfLevel(rEff.fConfLevel),
990fDirectory(0),
991fFunctions(0),
992fPaintGraph(0),
993fPaintHisto(0),
994fWeight(rEff.fWeight)
995{
996 // copy TObject bits
997 ((TObject&)rEff).Copy(*this);
998
999 // do not add cloned histograms to gDirectory
1000 {
1001 TDirectory::TContext ctx(nullptr);
1002 fTotalHistogram = (TH1*)((rEff.fTotalHistogram)->Clone());
1003 fPassedHistogram = (TH1*)((rEff.fPassedHistogram)->Clone());
1004 }
1005
1006 TString name = rEff.GetName();
1007 name += "_copy";
1008 SetName(name);
1009 TString title = "[copy] ";
1010 title += rEff.GetTitle();
1011 SetTitle(title);
1012
1014
1015 SetDirectory(0);
1016
1017 //copy style
1018 rEff.TAttLine::Copy(*this);
1019 rEff.TAttFill::Copy(*this);
1020 rEff.TAttMarker::Copy(*this);
1021}
1022
1023////////////////////////////////////////////////////////////////////////////////
1024///default destructor
1025
1027{
1028 //delete all function in fFunctions
1029 // use same logic as in TH1 destructor
1030 // (see TH1::~TH1 code in TH1.cxx)
1031 if(fFunctions) {
1033 TObject* obj = 0;
1034 while ((obj = fFunctions->First())) {
1035 while(fFunctions->Remove(obj)) { }
1036 if (!obj->TestBit(kNotDeleted)) {
1037 break;
1038 }
1039 delete obj;
1040 obj = 0;
1041 }
1042 delete fFunctions;
1043 fFunctions = 0;
1044 }
1045
1046 if(fDirectory)
1047 fDirectory->Remove(this);
1048
1049 delete fTotalHistogram;
1050 delete fPassedHistogram;
1051 delete fPaintGraph;
1052 delete fPaintHisto;
1053}
1054
1055////////////////////////////////////////////////////////////////////////////////
1056/**
1057 Calculates the boundaries for the frequentist Agresti-Coull interval
1058
1059 \param total number of total events
1060 \param passed 0 <= number of passed events <= total
1061 \param level confidence level
1062 \param bUpper true - upper boundary is returned
1063 false - lower boundary is returned
1064
1065
1066 \f{eqnarray*}{
1067 \alpha &=& 1 - \frac{level}{2} \\
1068 \kappa &=& \Phi^{-1}(1 - \alpha,1)\ ... normal\ quantile\ function\\
1069 mode &=& \frac{passed + \frac{\kappa^{2}}{2}}{total + \kappa^{2}}\\
1070 \Delta &=& \kappa * \sqrt{\frac{mode * (1 - mode)}{total + \kappa^{2}}}\\
1071 return &=& max(0,mode - \Delta)\ or\ min(1,mode + \Delta)
1072 \f}
1073
1074*/
1075
1077{
1078 Double_t alpha = (1.0 - level)/2;
1079 Double_t kappa = ROOT::Math::normal_quantile(1 - alpha,1);
1080
1081 Double_t mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa);
1082 Double_t delta = kappa * std::sqrt(mode * (1 - mode) / (total + kappa * kappa));
1083
1084 if(bUpper)
1085 return ((mode + delta) > 1) ? 1.0 : (mode + delta);
1086 else
1087 return ((mode - delta) < 0) ? 0.0 : (mode - delta);
1088}
1089
1090////////////////////////////////////////////////////////////////////////////////
1091/// Calculates the boundaries for the frequentist Feldman-Cousins interval
1092///
1093/// \param total number of total events
1094/// \param passed 0 <= number of passed events <= total
1095/// \param level confidence level
1096/// \param bUpper: true - upper boundary is returned
1097/// false - lower boundary is returned
1098
1100{
1101 Double_t lower = 0;
1102 Double_t upper = 1;
1103 if (!FeldmanCousinsInterval(total,passed,level, lower, upper)) {
1104 ::Error("FeldmanCousins","Error running FC method - return 0 or 1");
1105 }
1106 return (bUpper) ? upper : lower;
1107}
1108
1109////////////////////////////////////////////////////////////////////////////////
1110/// Calculates the interval boundaries using the frequentist methods of Feldman-Cousins
1111///
1112/// \param[in] total number of total events
1113/// \param[in] passed 0 <= number of passed events <= total
1114/// \param[in] level confidence level
1115/// \param[out] lower lower boundary returned on exit
1116/// \param[out] upper lower boundary returned on exit
1117/// \return a flag with the status of the calculation
1118///
1119/// Calculation:
1120///
1121/// The Feldman-Cousins is a frequentist method where the interval is estimated using a Neyman construction where the ordering
1122/// is based on the likelihood ratio:
1123/// \f[
1124/// LR = \frac{Binomial(k | N, \epsilon)}{Binomial(k | N, \hat{\epsilon} ) }
1125/// \f]
1126/// See G. J. Feldman and R. D. Cousins, Phys. Rev. D57 (1998) 3873
1127/// and R. D. Cousins, K. E. Hymes, J. Tucker, Nuclear Instruments and Methods in Physics Research A 612 (2010) 388
1128///
1129/// Implemented using classes developed by Jordan Tucker and Luca Lista
1130/// See File hist/hist/src/TEfficiencyHelper.h
1131
1133{
1135 double alpha = 1.-level;
1136 fc.Init(alpha);
1137 fc.Calculate(passed, total);
1138 lower = fc.Lower();
1139 upper = fc.Upper();
1140 return true;
1141}
1142
1143////////////////////////////////////////////////////////////////////////////////
1144/// Calculates the boundaries using the mid-P binomial
1145/// interval (Lancaster method) from B. Cousing and J. Tucker.
1146/// See http://arxiv.org/abs/0905.3831 for a description and references for the method
1147///
1148/// Modify equal_tailed to get the kind of interval you want.
1149/// Can also be converted to interval on ratio of poisson means X/Y by the substitutions
1150/// ~~~ {.cpp}
1151/// X = passed
1152/// total = X + Y
1153/// lower_poisson = lower/(1 - lower)
1154/// upper_poisson = upper/(1 - upper)
1155/// ~~~
1156
1158{
1159 const double alpha = 1. - level;
1160 const bool equal_tailed = true; // change if you don;t want equal tailed interval
1161 const double alpha_min = equal_tailed ? alpha/2 : alpha;
1162 const double tol = 1e-9; // tolerance
1163 double pmin = 0;
1164 double pmax = 0;
1165 double p = 0;
1166
1167 pmin = 0; pmax = 1;
1168
1169
1170 // treat special case for 0<passed<1
1171 // do a linear interpolation of the upper limit values
1172 if ( passed > 0 && passed < 1) {
1173 double p0 = MidPInterval(total,0.0,level,bUpper);
1174 double p1 = MidPInterval(total,1.0,level,bUpper);
1175 p = (p1 - p0) * passed + p0;
1176 return p;
1177 }
1178
1179 while (std::abs(pmax - pmin) > tol) {
1180 p = (pmin + pmax)/2;
1181 //double v = 0.5 * ROOT::Math::binomial_pdf(int(passed), p, int(total));
1182 // make it work for non integer using the binomial - beta relationship
1183 double v = 0.5 * ROOT::Math::beta_pdf(p, passed+1., total-passed+1)/(total+1);
1184 //if (passed > 0) v += ROOT::Math::binomial_cdf(int(passed - 1), p, int(total));
1185 // compute the binomial cdf at passed -1
1186 if ( (passed-1) >= 0) v += ROOT::Math::beta_cdf_c(p, passed, total-passed+1);
1187
1188 double vmin = (bUpper) ? alpha_min : 1.- alpha_min;
1189 if (v > vmin)
1190 pmin = p;
1191 else
1192 pmax = p;
1193 }
1194
1195 return p;
1196}
1197
1198
1199////////////////////////////////////////////////////////////////////////////////
1200/**
1201Calculates the boundaries for a Bayesian confidence interval (shortest or central interval depending on the option)
1202
1203\param[in] total number of total events
1204\param[in] passed 0 <= number of passed events <= total
1205\param[in] level confidence level
1206\param[in] alpha shape parameter > 0 for the prior distribution (fBeta_alpha)
1207\param[in] beta shape parameter > 0 for the prior distribution (fBeta_beta)
1208\param[in] bUpper
1209 - true - upper boundary is returned
1210 - false - lower boundary is returned
1211\param[in] bShortest ??
1212
1213Note: In the case central confidence interval is calculated.
1214 when passed = 0 (or passed = total) the lower (or upper)
1215 interval values will be larger than 0 (or smaller than 1).
1216
1217Calculation:
1218
1219The posterior probability in bayesian statistics is given by:
1220\f[
1221 P(\varepsilon |k,N) \propto L(\varepsilon|k,N) \times Prior(\varepsilon)
1222\f]
1223As an efficiency can be interpreted as probability of a positive outcome of
1224a Bernoullli trial the likelihood function is given by the binomial
1225distribution:
1226\f[
1227 L(\varepsilon|k,N) = Binomial(N,k) \varepsilon ^{k} (1 - \varepsilon)^{N-k}
1228\f]
1229At the moment only beta distributions are supported as prior probabilities
1230of the efficiency (\f$ B(\alpha,\beta)\f$ is the beta function):
1231\f[
1232 Prior(\varepsilon) = \frac{1}{B(\alpha,\beta)} \varepsilon ^{\alpha - 1} (1 - \varepsilon)^{\beta - 1}
1233\f]
1234The posterior probability is therefore again given by a beta distribution:
1235\f[
1236 P(\varepsilon |k,N) \propto \varepsilon ^{k + \alpha - 1} (1 - \varepsilon)^{N - k + \beta - 1}
1237\f]
1238In case of central intervals
1239the lower boundary for the equal-tailed confidence interval is given by the
1240inverse cumulative (= quantile) function for the quantile \f$ \frac{1 - level}{2} \f$.
1241The upper boundary for the equal-tailed confidence interval is given by the
1242inverse cumulative (= quantile) function for the quantile \f$ \frac{1 + level}{2} \f$.
1243Hence it is the solution \f$ \varepsilon \f$ of the following equation:
1244\f[
1245 I_{\varepsilon}(k + \alpha,N - k + \beta) = \frac{1}{norm} \int_{0}^{\varepsilon} dt t^{k + \alpha - 1} (1 - t)^{N - k + \beta - 1} = \frac{1 \pm level}{2}
1246\f]
1247In the case of shortest interval the minimum interval around the mode is found by minimizing the length of all intervals width the
1248given probability content. See TEfficiency::BetaShortestInterval
1249*/
1250
1252{
1253 Double_t a = double(passed)+alpha;
1254 Double_t b = double(total-passed)+beta;
1255
1256 if (bShortest) {
1257 double lower = 0;
1258 double upper = 1;
1259 BetaShortestInterval(level,a,b,lower,upper);
1260 return (bUpper) ? upper : lower;
1261 }
1262 else
1263 return BetaCentralInterval(level, a, b, bUpper);
1264}
1265
1266////////////////////////////////////////////////////////////////////////////////
1267/// Calculates the boundaries for a central confidence interval for a Beta distribution
1268///
1269/// \param[in] level confidence level
1270/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1271/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1272/// \param[in] bUpper true - upper boundary is returned
1273/// false - lower boundary is returned
1274
1276{
1277 if(bUpper) {
1278 if((a > 0) && (b > 0))
1279 return ROOT::Math::beta_quantile((1+level)/2,a,b);
1280 else {
1281 gROOT->Error("TEfficiency::BayesianCentral","Invalid input parameters - return 1");
1282 return 1;
1283 }
1284 }
1285 else {
1286 if((a > 0) && (b > 0))
1287 return ROOT::Math::beta_quantile((1-level)/2,a,b);
1288 else {
1289 gROOT->Error("TEfficiency::BayesianCentral","Invalid input parameters - return 0");
1290 return 0;
1291 }
1292 }
1293}
1294
1297 fCL(level), fAlpha(alpha), fBeta(beta)
1298 {}
1299
1301 // max allowed value of lower given the interval size
1303 }
1304
1305 Double_t operator() (double lower) const {
1306 // return length of interval
1307 Double_t plow = ROOT::Math::beta_cdf(lower, fAlpha, fBeta);
1308 Double_t pup = plow + fCL;
1309 double upper = ROOT::Math::beta_quantile(pup, fAlpha,fBeta);
1310 return upper-lower;
1311 }
1312 Double_t fCL; // interval size (confidence level)
1313 Double_t fAlpha; // beta distribution alpha parameter
1314 Double_t fBeta; // beta distribution beta parameter
1315
1316};
1317
1318////////////////////////////////////////////////////////////////////////////////
1319/// Calculates the boundaries for a shortest confidence interval for a Beta distribution
1320///
1321/// \param[in] level confidence level
1322/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1323/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1324/// \param[out] upper upper boundary is returned
1325/// \param[out] lower lower boundary is returned
1326///
1327/// The lower/upper boundary are then obtained by finding the shortest interval of the beta distribution
1328/// contained the desired probability level.
1329/// The length of all possible intervals is minimized in order to find the shortest one
1330
1332{
1333 if (a <= 0 || b <= 0) {
1334 lower = 0; upper = 1;
1335 gROOT->Error("TEfficiency::BayesianShortest","Invalid input parameters - return [0,1]");
1336 return kFALSE;
1337 }
1338
1339 // treat here special cases when mode == 0 or 1
1340 double mode = BetaMode(a,b);
1341 if (mode == 0.0) {
1342 lower = 0;
1343 upper = ROOT::Math::beta_quantile(level, a, b);
1344 return kTRUE;
1345 }
1346 if (mode == 1.0) {
1347 lower = ROOT::Math::beta_quantile_c(level, a, b);
1348 upper = 1.0;
1349 return kTRUE;
1350 }
1351 // special case when the shortest interval is undefined return the central interval
1352 // can happen for a posterior when passed=total=0
1353 //
1354 if ( a==b && a<=1.0) {
1355 lower = BetaCentralInterval(level,a,b,kFALSE);
1356 upper = BetaCentralInterval(level,a,b,kTRUE);
1357 return kTRUE;
1358 }
1359
1360 // for the other case perform a minimization
1361 // make a function of the length of the posterior interval as a function of lower bound
1362 Beta_interval_length intervalLength(level,a,b);
1363 // minimize the interval length
1366 minim.SetFunction(func, 0, intervalLength.LowerMax() );
1367 minim.SetNpx(2); // no need to bracket with many iterations. Just do few times to estimate some better points
1368 bool ret = minim.Minimize(100, 1.E-10,1.E-10);
1369 if (!ret) {
1370 gROOT->Error("TEfficiency::BayesianShortes","Error finding the shortest interval");
1371 return kFALSE;
1372 }
1373 lower = minim.XMinimum();
1374 upper = lower + minim.FValMinimum();
1375 return kTRUE;
1376}
1377
1378////////////////////////////////////////////////////////////////////////////////
1379/// Compute the mean (average) of the beta distribution
1380///
1381/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1382/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1383///
1384
1386{
1387 if (a <= 0 || b <= 0 ) {
1388 gROOT->Error("TEfficiency::BayesianMean","Invalid input parameters - return 0");
1389 return 0;
1390 }
1391
1392 Double_t mean = a / (a + b);
1393 return mean;
1394}
1395
1396////////////////////////////////////////////////////////////////////////////////
1397/// Compute the mode of the beta distribution
1398///
1399/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1400/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1401///
1402/// note the mode is defined for a Beta(a,b) only if (a,b)>1 (a = passed+alpha; b = total-passed+beta)
1403/// return then the following in case (a,b) < 1:
1404/// - if (a==b) return 0.5 (it is really undefined)
1405/// - if (a < b) return 0;
1406/// - if (a > b) return 1;
1407
1409{
1410 if (a <= 0 || b <= 0 ) {
1411 gROOT->Error("TEfficiency::BayesianMode","Invalid input parameters - return 0");
1412 return 0;
1413 }
1414 if ( a <= 1 || b <= 1) {
1415 if ( a < b) return 0;
1416 if ( a > b) return 1;
1417 if (a == b) return 0.5; // cannot do otherwise
1418 }
1419
1420 // since a and b are > 1 here denominator cannot be 0 or < 0
1421 Double_t mode = (a - 1.0) / (a + b -2.0);
1422 return mode;
1423}
1424////////////////////////////////////////////////////////////////////////////////
1425/// Building standard data structure of a TEfficiency object
1426///
1427/// Notes:
1428/// - calls: SetName(name), SetTitle(title)
1429/// - set the statistic option to the default (kFCP)
1430/// - appends this object to the current directory SetDirectory(gDirectory)
1431
1432void TEfficiency::Build(const char* name,const char* title)
1433{
1434 SetName(name);
1435 SetTitle(title);
1436
1439
1440 SetBit(kPosteriorMode,false);
1442 SetBit(kUseWeights,false);
1443
1444 //set normalisation factors to 0, otherwise the += may not work properly
1447}
1448
1449////////////////////////////////////////////////////////////////////////////////
1450/// Checks binning for each axis
1451///
1452/// It is assumed that the passed histograms have the same dimension.
1453
1455{
1456
1457 const TAxis* ax1 = 0;
1458 const TAxis* ax2 = 0;
1459
1460 //check binning along axis
1461 for(Int_t j = 0; j < pass.GetDimension(); ++j) {
1462 switch(j) {
1463 case 0:
1464 ax1 = pass.GetXaxis();
1465 ax2 = total.GetXaxis();
1466 break;
1467 case 1:
1468 ax1 = pass.GetYaxis();
1469 ax2 = total.GetYaxis();
1470 break;
1471 case 2:
1472 ax1 = pass.GetZaxis();
1473 ax2 = total.GetZaxis();
1474 break;
1475 }
1476
1477 if(ax1->GetNbins() != ax2->GetNbins()) {
1478 gROOT->Info("TEfficiency::CheckBinning","Histograms are not consistent: they have different number of bins");
1479 return false;
1480 }
1481
1482 for(Int_t i = 1; i <= ax1->GetNbins() + 1; ++i)
1483 if(!TMath::AreEqualRel(ax1->GetBinLowEdge(i), ax2->GetBinLowEdge(i), 1.E-15)) {
1484 gROOT->Info("TEfficiency::CheckBinning","Histograms are not consistent: they have different bin edges");
1485 return false;
1486 }
1487
1488
1489 }
1490
1491 return true;
1492}
1493
1494////////////////////////////////////////////////////////////////////////////////
1495/// Checks the consistence of the given histograms
1496///
1497/// The histograms are considered as consistent if:
1498/// - both have the same dimension
1499/// - both have the same binning
1500/// - pass.GetBinContent(i) <= total.GetBinContent(i) for each bin i
1501///
1502
1504{
1505 if(pass.GetDimension() != total.GetDimension()) {
1506 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects have different dimensions");
1507 return false;
1508 }
1509
1510 if(!CheckBinning(pass,total)) {
1511 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects have different binning");
1512 return false;
1513 }
1514
1515 if(!CheckEntries(pass,total)) {
1516 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects do not have consistent bin contents");
1517 return false;
1518 }
1519
1520 return true;
1521}
1522
1523////////////////////////////////////////////////////////////////////////////////
1524/// Checks whether bin contents are compatible with binomial statistics
1525///
1526/// The following inequality has to be valid for each bin i:
1527/// total.GetBinContent(i) >= pass.GetBinContent(i)
1528///
1529///
1530///
1531/// Note:
1532///
1533/// - It is assumed that both histograms have the same dimension and binning.
1534
1536{
1537
1538 //check: pass <= total
1539 Int_t nbinsx, nbinsy, nbinsz, nbins;
1540
1541 nbinsx = pass.GetNbinsX();
1542 nbinsy = pass.GetNbinsY();
1543 nbinsz = pass.GetNbinsZ();
1544
1545 switch(pass.GetDimension()) {
1546 case 1: nbins = nbinsx + 2; break;
1547 case 2: nbins = (nbinsx + 2) * (nbinsy + 2); break;
1548 case 3: nbins = (nbinsx + 2) * (nbinsy + 2) * (nbinsz + 2); break;
1549 default: nbins = 0;
1550 }
1551
1552 for(Int_t i = 0; i < nbins; ++i) {
1553 if(pass.GetBinContent(i) > total.GetBinContent(i)) {
1554 gROOT->Info("TEfficiency::CheckEntries","Histograms are not consistent: passed bin content > total bin content");
1555 return false;
1556 }
1557 }
1558
1559 return true;
1560}
1561
1562////////////////////////////////////////////////////////////////////////////////
1563/// Check if both histogram are weighted. If they are weighted a true is returned
1564///
1566{
1567 if (pass.GetSumw2N() == 0 && total.GetSumw2N() == 0) return false;
1568
1569 // check also that the total sum of weight and weight squares are consistent
1570 Double_t statpass[TH1::kNstat];
1571 Double_t stattotal[TH1::kNstat];
1572
1573 pass.GetStats(statpass);
1574 total.GetStats(stattotal);
1575
1576 double tolerance = (total.IsA() == TH1F::Class() ) ? 1.E-5 : 1.E-12;
1577
1578 //require: sum of weights == sum of weights^2
1579 if(!TMath::AreEqualRel(statpass[0],statpass[1],tolerance) ||
1580 !TMath::AreEqualRel(stattotal[0],stattotal[1],tolerance) ) {
1581 return true;
1582 }
1583
1584 // histograms are not weighted
1585 return false;
1586
1587}
1588
1589
1590////////////////////////////////////////////////////////////////////////////////
1591/// Create the graph used be painted (for dim=1 TEfficiency)
1592/// The return object is managed by the caller
1593
1595{
1596 if (GetDimension() != 1) {
1597 Error("CreatePaintingGraph","Call this function only for dimension == 1");
1598 return 0;
1599 }
1600
1601
1602 Int_t npoints = fTotalHistogram->GetNbinsX();
1603 TGraphAsymmErrors * graph = new TGraphAsymmErrors(npoints);
1604 graph->SetName("eff_graph");
1605 FillGraph(graph,opt);
1606
1607 return graph;
1608}
1609
1610
1611////////////////////////////////////////////////////////////////////////////////
1612/// Fill the graph to be painted with information from TEfficiency
1613/// Internal method called by TEfficiency::Paint or TEfficiency::CreateGraph
1614
1616{
1617 TString option = opt;
1618 option.ToLower();
1619
1620 Bool_t plot0Bins = false;
1621 if (option.Contains("e0") ) plot0Bins = true;
1622
1623 Double_t x,y,xlow,xup,ylow,yup;
1624 //point i corresponds to bin i+1 in histogram
1625 // point j is point graph index
1626 // LM: cannot use TGraph::SetPoint because it deletes the underlying
1627 // histogram each time (see TGraph::SetPoint)
1628 // so use it only when extra points are added to the graph
1629 Int_t j = 0;
1630 double * px = graph->GetX();
1631 double * py = graph->GetY();
1632 double * exl = graph->GetEXlow();
1633 double * exh = graph->GetEXhigh();
1634 double * eyl = graph->GetEYlow();
1635 double * eyh = graph->GetEYhigh();
1636 Int_t npoints = fTotalHistogram->GetNbinsX();
1637 for (Int_t i = 0; i < npoints; ++i) {
1638 if (!plot0Bins && fTotalHistogram->GetBinContent(i+1) == 0 ) continue;
1640 y = GetEfficiency(i+1);
1642 xup = fTotalHistogram->GetBinWidth(i+1) - xlow;
1643 ylow = GetEfficiencyErrorLow(i+1);
1644 yup = GetEfficiencyErrorUp(i+1);
1645 // in the case the graph already existed and extra points have been added
1646 if (j >= graph->GetN() ) {
1647 graph->SetPoint(j,x,y);
1648 graph->SetPointError(j,xlow,xup,ylow,yup);
1649 }
1650 else {
1651 px[j] = x;
1652 py[j] = y;
1653 exl[j] = xlow;
1654 exh[j] = xup;
1655 eyl[j] = ylow;
1656 eyh[j] = yup;
1657 }
1658 j++;
1659 }
1660
1661 // tell the graph the effective number of points
1662 graph->Set(j);
1663 //refresh title before painting if changed
1664 TString oldTitle = graph->GetTitle();
1665 TString newTitle = GetTitle();
1666 if (oldTitle != newTitle ) {
1667 graph->SetTitle(newTitle);
1668 }
1669
1670 // set the axis labels
1673 if (xlabel) graph->GetXaxis()->SetTitle(xlabel);
1674 if (ylabel) graph->GetYaxis()->SetTitle(ylabel);
1675
1676 //copying style information
1680
1681 // copy axis labels if existing. Assume are there in the total histogram
1682 if (fTotalHistogram->GetXaxis()->GetLabels() != nullptr) {
1683 for (int ibin = 1; ibin <= fTotalHistogram->GetXaxis()->GetNbins(); ++ibin) {
1684 // we need to find the right bin for the Histogram representing the xaxis of the graph
1685 int grbin = graph->GetXaxis()->FindBin(fTotalHistogram->GetXaxis()->GetBinCenter(ibin));
1686 graph->GetXaxis()->SetBinLabel(grbin, fTotalHistogram->GetXaxis()->GetBinLabel(ibin));
1687 }
1688 }
1689 // this method forces the graph to compute correctly the axis
1690 // according to the given points
1691 graph->GetHistogram();
1692
1693}
1694
1695////////////////////////////////////////////////////////////////////////////////
1696/// Create the histogram used to be painted (for dim=2 TEfficiency)
1697/// The return object is managed by the caller
1698
1700{
1701 if (GetDimension() != 2) {
1702 Error("CreatePaintingistogram","Call this function only for dimension == 2");
1703 return 0;
1704 }
1705
1706 Int_t nbinsx = fTotalHistogram->GetNbinsX();
1707 Int_t nbinsy = fTotalHistogram->GetNbinsY();
1708 TAxis * xaxis = fTotalHistogram->GetXaxis();
1709 TAxis * yaxis = fTotalHistogram->GetYaxis();
1710 TH2 * hist = 0;
1711
1712 if (xaxis->IsVariableBinSize() && yaxis->IsVariableBinSize() )
1713 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXbins()->GetArray(),
1714 nbinsy,yaxis->GetXbins()->GetArray());
1715 else if (xaxis->IsVariableBinSize() && ! yaxis->IsVariableBinSize() )
1716 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXbins()->GetArray(),
1717 nbinsy,yaxis->GetXmin(), yaxis->GetXmax());
1718 else if (!xaxis->IsVariableBinSize() && yaxis->IsVariableBinSize() )
1719 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXmin(), xaxis->GetXmax(),
1720 nbinsy,yaxis->GetXbins()->GetArray());
1721 else
1722 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXmin(), xaxis->GetXmax(),
1723 nbinsy,yaxis->GetXmin(), yaxis->GetXmax());
1724
1725
1726 hist->SetDirectory(0);
1727
1728 FillHistogram(hist);
1729
1730 return hist;
1731}
1732
1733////////////////////////////////////////////////////////////////////////////////
1734/// Fill the 2d histogram to be painted with information from TEfficiency 2D
1735/// Internal method called by TEfficiency::Paint or TEfficiency::CreatePaintingGraph
1736
1738{
1739 //refresh title before each painting
1740 hist->SetTitle(GetTitle());
1741
1742 // set the axis labels
1746 if (xlabel) hist->GetXaxis()->SetTitle(xlabel);
1747 if (ylabel) hist->GetYaxis()->SetTitle(ylabel);
1748 if (zlabel) hist->GetZaxis()->SetTitle(zlabel);
1749
1750 Int_t bin;
1751 Int_t nbinsx = hist->GetNbinsX();
1752 Int_t nbinsy = hist->GetNbinsY();
1753 for(Int_t i = 0; i < nbinsx + 2; ++i) {
1754 for(Int_t j = 0; j < nbinsy + 2; ++j) {
1755 bin = GetGlobalBin(i,j);
1756 hist->SetBinContent(bin,GetEfficiency(bin));
1757 }
1758 }
1759
1760 // copy axis labels if existing. Assume are there in the total histogram
1761 if (fTotalHistogram->GetXaxis()->GetLabels() != nullptr) {
1762 for (int ibinx = 1; ibinx <= fTotalHistogram->GetXaxis()->GetNbins(); ++ibinx)
1763 hist->GetXaxis()->SetBinLabel(ibinx, fTotalHistogram->GetXaxis()->GetBinLabel(ibinx));
1764 }
1765 if (fTotalHistogram->GetYaxis()->GetLabels() != nullptr) {
1766 for (int ibiny = 1; ibiny <= fTotalHistogram->GetYaxis()->GetNbins(); ++ibiny)
1767 hist->GetYaxis()->SetBinLabel(ibiny, fTotalHistogram->GetYaxis()->GetBinLabel(ibiny));
1768 }
1769
1770 //copying style information
1771 TAttLine::Copy(*hist);
1772 TAttFill::Copy(*hist);
1773 TAttMarker::Copy(*hist);
1774 hist->SetStats(0);
1775
1776 return;
1777
1778}
1779////////////////////////////////////////////////////////////////////////////////
1780/**
1781Calculates the boundaries for the frequentist Clopper-Pearson interval
1782
1783This interval is recommended by the PDG.
1784
1785\param[in] total number of total events
1786\param[in] passed 0 <= number of passed events <= total
1787\param[in] level confidence level
1788\param[in] bUpper true - upper boundary is returned
1789 ;false - lower boundary is returned
1790
1791Calculation:
1792
1793The lower boundary of the Clopper-Pearson interval is the "exact" inversion
1794of the test:
1795 \f{eqnarray*}{
1796 P(x \geq passed; total) &=& \frac{1 - level}{2}\\
1797 P(x \geq passed; total) &=& 1 - P(x \leq passed - 1; total)\\
1798 &=& 1 - \frac{1}{norm} * \int_{0}^{1 - \varepsilon} t^{total - passed} (1 - t)^{passed - 1} dt\\
1799 &=& 1 - \frac{1}{norm} * \int_{\varepsilon}^{1} t^{passed - 1} (1 - t)^{total - passed} dt\\
1800 &=& \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed - 1} (1 - t)^{total - passed} dt\\
1801 &=& I_{\varepsilon}(passed,total - passed + 1)
1802 \f}
1803The lower boundary is therefore given by the \f$ \frac{1 - level}{2}\f$ quantile
1804of the beta distribution.
1805
1806The upper boundary of the Clopper-Pearson interval is the "exact" inversion
1807of the test:
1808 \f{eqnarray*}{
1809 P(x \leq passed; total) &=& \frac{1 - level}{2}\\
1810 P(x \leq passed; total) &=& \frac{1}{norm} * \int_{0}^{1 - \varepsilon} t^{total - passed - 1} (1 - t)^{passed} dt\\
1811 &=& \frac{1}{norm} * \int_{\varepsilon}^{1} t^{passed} (1 - t)^{total - passed - 1} dt\\
1812 &=& 1 - \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed} (1 - t)^{total - passed - 1} dt\\
1813 \Rightarrow 1 - \frac{1 - level}{2} &=& \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed} (1 - t)^{total - passed -1} dt\\
1814 \frac{1 + level}{2} &=& I_{\varepsilon}(passed + 1,total - passed)
1815 \f}
1816The upper boundary is therefore given by the \f$\frac{1 + level}{2}\f$ quantile
1817of the beta distribution.
1818
1819Note: The connection between the binomial distribution and the regularized
1820 incomplete beta function \f$ I_{\varepsilon}(\alpha,\beta)\f$ has been used.
1821*/
1822
1824{
1825 Double_t alpha = (1.0 - level) / 2;
1826 if(bUpper)
1827 return ((passed == total) ? 1.0 : ROOT::Math::beta_quantile(1 - alpha,passed + 1,total-passed));
1828 else
1829 return ((passed == 0) ? 0.0 : ROOT::Math::beta_quantile(alpha,passed,total-passed+1.0));
1830}
1831////////////////////////////////////////////////////////////////////////////////
1832/**
1833 Calculates the combined efficiency and its uncertainties
1834
1835 This method does a bayesian combination of the given samples.
1836
1837 \param[in] up contains the upper limit of the confidence interval afterwards
1838 \param[in] low contains the lower limit of the confidence interval afterwards
1839 \param[in] n number of samples which are combined
1840 \param[in] pass array of length n containing the number of passed events
1841 \param[in] total array of length n containing the corresponding numbers of total events
1842 \param[in] alpha shape parameters for the beta distribution as prior
1843 \param[in] beta shape parameters for the beta distribution as prior
1844 \param[in] level desired confidence level
1845 \param[in] w weights for each sample; if not given, all samples get the weight 1
1846 The weights do not need to be normalized, since they are internally renormalized
1847 to the number of effective entries.
1848 \param[in] opt
1849 - mode : The mode is returned instead of the mean of the posterior as best value
1850 When using the mode the shortest interval is also computed instead of the central one
1851 - shortest: compute shortest interval (done by default if mode option is set)
1852 - central: compute central interval (done by default if mode option is NOT set)
1853
1854 Calculation:
1855
1856 The combined posterior distributions is calculated from the Bayes theorem assuming a common prior Beta distribution.
1857 It is easy to proof that the combined posterior is then:
1858 \f{eqnarray*}{
1859 P_{comb}(\epsilon |{w_{i}}; {k_{i}}; {N_{i}}) &=& B(\epsilon, \sum_{i}{ w_{i} k_{i}} + \alpha, \sum_{i}{ w_{i}(n_{i}-k_{i})}+\beta)\\
1860 w_{i} &=& weight\ for\ each\ sample\ renormalized\ to\ the\ effective\ entries\\
1861 w^{'}_{i} &=& w_{i} \frac{ \sum_{i} {w_{i} } } { \sum_{i} {w_{i}^{2} } }
1862 \f}
1863
1864 The estimated efficiency is the mode (or the mean) of the obtained posterior distribution
1865
1866 The boundaries of the confidence interval for a confidence level (1 - a)
1867 are given by the a/2 and 1-a/2 quantiles of the resulting cumulative
1868 distribution.
1869
1870 Example (uniform prior distribution):
1871
1872Begin_Macro(source)
1873{
1874 TCanvas* c1 = new TCanvas("c1","",600,800);
1875 c1->Divide(1,2);
1876 c1->SetFillStyle(1001);
1877 c1->SetFillColor(kWhite);
1878
1879 TF1* p1 = new TF1("p1","TMath::BetaDist(x,19,9)",0,1);
1880 TF1* p2 = new TF1("p2","TMath::BetaDist(x,4,8)",0,1);
1881 TF1* comb = new TF1("comb2","TMath::BetaDist(x,[0],[1])",0,1);
1882 double nrm = 1./(0.6*0.6+0.4*0.4); // weight normalization
1883 double a = 0.6*18.0 + 0.4*3.0 + 1.0; // new alpha parameter of combined beta dist.
1884 double b = 0.6*10+0.4*7+1.0; // new beta parameter of combined beta dist.
1885 comb->SetParameters(nrm*a ,nrm *b );
1886 TF1* const1 = new TF1("const1","0.05",0,1);
1887 TF1* const2 = new TF1("const2","0.95",0,1);
1888
1889 p1->SetLineColor(kRed);
1890 p1->SetTitle("combined posteriors;#epsilon;P(#epsilon|k,N)");
1891 p2->SetLineColor(kBlue);
1892 comb->SetLineColor(kGreen+2);
1893
1894 TLegend* leg1 = new TLegend(0.12,0.65,0.5,0.85);
1895 leg1->AddEntry(p1,"k1 = 18, N1 = 26","l");
1896 leg1->AddEntry(p2,"k2 = 3, N2 = 10","l");
1897 leg1->AddEntry(comb,"combined: p1 = 0.6, p2=0.4","l");
1898
1899 c1->cd(1);
1900 comb->Draw();
1901 p1->Draw("same");
1902 p2->Draw("same");
1903 leg1->Draw("same");
1904 c1->cd(2);
1905 const1->SetLineWidth(1);
1906 const2->SetLineWidth(1);
1907 TGraph* gr = (TGraph*)comb->DrawIntegral();
1908 gr->SetTitle("cumulative function of combined posterior with boundaries for cl = 95%;#epsilon;CDF");
1909 const1->Draw("same");
1910 const2->Draw("same");
1911
1912 c1->cd(0);
1913 return c1;
1914}
1915End_Macro
1916
1917**/
1918////////////////////////////////////////////////////////////////////
1920 const Int_t* pass,const Int_t* total,
1921 Double_t alpha, Double_t beta,
1922 Double_t level,const Double_t* w,Option_t* opt)
1923{
1924 TString option(opt);
1925 option.ToLower();
1926
1927 //LM: new formula for combination
1928 // works only if alpha beta are the same always
1929 // the weights are normalized to w(i) -> N_eff w(i)/ Sum w(i)
1930 // i.e. w(i) -> Sum (w(i) / Sum (w(i)^2) * w(i)
1931 // norm = Sum (w(i) / Sum (w(i)^2)
1932 double ntot = 0;
1933 double ktot = 0;
1934 double sumw = 0;
1935 double sumw2 = 0;
1936 for (int i = 0; i < n ; ++i) {
1937 if(pass[i] > total[i]) {
1938 ::Error("TEfficiency::Combine","total events = %i < passed events %i",total[i],pass[i]);
1939 ::Info("TEfficiency::Combine","stop combining");
1940 return -1;
1941 }
1942
1943 ntot += w[i] * total[i];
1944 ktot += w[i] * pass[i];
1945 sumw += w[i];
1946 sumw2 += w[i]*w[i];
1947 //mean += w[i] * (pass[i] + alpha[i])/(total[i] + alpha[i] + beta[i]);
1948 }
1949 double norm = sumw/sumw2;
1950 ntot *= norm;
1951 ktot *= norm;
1952 if(ktot > ntot) {
1953 ::Error("TEfficiency::Combine","total = %f < passed %f",ntot,ktot);
1954 ::Info("TEfficiency::Combine","stop combining");
1955 return -1;
1956 }
1957
1958 double a = ktot + alpha;
1959 double b = ntot - ktot + beta;
1960
1961 double mean = a/(a+b);
1962 double mode = BetaMode(a,b);
1963
1964
1965 Bool_t shortestInterval = option.Contains("sh") || ( option.Contains("mode") && !option.Contains("cent") );
1966
1967 if (shortestInterval)
1968 BetaShortestInterval(level, a, b, low, up);
1969 else {
1970 low = BetaCentralInterval(level, a, b, false);
1971 up = BetaCentralInterval(level, a, b, true);
1972 }
1973
1974 if (option.Contains("mode")) return mode;
1975 return mean;
1976
1977}
1978////////////////////////////////////////////////////////////////////////////////
1979/// Combines a list of 1-dimensional TEfficiency objects
1980///
1981/// A TGraphAsymmErrors object is returned which contains the estimated
1982/// efficiency and its uncertainty for each bin.
1983/// If the combination fails, a zero pointer is returned.
1984///
1985/// At the moment the combining is only implemented for bayesian statistics.
1986///
1987/// \param[in] pList list containing TEfficiency objects which should be combined
1988/// only one-dimensional efficiencies are taken into account
1989/// \param[in] option
1990/// - s : strict combining; only TEfficiency objects with the same beta
1991/// prior and the flag kIsBayesian == true are combined
1992/// If not specified the prior parameter of the first TEfficiency object is used
1993/// - v : verbose mode; print information about combining
1994/// - cl=x : set confidence level (0 < cl < 1). If not specified, the
1995/// confidence level of the first TEfficiency object is used.
1996/// - mode Use mode of combined posterior as estimated value for the efficiency
1997/// - shortest: compute shortest interval (done by default if mode option is set)
1998/// - central: compute central interval (done by default if mode option is NOT set)
1999/// \param[in] n number of weights (has to be the number of one-dimensional
2000/// TEfficiency objects in pList)
2001/// If no weights are passed, the internal weights GetWeight() of
2002/// the given TEfficiency objects are used.
2003/// \param[in] w array of length n with weights for each TEfficiency object in
2004/// pList (w[0] correspond to pList->First ... w[n-1] -> pList->Last)
2005/// The weights do not have to be normalised.
2006///
2007/// For each bin the calculation is done by the Combine(double&, double& ...) method.
2008
2010 Int_t n,const Double_t* w)
2011{
2012 TString opt = option;
2013 opt.ToLower();
2014
2015 //parameter of prior distribution, confidence level and normalisation factor
2016 Double_t alpha = -1;
2017 Double_t beta = -1;
2018 Double_t level = 0;
2019
2020 //flags for combining
2021 Bool_t bStrict = false;
2022 Bool_t bOutput = false;
2023 Bool_t bWeights = false;
2024 //list of all information needed to weight and combine efficiencies
2025 std::vector<TH1*> vTotal; vTotal.reserve(n);
2026 std::vector<TH1*> vPassed; vPassed.reserve(n);
2027 std::vector<Double_t> vWeights; vWeights.reserve(n);
2028 // std::vector<Double_t> vAlpha;
2029 // std::vector<Double_t> vBeta;
2030
2031 if(opt.Contains("s")) {
2032 opt.ReplaceAll("s","");
2033 bStrict = true;
2034 }
2035
2036 if(opt.Contains("v")) {
2037 opt.ReplaceAll("v","");
2038 bOutput = true;
2039 }
2040
2041 if(opt.Contains("cl=")) {
2042 Ssiz_t pos = opt.Index("cl=") + 3;
2043 level = atof( opt(pos,opt.Length() ).Data() );
2044 if((level <= 0) || (level >= 1))
2045 level = 0;
2046 opt.ReplaceAll("cl=","");
2047 }
2048
2049 //are weights explicitly given
2050 if(n && w) {
2051 bWeights = true;
2052 for(Int_t k = 0; k < n; ++k) {
2053 if(w[k] > 0)
2054 vWeights.push_back(w[k]);
2055 else {
2056 gROOT->Error("TEfficiency::Combine","invalid custom weight found w = %.2lf",w[k]);
2057 gROOT->Info("TEfficiency::Combine","stop combining");
2058 return 0;
2059 }
2060 }
2061 }
2062
2063 TIter next(pList);
2064 TObject* obj = 0;
2065 TEfficiency* pEff = 0;
2066 while((obj = next())) {
2067 pEff = dynamic_cast<TEfficiency*>(obj);
2068 //is object a TEfficiency object?
2069 if(pEff) {
2070 if(pEff->GetDimension() > 1)
2071 continue;
2072 if(!level) level = pEff->GetConfidenceLevel();
2073
2074 if(alpha<1) alpha = pEff->GetBetaAlpha();
2075 if(beta<1) beta = pEff->GetBetaBeta();
2076
2077 //if strict combining, check priors, confidence level and statistic
2078 if(bStrict) {
2079 if(alpha != pEff->GetBetaAlpha())
2080 continue;
2081 if(beta != pEff->GetBetaBeta())
2082 continue;
2083 if(!pEff->UsesBayesianStat())
2084 continue;
2085 }
2086
2087 vTotal.push_back(pEff->fTotalHistogram);
2088 vPassed.push_back(pEff->fPassedHistogram);
2089
2090 //no weights given -> use weights of TEfficiency objects
2091 if(!bWeights)
2092 vWeights.push_back(pEff->fWeight);
2093
2094 //strict combining -> using global prior
2095 // if(bStrict) {
2096 // vAlpha.push_back(alpha);
2097 // vBeta.push_back(beta);
2098 // }
2099 // else {
2100 // vAlpha.push_back(pEff->GetBetaAlpha());
2101 // vBeta.push_back(pEff->GetBetaBeta());
2102 // }
2103 }
2104 }
2105
2106 //no TEfficiency objects found
2107 if(vTotal.empty()) {
2108 gROOT->Error("TEfficiency::Combine","no TEfficiency objects in given list");
2109 gROOT->Info("TEfficiency::Combine","stop combining");
2110 return 0;
2111 }
2112
2113 //invalid number of custom weights
2114 if(bWeights && (n != (Int_t)vTotal.size())) {
2115 gROOT->Error("TEfficiency::Combine","number of weights n=%i differs from number of TEfficiency objects k=%i which should be combined",n,(Int_t)vTotal.size());
2116 gROOT->Info("TEfficiency::Combine","stop combining");
2117 return 0;
2118 }
2119
2120 Int_t nbins_max = vTotal.at(0)->GetNbinsX();
2121 //check binning of all histograms
2122 for(UInt_t i=0; i<vTotal.size(); ++i) {
2123 if (!TEfficiency::CheckBinning(*vTotal.at(0),*vTotal.at(i)) )
2124 gROOT->Warning("TEfficiency::Combine","histograms have not the same binning -> results may be useless");
2125 if(vTotal.at(i)->GetNbinsX() < nbins_max) nbins_max = vTotal.at(i)->GetNbinsX();
2126 }
2127
2128 //display information about combining
2129 if(bOutput) {
2130 gROOT->Info("TEfficiency::Combine","combining %i TEfficiency objects",(Int_t)vTotal.size());
2131 if(bWeights)
2132 gROOT->Info("TEfficiency::Combine","using custom weights");
2133 if(bStrict) {
2134 gROOT->Info("TEfficiency::Combine","using the following prior probability for the efficiency: P(e) ~ Beta(e,%.3lf,%.3lf)",alpha,beta);
2135 }
2136 else
2137 gROOT->Info("TEfficiency::Combine","using individual priors of each TEfficiency object");
2138 gROOT->Info("TEfficiency::Combine","confidence level = %.2lf",level);
2139 }
2140
2141 //create TGraphAsymmErrors with efficiency
2142 std::vector<Double_t> x(nbins_max);
2143 std::vector<Double_t> xlow(nbins_max);
2144 std::vector<Double_t> xhigh(nbins_max);
2145 std::vector<Double_t> eff(nbins_max);
2146 std::vector<Double_t> efflow(nbins_max);
2147 std::vector<Double_t> effhigh(nbins_max);
2148
2149 //parameters for combining:
2150 //number of objects
2151 Int_t num = vTotal.size();
2152 std::vector<Int_t> pass(num);
2153 std::vector<Int_t> total(num);
2154
2155 //loop over all bins
2156 Double_t low = 0;
2157 Double_t up = 0;
2158 for(Int_t i=1; i <= nbins_max; ++i) {
2159 //the binning of the x-axis is taken from the first total histogram
2160 x[i-1] = vTotal.at(0)->GetBinCenter(i);
2161 xlow[i-1] = x[i-1] - vTotal.at(0)->GetBinLowEdge(i);
2162 xhigh[i-1] = vTotal.at(0)->GetBinWidth(i) - xlow[i-1];
2163
2164 for(Int_t j = 0; j < num; ++j) {
2165 pass[j] = (Int_t)(vPassed.at(j)->GetBinContent(i) + 0.5);
2166 total[j] = (Int_t)(vTotal.at(j)->GetBinContent(i) + 0.5);
2167 }
2168
2169 //fill efficiency and errors
2170 eff[i-1] = Combine(up,low,num,&pass[0],&total[0],alpha,beta,level,&vWeights[0],opt.Data());
2171 //did an error occurred ?
2172 if(eff[i-1] == -1) {
2173 gROOT->Error("TEfficiency::Combine","error occurred during combining");
2174 gROOT->Info("TEfficiency::Combine","stop combining");
2175 return 0;
2176 }
2177 efflow[i-1]= eff[i-1] - low;
2178 effhigh[i-1]= up - eff[i-1];
2179 }//loop over all bins
2180
2181 TGraphAsymmErrors* gr = new TGraphAsymmErrors(nbins_max,&x[0],&eff[0],&xlow[0],&xhigh[0],&efflow[0],&effhigh[0]);
2182
2183 return gr;
2184}
2185
2186////////////////////////////////////////////////////////////////////////////////
2187/// Compute distance from point px,py to a graph.
2188///
2189/// Compute the closest distance of approach from point px,py to this line.
2190/// The distance is computed in pixels units.
2191///
2192/// Forward the call to the painted graph
2193
2195{
2196 if (fPaintGraph) return fPaintGraph->DistancetoPrimitive(px,py);
2197 if (fPaintHisto) return fPaintHisto->DistancetoPrimitive(px,py);
2198 return 0;
2199}
2200
2201
2202////////////////////////////////////////////////////////////////////////////////
2203/// Draws the current TEfficiency object
2204///
2205/// \param[in] opt
2206/// - 1-dimensional case: same options as TGraphAsymmErrors::Draw()
2207/// but as default "AP" is used
2208/// - 2-dimensional case: same options as TH2::Draw()
2209/// - 3-dimensional case: not yet supported
2210///
2211/// Specific TEfficiency drawing options:
2212/// - E0 - plot bins where the total number of passed events is zero
2213/// (the error interval will be [0,1] )
2214
2216{
2217 //check options
2218 TString option = opt;
2219 option.ToLower();
2220 // use by default "AP"
2221 if (option.IsNull() ) option = "ap";
2222
2223 if(gPad && !option.Contains("same"))
2224 gPad->Clear();
2225 else {
2226 // add always "a" if not present
2227 if (!option.Contains("a") ) option += "a";
2228 }
2229
2230 // add always p to the option
2231 if (!option.Contains("p") ) option += "p";
2232
2233
2234 AppendPad(option.Data());
2235}
2236
2237////////////////////////////////////////////////////////////////////////////////
2238/// Execute action corresponding to one event.
2239///
2240/// This member function is called when the drawn class is clicked with the locator
2241/// If Left button clicked on one of the line end points, this point
2242/// follows the cursor until button is released.
2243///
2244/// if Middle button clicked, the line is moved parallel to itself
2245/// until the button is released.
2246/// Forward the call to the underlying graph
2247
2249{
2251 else if (fPaintHisto) fPaintHisto->ExecuteEvent(event,px,py);
2252}
2253
2254////////////////////////////////////////////////////////////////////////////////
2255/// This function is used for filling the two histograms.
2256///
2257/// \param[in] bPassed flag whether the current event passed the selection
2258/// - true: both histograms are filled
2259/// - false: only the total histogram is filled
2260/// \param[in] x x-value
2261/// \param[in] y y-value (use default=0 for 1-D efficiencies)
2262/// \param[in] z z-value (use default=0 for 2-D or 1-D efficiencies)
2263
2265{
2266 switch(GetDimension()) {
2267 case 1:
2269 if(bPassed)
2271 break;
2272 case 2:
2273 ((TH2*)(fTotalHistogram))->Fill(x,y);
2274 if(bPassed)
2275 ((TH2*)(fPassedHistogram))->Fill(x,y);
2276 break;
2277 case 3:
2278 ((TH3*)(fTotalHistogram))->Fill(x,y,z);
2279 if(bPassed)
2280 ((TH3*)(fPassedHistogram))->Fill(x,y,z);
2281 break;
2282 }
2283}
2284
2285////////////////////////////////////////////////////////////////////////////////
2286///This function is used for filling the two histograms with a weight.
2287///
2288/// \param[in] bPassed flag whether the current event passed the selection
2289/// - true: both histograms are filled
2290/// - false: only the total histogram is filled
2291/// \param[in] weight weight for the event
2292/// \param[in] x x-value
2293/// \param[in] y y-value (use default=0 for 1-D efficiencies)
2294/// \param[in] z z-value (use default=0 for 2-D or 1-D efficiencies)
2295///
2296/// Note: - this function will call SetUseWeightedEvents if it was not called by the user before
2297
2299{
2300 if(!TestBit(kUseWeights))
2301 {
2302 // Info("FillWeighted","call SetUseWeightedEvents() manually to ensure correct storage of sum of weights squared");
2304 }
2305
2306 switch(GetDimension()) {
2307 case 1:
2308 fTotalHistogram->Fill(x,weight);
2309 if(bPassed)
2310 fPassedHistogram->Fill(x,weight);
2311 break;
2312 case 2:
2313 ((TH2*)(fTotalHistogram))->Fill(x,y,weight);
2314 if(bPassed)
2315 ((TH2*)(fPassedHistogram))->Fill(x,y,weight);
2316 break;
2317 case 3:
2318 ((TH3*)(fTotalHistogram))->Fill(x,y,z,weight);
2319 if(bPassed)
2320 ((TH3*)(fPassedHistogram))->Fill(x,y,z,weight);
2321 break;
2322 }
2323}
2324
2325////////////////////////////////////////////////////////////////////////////////
2326/// Returns the global bin number containing the given values
2327///
2328/// Note:
2329///
2330/// - values which belong to dimensions higher than the current dimension
2331/// of the TEfficiency object are ignored (i.e. for 1-dimensional
2332/// efficiencies only the x-value is considered)
2333
2335{
2337 Int_t ny = 0;
2338 Int_t nz = 0;
2339
2340 switch(GetDimension()) {
2341 case 3: nz = fTotalHistogram->GetZaxis()->FindFixBin(z);
2342 case 2: ny = fTotalHistogram->GetYaxis()->FindFixBin(y);break;
2343 }
2344
2345 return GetGlobalBin(nx,ny,nz);
2346}
2347
2348////////////////////////////////////////////////////////////////////////////////
2349/// Fits the efficiency using the TBinomialEfficiencyFitter class
2350///
2351/// The resulting fit function is added to the list of associated functions.
2352///
2353/// Options:
2354/// - "+": previous fitted functions in the list are kept, by default
2355/// all functions in the list are deleted
2356/// - for more fitting options see TBinomialEfficiencyFitter::Fit
2357
2359{
2360 TString option = opt;
2361 option.ToLower();
2362
2363 //replace existing functions in list with same name
2364 Bool_t bDeleteOld = true;
2365 if(option.Contains("+")) {
2366 option.ReplaceAll("+","");
2367 bDeleteOld = false;
2368 }
2369
2371
2372 TFitResultPtr result = Fitter.Fit(f1,option.Data());
2373
2374 //create copy which is appended to the list
2375 TF1* pFunc = new TF1(*f1);
2376
2377 if(bDeleteOld) {
2378 TIter next(fFunctions);
2379 TObject* obj = 0;
2380 while((obj = next())) {
2381 if(obj->InheritsFrom(TF1::Class())) {
2382 fFunctions->Remove(obj);
2383 delete obj;
2384 }
2385 }
2386 }
2387
2388 // create list if necessary
2389 if(!fFunctions)
2390 fFunctions = new TList();
2391
2392 fFunctions->Add(pFunc);
2393
2394 return result;
2395}
2396
2397////////////////////////////////////////////////////////////////////////////////
2398/// Returns a cloned version of fPassedHistogram
2399///
2400/// Notes:
2401/// - The histogram is filled with unit weights. You might want to scale
2402/// it with the global weight GetWeight().
2403/// - The returned object is owned by the user who has to care about the
2404/// deletion of the new TH1 object.
2405/// - This histogram is by default NOT attached to the current directory
2406/// to avoid duplication of data. If you want to store it automatically
2407/// during the next TFile::Write() command, you have to attach it to
2408/// the corresponding directory.
2409///
2410/// ~~~~~~~{.cpp}
2411/// TFile* pFile = new TFile("passed.root","update");
2412/// TEfficiency* pEff = (TEfficiency*)gDirectory->Get("my_eff");
2413/// TH1* copy = pEff->GetCopyPassedHisto();
2414/// copy->SetDirectory(gDirectory);
2415/// pFile->Write();
2416/// ~~~~~~~
2417
2419{
2420 // do not add cloned histogram to gDirectory
2421 TDirectory::TContext ctx(nullptr);
2422 TH1* tmp = (TH1*)(fPassedHistogram->Clone());
2423
2424 return tmp;
2425}
2426
2427////////////////////////////////////////////////////////////////////////////////
2428/// Returns a cloned version of fTotalHistogram
2429///
2430/// Notes:
2431/// - The histogram is filled with unit weights. You might want to scale
2432/// it with the global weight GetWeight().
2433/// - The returned object is owned by the user who has to care about the
2434/// deletion of the new TH1 object.
2435/// - This histogram is by default NOT attached to the current directory
2436/// to avoid duplication of data. If you want to store it automatically
2437/// during the next TFile::Write() command, you have to attach it to
2438/// the corresponding directory.
2439///
2440/// ~~~~~~~{.cpp}
2441/// TFile* pFile = new TFile("total.root","update");
2442/// TEfficiency* pEff = (TEfficiency*)gDirectory->Get("my_eff");
2443/// TH1* copy = pEff->GetCopyTotalHisto();
2444/// copy->SetDirectory(gDirectory);
2445/// pFile->Write();
2446/// ~~~~~~~
2447
2449{
2450 // do not add cloned histogram to gDirectory
2451 TDirectory::TContext ctx(nullptr);
2452 TH1* tmp = (TH1*)(fTotalHistogram->Clone());
2453
2454 return tmp;
2455}
2456
2457////////////////////////////////////////////////////////////////////////////////
2458///returns the dimension of the current TEfficiency object
2459
2461{
2462 return fTotalHistogram->GetDimension();
2463}
2464
2465////////////////////////////////////////////////////////////////////////////////
2466/// Returns the efficiency in the given global bin
2467///
2468/// Note:
2469/// - The estimated efficiency depends on the chosen statistic option:
2470/// for frequentist ones:
2471/// \f$ \hat{\varepsilon} = \frac{passed}{total} \f$
2472/// for bayesian ones the expectation value of the resulting posterior
2473/// distribution is returned:
2474/// \f$ \hat{\varepsilon} = \frac{passed + \alpha}{total + \alpha + \beta} \f$
2475/// If the bit kPosteriorMode is set (or the method TEfficiency::UsePosteriorMode() has been called ) the
2476/// mode (most probable value) of the posterior is returned:
2477/// \f$ \hat{\varepsilon} = \frac{passed + \alpha -1}{total + \alpha + \beta -2} \f$
2478/// - If the denominator is equal to 0, an efficiency of 0 is returned.
2479/// - When \f$ passed + \alpha < 1 \f$ or \f$ total - passed + \beta < 1 \f$ the above
2480/// formula for the mode is not valid. In these cases values the estimated efficiency is 0 or 1.
2481
2483{
2486
2487 if(TestBit(kIsBayesian)) {
2488
2489 // parameters for the beta prior distribution
2492
2493 Double_t aa,bb;
2494 if(TestBit(kUseWeights))
2495 {
2497 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2499
2500 if (tw2 <= 0 ) return pw/tw;
2501
2502 // tw/tw2 renormalize the weights
2503 double norm = tw/tw2;
2504 aa = pw * norm + alpha;
2505 bb = (tw - pw) * norm + beta;
2506 }
2507 else
2508 {
2509 aa = passed + alpha;
2510 bb = total - passed + beta;
2511 }
2512
2513 if (!TestBit(kPosteriorMode) )
2514 return BetaMean(aa,bb);
2515 else
2516 return BetaMode(aa,bb);
2517
2518 }
2519 else
2520 return (total)? ((Double_t)passed)/total : 0;
2521}
2522
2523////////////////////////////////////////////////////////////////////////////////
2524/// Returns the lower error on the efficiency in the given global bin
2525///
2526/// The result depends on the current confidence level fConfLevel and the
2527/// chosen statistic option fStatisticOption. See SetStatisticOption(Int_t) for
2528/// more details.
2529///
2530/// Note: If the histograms are filled with weights, only bayesian methods and the
2531/// normal approximation are supported.
2532
2534{
2537
2538 Double_t eff = GetEfficiency(bin);
2539
2540 // check whether weights have been used
2541 if(TestBit(kUseWeights))
2542 {
2544 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2546 Double_t pw2 = fPassedHistogram->GetSumw2()->At(bin);
2547
2548 if(TestBit(kIsBayesian))
2549 {
2552
2553 if (tw2 <= 0) return 0;
2554
2555 // tw/tw2 renormalize the weights
2556 Double_t norm = tw/tw2;
2557 Double_t aa = pw * norm + alpha;
2558 Double_t bb = (tw - pw) * norm + beta;
2559 Double_t low = 0;
2560 Double_t upper = 1;
2563 }
2564 else {
2566 }
2567
2568 return eff - low;
2569 }
2570 else
2571 {
2573 {
2574 Warning("GetEfficiencyErrorLow","frequentist confidence intervals for weights are only supported by the normal approximation");
2575 Info("GetEfficiencyErrorLow","setting statistic option to kFNormal");
2576 const_cast<TEfficiency*>(this)->SetStatisticOption(kFNormal);
2577 }
2578
2579 Double_t variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
2580 Double_t sigma = sqrt(variance);
2581
2582 Double_t prob = 0.5 * (1.- fConfLevel);
2584
2585 // avoid to return errors which makes eff-err < 0
2586 return (eff - delta < 0) ? eff : delta;
2587 }
2588 }
2589 else
2590 {
2591 if(TestBit(kIsBayesian))
2592 {
2593 // parameters for the beta prior distribution
2596 return (eff - Bayesian(total,passed,fConfLevel,alpha,beta,false,TestBit(kShortestInterval)));
2597 }
2598 else
2599 return (eff - fBoundary(total,passed,fConfLevel,false));
2600 }
2601}
2602
2603////////////////////////////////////////////////////////////////////////////////
2604/// Returns the upper error on the efficiency in the given global bin
2605///
2606/// The result depends on the current confidence level fConfLevel and the
2607/// chosen statistic option fStatisticOption. See SetStatisticOption(Int_t) for
2608/// more details.
2609///
2610/// Note: If the histograms are filled with weights, only bayesian methods and the
2611/// normal approximation are supported.
2612
2614{
2617
2618 Double_t eff = GetEfficiency(bin);
2619
2620 // check whether weights have been used
2621 if(TestBit(kUseWeights))
2622 {
2624 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2626 Double_t pw2 = fPassedHistogram->GetSumw2()->At(bin);
2627
2628 if(TestBit(kIsBayesian))
2629 {
2632
2633 if (tw2 <= 0) return 0;
2634
2635 // tw/tw2 renormalize the weights
2636 Double_t norm = tw/tw2;
2637 Double_t aa = pw * norm + alpha;
2638 Double_t bb = (tw - pw) * norm + beta;
2639 Double_t low = 0;
2640 Double_t upper = 1;
2643 }
2644 else {
2645 upper = TEfficiency::BetaCentralInterval(fConfLevel,aa,bb,true);
2646 }
2647
2648 return upper - eff;
2649 }
2650 else
2651 {
2653 {
2654 Warning("GetEfficiencyErrorUp","frequentist confidence intervals for weights are only supported by the normal approximation");
2655 Info("GetEfficiencyErrorUp","setting statistic option to kFNormal");
2656 const_cast<TEfficiency*>(this)->SetStatisticOption(kFNormal);
2657 }
2658
2659 Double_t variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
2660 Double_t sigma = sqrt(variance);
2661
2662 Double_t prob = 0.5 * (1.- fConfLevel);
2664
2665 return (eff + delta > 1) ? 1.-eff : delta;
2666 }
2667 }
2668 else
2669 {
2670 if(TestBit(kIsBayesian))
2671 {
2672 // parameters for the beta prior distribution
2675 return (Bayesian(total,passed,fConfLevel,alpha,beta,true,TestBit(kShortestInterval)) - eff);
2676 }
2677 else
2678 return fBoundary(total,passed,fConfLevel,true) - eff;
2679 }
2680}
2681
2682////////////////////////////////////////////////////////////////////////////////
2683/// Returns the global bin number which can be used as argument for the
2684/// following functions:
2685///
2686/// - GetEfficiency(bin), GetEfficiencyErrorLow(bin), GetEfficiencyErrorUp(bin)
2687/// - SetPassedEvents(bin), SetTotalEvents(bin)
2688///
2689/// see TH1::GetBin() for conventions on numbering bins
2690
2692{
2693 return fTotalHistogram->GetBin(binx,biny,binz);
2694}
2695
2696////////////////////////////////////////////////////////////////////////////////
2697
2699{
2700 return (fFunctions) ? fFunctions : fFunctions = new TList();
2701}
2702
2703////////////////////////////////////////////////////////////////////////////////
2704/// Merges the TEfficiency objects in the given list to the given
2705/// TEfficiency object using the operator+=(TEfficiency&)
2706///
2707/// The merged result is stored in the current object. The statistic options and
2708/// the confidence level are taken from the current object.
2709///
2710/// This function should be used when all TEfficiency objects correspond to
2711/// the same process.
2712///
2713/// The new weight is set according to:
2714/// \f$ \frac{1}{w_{new}} = \sum_{i} \frac{1}{w_{i}} \f$
2715
2717{
2718 if(!pList->IsEmpty()) {
2719 TIter next(pList);
2720 TObject* obj = 0;
2721 TEfficiency* pEff = 0;
2722 while((obj = next())) {
2723 pEff = dynamic_cast<TEfficiency*>(obj);
2724 if(pEff) {
2725 *this += *pEff;
2726 }
2727 }
2728 }
2730}
2731
2732////////////////////////////////////////////////////////////////////////////////
2733/**
2734Returns the confidence limits for the efficiency supposing that the
2735efficiency follows a normal distribution with the rms below
2736
2737\param[in] total number of total events
2738\param[in] passed 0 <= number of passed events <= total
2739\param[in] level confidence level
2740\param[in] bUpper
2741 - true - upper boundary is returned
2742 - false - lower boundary is returned
2743
2744Calculation:
2745
2746\f{eqnarray*}{
2747 \hat{\varepsilon} &=& \frac{passed}{total}\\
2748 \sigma_{\varepsilon} &=& \sqrt{\frac{\hat{\varepsilon} (1 - \hat{\varepsilon})}{total}}\\
2749 \varepsilon_{low} &=& \hat{\varepsilon} \pm \Phi^{-1}(\frac{level}{2},\sigma_{\varepsilon})
2750\f}
2751*/
2752
2754{
2755 Double_t alpha = (1.0 - level)/2;
2756 if (total == 0) return (bUpper) ? 1 : 0;
2757 Double_t average = passed / total;
2758 Double_t sigma = std::sqrt(average * (1 - average) / total);
2759 Double_t delta = ROOT::Math::normal_quantile(1 - alpha,sigma);
2760
2761 if(bUpper)
2762 return ((average + delta) > 1) ? 1.0 : (average + delta);
2763 else
2764 return ((average - delta) < 0) ? 0.0 : (average - delta);
2765}
2766
2767////////////////////////////////////////////////////////////////////////////////
2768/// Adds the histograms of another TEfficiency object to current histograms
2769///
2770/// The statistic options and the confidence level remain unchanged.
2771///
2772/// fTotalHistogram += rhs.fTotalHistogram;
2773/// fPassedHistogram += rhs.fPassedHistogram;
2774///
2775/// calculates a new weight:
2776/// current weight of this TEfficiency object = \f$ w_{1} \f$
2777/// weight of rhs = \f$ w_{2} \f$
2778/// \f$ w_{new} = \frac{w_{1} \times w_{2}}{w_{1} + w_{2}} \f$
2779
2781{
2782
2783 if (fTotalHistogram == 0 && fPassedHistogram == 0) {
2784 // efficiency is empty just copy it over
2785 *this = rhs;
2786 return *this;
2787 }
2788 else if (fTotalHistogram == 0 || fPassedHistogram == 0) {
2789 Fatal("operator+=","Adding to a non consistent TEfficiency object which has not a total or a passed histogram ");
2790 return *this;
2791 }
2792
2793 if (rhs.fTotalHistogram == 0 && rhs.fPassedHistogram == 0 ) {
2794 Warning("operator+=","no operation: adding an empty object");
2795 return *this;
2796 }
2797 else if (rhs.fTotalHistogram == 0 || rhs.fPassedHistogram == 0 ) {
2798 Fatal("operator+=","Adding a non consistent TEfficiency object which has not a total or a passed histogram ");
2799 return *this;
2800 }
2801
2804
2807
2808 SetWeight((fWeight * rhs.GetWeight())/(fWeight + rhs.GetWeight()));
2809
2810 return *this;
2811}
2812
2813////////////////////////////////////////////////////////////////////////////////
2814/// Assignment operator
2815///
2816/// The histograms, statistic option, confidence level, weight and paint styles
2817/// of rhs are copied to the this TEfficiency object.
2818///
2819/// Note: - The list of associated functions is not copied. After this
2820/// operation the list of associated functions is empty.
2821
2823{
2824 if(this != &rhs)
2825 {
2826 //statistic options
2830 SetBetaBeta(rhs.GetBetaBeta());
2831 SetWeight(rhs.GetWeight());
2832
2833 //associated list of functions
2834 if(fFunctions)
2835 fFunctions->Delete();
2836
2837 //copy histograms
2838 delete fTotalHistogram;
2839 delete fPassedHistogram;
2840
2841 // do not add cloned histogram to gDirectory
2842 {
2843 TDirectory::TContext ctx(nullptr);
2846 }
2847 //delete temporary paint objects
2848 delete fPaintHisto;
2849 delete fPaintGraph;
2850 fPaintHisto = 0;
2851 fPaintGraph = 0;
2852
2853 //copy style
2854 rhs.TAttLine::Copy(*this);
2855 rhs.TAttFill::Copy(*this);
2856 rhs.TAttMarker::Copy(*this);
2857 }
2858
2859 return *this;
2860}
2861
2862////////////////////////////////////////////////////////////////////////////////
2863/// Paints this TEfficiency object
2864///
2865/// For details on the possible option see Draw(Option_t*)
2866///
2867/// Note for 1D classes
2868/// In 1D the TEfficiency uses a TGraphAsymmErrors for drawing
2869/// The TGraph is created only the first time Paint is used. The user can manipulate the
2870/// TGraph via the method TEfficiency::GetPaintedGraph()
2871/// The TGraph creates behing an histogram for the axis. The histogram is created also only the first time.
2872/// If the axis needs to be updated because in the meantime the class changed use this trick
2873/// which will trigger a re-calculation of the axis of the graph
2874/// TEfficiency::GetPaintedGraph()->Set(0)
2875///
2876/// Note that in order to access the painted graph via GetPaintedGraph() you need either to call Paint or better
2877/// gPad->Update();
2878///
2879
2881{
2882
2883
2884 if(!gPad)
2885 return;
2886
2887
2888 //use TGraphAsymmErrors for painting
2889 if(GetDimension() == 1) {
2890 if(!fPaintGraph) {
2891 fPaintGraph = CreateGraph(opt);
2892 }
2893 else
2894 // update existing graph already created
2895 FillGraph(fPaintGraph, opt);
2896
2897 //paint graph
2898
2899 fPaintGraph->Paint(opt);
2900
2901 //paint all associated functions
2902 if(fFunctions) {
2903 //paint box with fit parameters
2904 //the fit statistics will be painted if gStyle->SetOptFit(1) has been
2905 // called by the user
2906 TIter next(fFunctions);
2907 TObject* obj = 0;
2908 while((obj = next())) {
2909 if(obj->InheritsFrom(TF1::Class())) {
2910 fPaintGraph->PaintStats((TF1*)obj);
2911 ((TF1*)obj)->Paint("sameC");
2912 }
2913 }
2914 }
2915
2916 return;
2917 }
2918
2919 //use TH2 for painting
2920 if(GetDimension() == 2) {
2921 if(!fPaintHisto) {
2923 }
2924 else
2926
2927 //paint histogram
2928 fPaintHisto->Paint(opt);
2929 return;
2930 }
2931 Warning("Paint","Painting 3D efficiency is not implemented");
2932}
2933
2934////////////////////////////////////////////////////////////////////////////////
2935/// Have histograms fixed bins along each axis?
2936
2937void TEfficiency::SavePrimitive(std::ostream& out,Option_t* opt)
2938{
2939 Bool_t equi_bins = true;
2940
2941 //indentation
2942 TString indent = " ";
2943 //names for arrays containing the bin edges
2944 //static counter needed if more objects are saved
2945 static Int_t naxis = 0;
2946 TString sxaxis="xAxis",syaxis="yAxis",szaxis="zAxis";
2947
2948 //note the missing break statements!
2949 switch(GetDimension()) {
2950 case 3:
2951 equi_bins = equi_bins && !fTotalHistogram->GetZaxis()->GetXbins()->fArray
2953 case 2:
2954 equi_bins = equi_bins && !fTotalHistogram->GetYaxis()->GetXbins()->fArray
2956 case 1:
2957 equi_bins = equi_bins && !fTotalHistogram->GetXaxis()->GetXbins()->fArray
2959 }
2960
2961 //create arrays containing the variable binning
2962 if(!equi_bins) {
2963 Int_t i;
2964 ++naxis;
2965 sxaxis += naxis;
2966 syaxis += naxis;
2967 szaxis += naxis;
2968 //x axis
2969 out << indent << "Double_t " << sxaxis << "["
2970 << fTotalHistogram->GetXaxis()->GetXbins()->fN << "] = {";
2971 for (i = 0; i < fTotalHistogram->GetXaxis()->GetXbins()->fN; ++i) {
2972 if (i != 0) out << ", ";
2973 out << fTotalHistogram->GetXaxis()->GetXbins()->fArray[i];
2974 }
2975 out << "}; " << std::endl;
2976 //y axis
2977 if(GetDimension() > 1) {
2978 out << indent << "Double_t " << syaxis << "["
2979 << fTotalHistogram->GetYaxis()->GetXbins()->fN << "] = {";
2980 for (i = 0; i < fTotalHistogram->GetYaxis()->GetXbins()->fN; ++i) {
2981 if (i != 0) out << ", ";
2982 out << fTotalHistogram->GetYaxis()->GetXbins()->fArray[i];
2983 }
2984 out << "}; " << std::endl;
2985 }
2986 //z axis
2987 if(GetDimension() > 2) {
2988 out << indent << "Double_t " << szaxis << "["
2989 << fTotalHistogram->GetZaxis()->GetXbins()->fN << "] = {";
2990 for (i = 0; i < fTotalHistogram->GetZaxis()->GetXbins()->fN; ++i) {
2991 if (i != 0) out << ", ";
2992 out << fTotalHistogram->GetZaxis()->GetXbins()->fArray[i];
2993 }
2994 out << "}; " << std::endl;
2995 }
2996 }//creating variable binning
2997
2998 //TEfficiency pointer has efficiency name + counter
2999 static Int_t eff_count = 0;
3000 ++eff_count;
3001 TString eff_name = GetName();
3002 eff_name += eff_count;
3003
3004 const char* name = eff_name.Data();
3005
3006 //construct TEfficiency object
3007 const char quote = '"';
3008 out << indent << std::endl;
3009 out << indent << ClassName() << " * " << name << " = new " << ClassName()
3010 << "(" << quote << GetName() << quote << "," << quote
3011 << GetTitle() << quote <<",";
3012 //fixed bin size -> use n,min,max constructor
3013 if(equi_bins) {
3014 out << fTotalHistogram->GetXaxis()->GetNbins() << ","
3015 << fTotalHistogram->GetXaxis()->GetXmin() << ","
3017 if(GetDimension() > 1) {
3018 out << "," << fTotalHistogram->GetYaxis()->GetNbins() << ","
3019 << fTotalHistogram->GetYaxis()->GetXmin() << ","
3021 }
3022 if(GetDimension() > 2) {
3023 out << "," << fTotalHistogram->GetZaxis()->GetNbins() << ","
3024 << fTotalHistogram->GetZaxis()->GetXmin() << ","
3026 }
3027 }
3028 //variable bin size -> use n,*bins constructor
3029 else {
3030 out << fTotalHistogram->GetXaxis()->GetNbins() << "," << sxaxis;
3031 if(GetDimension() > 1)
3032 out << "," << fTotalHistogram->GetYaxis()->GetNbins() << ","
3033 << syaxis;
3034 if(GetDimension() > 2)
3035 out << "," << fTotalHistogram->GetZaxis()->GetNbins() << ","
3036 << szaxis;
3037 }
3038 out << ");" << std::endl;
3039 out << indent << std::endl;
3040
3041 //set statistic options
3042 out << indent << name << "->SetConfidenceLevel(" << fConfLevel << ");"
3043 << std::endl;
3044 out << indent << name << "->SetBetaAlpha(" << fBeta_alpha << ");"
3045 << std::endl;
3046 out << indent << name << "->SetBetaBeta(" << fBeta_beta << ");" << std::endl;
3047 out << indent << name << "->SetWeight(" << fWeight << ");" << std::endl;
3048 out << indent << name << "->SetStatisticOption(" << fStatisticOption << ");"
3049 << std::endl;
3050 out << indent << name << "->SetPosteriorMode(" << TestBit(kPosteriorMode) << ");" << std::endl;
3051 out << indent << name << "->SetShortestInterval(" << TestBit(kShortestInterval) << ");" << std::endl;
3052 if(TestBit(kUseWeights))
3053 out << indent << name << "->SetUseWeightedEvents();" << std::endl;
3054
3055 // save bin-by-bin prior parameters
3056 for(unsigned int i = 0; i < fBeta_bin_params.size(); ++i)
3057 {
3058 out << indent << name << "->SetBetaBinParameters(" << i << "," << fBeta_bin_params.at(i).first
3059 << "," << fBeta_bin_params.at(i).second << ");" << std::endl;
3060 }
3061
3062 //set bin contents
3063 Int_t nbins = fTotalHistogram->GetNbinsX() + 2;
3064 if(GetDimension() > 1)
3065 nbins *= fTotalHistogram->GetNbinsY() + 2;
3066 if(GetDimension() > 2)
3067 nbins *= fTotalHistogram->GetNbinsZ() + 2;
3068
3069 //important: set first total number than passed number
3070 for(Int_t i = 0; i < nbins; ++i) {
3071 out << indent << name <<"->SetTotalEvents(" << i << "," <<
3072 fTotalHistogram->GetBinContent(i) << ");" << std::endl;
3073 out << indent << name <<"->SetPassedEvents(" << i << "," <<
3074 fPassedHistogram->GetBinContent(i) << ");" << std::endl;
3075 }
3076
3077 //save list of functions
3078 TIter next(fFunctions);
3079 TObject* obj = 0;
3080 while((obj = next())) {
3081 obj->SavePrimitive(out,"nodraw");
3082 if(obj->InheritsFrom(TF1::Class())) {
3083 out << indent << name << "->GetListOfFunctions()->Add("
3084 << obj->GetName() << ");" << std::endl;
3085 }
3086 }
3087
3088 //set style
3092
3093 //draw TEfficiency object
3094 TString option = opt;
3095 option.ToLower();
3096 if (!option.Contains("nodraw"))
3097 out<< indent << name<< "->Draw(" << quote << opt << quote << ");"
3098 << std::endl;
3099}
3100
3101////////////////////////////////////////////////////////////////////////////////
3102/// Sets the shape parameter &alpha;
3103///
3104/// The prior probability of the efficiency is given by the beta distribution:
3105/// \f[
3106/// f(\varepsilon;\alpha;\beta) = \frac{1}{B(\alpha,\beta)} \varepsilon^{\alpha-1} (1 - \varepsilon)^{\beta-1}
3107/// \f]
3108///
3109/// Note: - both shape parameters have to be positive (i.e. > 0)
3110
3112{
3113 if(alpha > 0)
3114 fBeta_alpha = alpha;
3115 else
3116 Warning("SetBetaAlpha(Double_t)","invalid shape parameter %.2lf",alpha);
3117}
3118
3119////////////////////////////////////////////////////////////////////////////////
3120/// Sets the shape parameter &beta;
3121///
3122/// The prior probability of the efficiency is given by the beta distribution:
3123/// \f[
3124/// f(\varepsilon;\alpha,\beta) = \frac{1}{B(\alpha,\beta)} \varepsilon^{\alpha-1} (1 - \varepsilon)^{\beta-1}
3125/// \f]
3126///
3127/// Note: - both shape parameters have to be positive (i.e. > 0)
3128
3130{
3131 if(beta > 0)
3132 fBeta_beta = beta;
3133 else
3134 Warning("SetBetaBeta(Double_t)","invalid shape parameter %.2lf",beta);
3135}
3136
3137////////////////////////////////////////////////////////////////////////////////
3138/// Sets different shape parameter &alpha; and &beta;
3139/// for the prior distribution for each bin. By default the global parameter are used if they are not set
3140/// for the specific bin
3141/// The prior probability of the efficiency is given by the beta distribution:
3142/// \f[
3143/// f(\varepsilon;\alpha;\beta) = \frac{1}{B(\alpha,\beta)} \varepsilon^{\alpha-1} (1 - \varepsilon)^{\beta-1}
3144/// \f]
3145///
3146/// Note:
3147/// - both shape parameters have to be positive (i.e. > 0)
3148/// - bin gives the global bin number (cf. GetGlobalBin)
3149
3151{
3152 if (!fPassedHistogram || !fTotalHistogram) return;
3154 // doing this I get h1->fN which is available only for a TH1D
3155 UInt_t n = h1->GetBin(h1->GetNbinsX()+1, h1->GetNbinsY()+1, h1->GetNbinsZ()+1 ) + 1;
3156
3157 // in case vector is not created do with default alpha, beta params
3158 if (fBeta_bin_params.size() != n )
3159 fBeta_bin_params = std::vector<std::pair<Double_t, Double_t> >(n, std::make_pair(fBeta_alpha, fBeta_beta) );
3160
3161 // vector contains also values for under/overflows
3162 fBeta_bin_params[bin] = std::make_pair(alpha,beta);
3163 SetBit(kUseBinPrior,true);
3164
3165}
3166
3167////////////////////////////////////////////////////////////////////////////////
3168/// Set the bins for the underlined passed and total histograms
3169/// If the class have been already filled the previous contents will be lost
3170
3172{
3173 if (GetDimension() != 1) {
3174 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3175 return kFALSE;
3176 }
3177 if (fTotalHistogram->GetEntries() != 0 ) {
3178 Warning("SetBins","Histogram entries will be lost after SetBins");
3181 }
3184 return kTRUE;
3185}
3186
3187////////////////////////////////////////////////////////////////////////////////
3188/// Set the bins for the underlined passed and total histograms
3189/// If the class have been already filled the previous contents will be lost
3190
3192{
3193 if (GetDimension() != 1) {
3194 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3195 return kFALSE;
3196 }
3197 if (fTotalHistogram->GetEntries() != 0 ) {
3198 Warning("SetBins","Histogram entries will be lost after SetBins");
3201 }
3202 fPassedHistogram->SetBins(nx,xBins);
3203 fTotalHistogram->SetBins(nx,xBins);
3204 return kTRUE;
3205}
3206
3207////////////////////////////////////////////////////////////////////////////////
3208/// Set the bins for the underlined passed and total histograms
3209/// If the class have been already filled the previous contents will be lost
3210
3212{
3213 if (GetDimension() != 2) {
3214 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3215 return kFALSE;
3216 }
3217 if (fTotalHistogram->GetEntries() != 0 ) {
3218 Warning("SetBins","Histogram entries will be lost after SetBins");
3221 }
3224 return kTRUE;
3225}
3226
3227////////////////////////////////////////////////////////////////////////////////
3228/// Set the bins for the underlined passed and total histograms
3229/// If the class have been already filled the previous contents will be lost
3230
3231Bool_t TEfficiency::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins)
3232{
3233 if (GetDimension() != 2) {
3234 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3235 return kFALSE;
3236 }
3237 if (fTotalHistogram->GetEntries() != 0 ) {
3238 Warning("SetBins","Histogram entries will be lost after SetBins");
3241 }
3242 fPassedHistogram->SetBins(nx,xBins,ny,yBins);
3243 fTotalHistogram->SetBins(nx,xBins,ny,yBins);
3244 return kTRUE;
3245}
3246
3247////////////////////////////////////////////////////////////////////////////////
3248/// Set the bins for the underlined passed and total histograms
3249/// If the class have been already filled the previous contents will be lost
3250
3252 Int_t nz, Double_t zmin, Double_t zmax)
3253{
3254 if (GetDimension() != 3) {
3255 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3256 return kFALSE;
3257 }
3258 if (fTotalHistogram->GetEntries() != 0 ) {
3259 Warning("SetBins","Histogram entries will be lost after SetBins");
3262 }
3263 fPassedHistogram->SetBins(nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
3264 fTotalHistogram->SetBins (nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
3265 return kTRUE;
3266}
3267
3268////////////////////////////////////////////////////////////////////////////////
3269/// Set the bins for the underlined passed and total histograms
3270/// If the class have been already filled the previous contents will be lost
3271
3272Bool_t TEfficiency::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins, Int_t nz,
3273 const Double_t *zBins )
3274{
3275 if (GetDimension() != 3) {
3276 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3277 return kFALSE;
3278 }
3279 if (fTotalHistogram->GetEntries() != 0 ) {
3280 Warning("SetBins","Histogram entries will be lost after SetBins");
3283 }
3284 fPassedHistogram->SetBins(nx,xBins,ny,yBins,nz,zBins);
3285 fTotalHistogram->SetBins(nx,xBins,ny,yBins,nz,zBins);
3286 return kTRUE;
3287}
3288
3289////////////////////////////////////////////////////////////////////////////////
3290/// Sets the confidence level (0 < level < 1)
3291/// The default value is 1-sigma :~ 0.683
3292
3294{
3295 if((level > 0) && (level < 1))
3296 fConfLevel = level;
3297 else
3298 Warning("SetConfidenceLevel(Double_t)","invalid confidence level %.2lf",level);
3299}
3300
3301////////////////////////////////////////////////////////////////////////////////
3302/// Sets the directory holding this TEfficiency object
3303///
3304/// A reference to this TEfficiency object is removed from the current
3305/// directory (if it exists) and a new reference to this TEfficiency object is
3306/// added to the given directory.
3307///
3308/// Notes: - If the given directory is 0, the TEfficiency object does not
3309/// belong to any directory and will not be written to file during the
3310/// next TFile::Write() command.
3311
3313{
3314 if(fDirectory == dir)
3315 return;
3316 if(fDirectory)
3317 fDirectory->Remove(this);
3318 fDirectory = dir;
3319 if(fDirectory)
3320 fDirectory->Append(this);
3321}
3322
3323////////////////////////////////////////////////////////////////////////////////
3324/// Sets the name
3325///
3326/// Note: The names of the internal histograms are set to "name + _total" and
3327/// "name + _passed" respectively.
3328
3330{
3332
3333 //setting the names (appending the correct ending)
3334 TString name_total = name + TString("_total");
3335 TString name_passed = name + TString("_passed");
3336 fTotalHistogram->SetName(name_total);
3337 fPassedHistogram->SetName(name_passed);
3338}
3339
3340////////////////////////////////////////////////////////////////////////////////
3341/// Sets the number of passed events in the given global bin
3342///
3343/// returns "true" if the number of passed events has been updated
3344/// otherwise "false" ist returned
3345///
3346/// Note: - requires: 0 <= events <= fTotalHistogram->GetBinContent(bin)
3347
3349{
3350 if(events <= fTotalHistogram->GetBinContent(bin)) {
3351 fPassedHistogram->SetBinContent(bin,events);
3352 return true;
3353 }
3354 else {
3355 Error("SetPassedEvents(Int_t,Int_t)","total number of events (%.1lf) in bin %i is less than given number of passed events %i",fTotalHistogram->GetBinContent(bin),bin,events);
3356 return false;
3357 }
3358}
3359
3360////////////////////////////////////////////////////////////////////////////////
3361/// Sets the histogram containing the passed events
3362///
3363/// The given histogram is cloned and stored internally as histogram containing
3364/// the passed events. The given histogram has to be consistent with the current
3365/// fTotalHistogram (see CheckConsistency(const TH1&,const TH1&)).
3366/// The method returns whether the fPassedHistogram has been replaced (true) or
3367/// not (false).
3368///
3369/// Note: The list of associated functions fFunctions is cleared.
3370///
3371/// Option:
3372/// - "f": force the replacement without checking the consistency
3373/// This can lead to inconsistent histograms and useless results
3374/// or unexpected behaviour. But sometimes it might be the only
3375/// way to change the histograms. If you use this option, you
3376/// should ensure that the fTotalHistogram is replaced by a
3377/// consistent one (with respect to rPassed) as well.
3378
3380{
3381 TString option = opt;
3382 option.ToLower();
3383
3384 Bool_t bReplace = option.Contains("f");
3385
3386 if(!bReplace)
3387 bReplace = CheckConsistency(rPassed,*fTotalHistogram);
3388
3389 if(bReplace) {
3390 delete fPassedHistogram;
3391 // do not add cloned histogram to gDirectory
3392 {
3393 TDirectory::TContext ctx(nullptr);
3394 fPassedHistogram = (TH1*)(rPassed.Clone());
3396 }
3397
3398 if(fFunctions)
3399 fFunctions->Delete();
3400
3401 //check whether both histograms are filled with weights
3402 bool useWeights = CheckWeights(rPassed,*fTotalHistogram);
3403
3404 SetUseWeightedEvents(useWeights);
3405
3406 return true;
3407 }
3408 else
3409 return false;
3410}
3411
3412////////////////////////////////////////////////////////////////////////////////
3413/// Sets the statistic option which affects the calculation of the confidence interval
3414///
3415/// Options:
3416/// - kFCP (=0)(default): using the Clopper-Pearson interval (recommended by PDG)
3417/// sets kIsBayesian = false
3418/// see also ClopperPearson
3419/// - kFNormal (=1) : using the normal approximation
3420/// sets kIsBayesian = false
3421/// see also Normal
3422/// - kFWilson (=2) : using the Wilson interval
3423/// sets kIsBayesian = false
3424/// see also Wilson
3425/// - kFAC (=3) : using the Agresti-Coull interval
3426/// sets kIsBayesian = false
3427/// see also AgrestiCoull
3428/// - kFFC (=4) : using the Feldman-Cousins frequentist method
3429/// sets kIsBayesian = false
3430/// see also FeldmanCousins
3431/// - kBJeffrey (=5) : using the Jeffrey interval
3432/// sets kIsBayesian = true, fBeta_alpha = 0.5 and fBeta_beta = 0.5
3433/// see also Bayesian
3434/// - kBUniform (=6) : using a uniform prior
3435/// sets kIsBayesian = true, fBeta_alpha = 1 and fBeta_beta = 1
3436/// see also Bayesian
3437/// - kBBayesian (=7) : using a custom prior defined by fBeta_alpha and fBeta_beta
3438/// sets kIsBayesian = true
3439/// see also Bayesian
3440/// - kMidP (=8) : using the Lancaster Mid-P method
3441/// sets kIsBayesian = false
3442
3443
3445{
3446 fStatisticOption = option;
3447
3448 switch(option)
3449 {
3450 case kFCP:
3452 SetBit(kIsBayesian,false);
3453 break;
3454 case kFNormal:
3455 fBoundary = &Normal;
3456 SetBit(kIsBayesian,false);
3457 break;
3458 case kFWilson:
3459 fBoundary = &Wilson;
3460 SetBit(kIsBayesian,false);
3461 break;
3462 case kFAC:
3464 SetBit(kIsBayesian,false);
3465 break;
3466 case kFFC:
3468 SetBit(kIsBayesian,false);
3469 break;
3470 case kMidP:
3472 SetBit(kIsBayesian,false);
3473 break;
3474 case kBJeffrey:
3475 fBeta_alpha = 0.5;
3476 fBeta_beta = 0.5;
3477 SetBit(kIsBayesian,true);
3478 SetBit(kUseBinPrior,false);
3479 break;
3480 case kBUniform:
3481 fBeta_alpha = 1;
3482 fBeta_beta = 1;
3483 SetBit(kIsBayesian,true);
3484 SetBit(kUseBinPrior,false);
3485 break;
3486 case kBBayesian:
3487 SetBit(kIsBayesian,true);
3488 break;
3489 default:
3492 SetBit(kIsBayesian,false);
3493 }
3494}
3495
3496////////////////////////////////////////////////////////////////////////////////
3497/// Sets the title
3498///
3499/// Notes:
3500/// - The titles of the internal histograms are set to "title + (total)"
3501/// or "title + (passed)" respectively.
3502/// - It is possible to label the axis of the histograms as usual (see
3503/// TH1::SetTitle).
3504///
3505/// Example: Setting the title to "My Efficiency" and label the axis
3506/// pEff->SetTitle("My Efficiency;x label;eff");
3507
3508void TEfficiency::SetTitle(const char* title)
3509{
3510
3511 //setting the titles (looking for the first semicolon and insert the tokens there)
3512 TString title_passed = title;
3513 TString title_total = title;
3514 Ssiz_t pos = title_passed.First(";");
3515 if (pos != kNPOS) {
3516 title_passed.Insert(pos," (passed)");
3517 title_total.Insert(pos," (total)");
3518 }
3519 else {
3520 title_passed.Append(" (passed)");
3521 title_total.Append(" (total)");
3522 }
3523 fPassedHistogram->SetTitle(title_passed);
3524 fTotalHistogram->SetTitle(title_total);
3525
3526 // strip (total) for the TEfficiency title
3527 // HIstogram SetTitle has already stripped the axis
3528 TString teffTitle = fTotalHistogram->GetTitle();
3529 teffTitle.ReplaceAll(" (total)","");
3530 TNamed::SetTitle(teffTitle);
3531
3532}
3533
3534////////////////////////////////////////////////////////////////////////////////
3535/// Sets the number of total events in the given global bin
3536///
3537/// returns "true" if the number of total events has been updated
3538/// otherwise "false" ist returned
3539///
3540/// Note: - requires: fPassedHistogram->GetBinContent(bin) <= events
3541
3543{
3544 if(events >= fPassedHistogram->GetBinContent(bin)) {
3545 fTotalHistogram->SetBinContent(bin,events);
3546 return true;
3547 }
3548 else {
3549 Error("SetTotalEvents(Int_t,Int_t)","passed number of events (%.1lf) in bin %i is bigger than given number of total events %i",fPassedHistogram->GetBinContent(bin),bin,events);
3550 return false;
3551 }
3552}
3553
3554////////////////////////////////////////////////////////////////////////////////
3555/// Sets the histogram containing all events
3556///
3557/// The given histogram is cloned and stored internally as histogram containing
3558/// all events. The given histogram has to be consistent with the current
3559/// fPassedHistogram (see CheckConsistency(const TH1&,const TH1&)).
3560/// The method returns whether the fTotalHistogram has been replaced (true) or
3561/// not (false).
3562///
3563/// Note: The list of associated functions fFunctions is cleared.
3564///
3565/// Option:
3566/// - "f": force the replacement without checking the consistency
3567/// This can lead to inconsistent histograms and useless results
3568/// or unexpected behaviour. But sometimes it might be the only
3569/// way to change the histograms. If you use this option, you
3570/// should ensure that the fPassedHistogram is replaced by a
3571/// consistent one (with respect to rTotal) as well.
3572
3574{
3575 TString option = opt;
3576 option.ToLower();
3577
3578 Bool_t bReplace = option.Contains("f");
3579
3580 if(!bReplace)
3581 bReplace = CheckConsistency(*fPassedHistogram,rTotal);
3582
3583 if(bReplace) {
3584 delete fTotalHistogram;
3585 // do not add cloned histogram to gDirectory
3586 {
3587 TDirectory::TContext ctx(nullptr);
3588 fTotalHistogram = (TH1*)(rTotal.Clone());
3589 }
3591
3592 if(fFunctions)
3593 fFunctions->Delete();
3594
3595 //check whether both histograms are filled with weights
3596 bool useWeights = CheckWeights(*fPassedHistogram,rTotal);
3597 SetUseWeightedEvents(useWeights);
3598
3599 return true;
3600 }
3601 else
3602 return false;
3603}
3604
3605////////////////////////////////////////////////////////////////////////////////
3606
3608{
3609 if (on && !TestBit(kUseWeights) )
3610 gROOT->Info("TEfficiency::SetUseWeightedEvents","Handle weighted events for computing efficiency");
3611
3612 SetBit(kUseWeights,on);
3613
3618}
3619
3620////////////////////////////////////////////////////////////////////////////////
3621/// Sets the global weight for this TEfficiency object
3622///
3623/// Note: - weight has to be positive ( > 0)
3624
3626{
3627 if(weight > 0)
3628 fWeight = weight;
3629 else
3630 Warning("SetWeight","invalid weight %.2lf",weight);
3631}
3632
3633////////////////////////////////////////////////////////////////////////////////
3634/**
3635Calculates the boundaries for the frequentist Wilson interval
3636
3637\param[in] total number of total events
3638\param[in] passed 0 <= number of passed events <= total
3639\param[in] level confidence level
3640\param[in] bUpper
3641 - true - upper boundary is returned
3642 - false - lower boundary is returned
3643
3644Calculation:
3645\f{eqnarray*}{
3646 \alpha &=& 1 - \frac{level}{2}\\
3647 \kappa &=& \Phi^{-1}(1 - \alpha,1) ...\ normal\ quantile\ function\\
3648 mode &=& \frac{passed + \frac{\kappa^{2}}{2}}{total + \kappa^{2}}\\
3649 \Delta &=& \frac{\kappa}{total + \kappa^{2}} * \sqrt{passed (1 - \frac{passed}{total}) + \frac{\kappa^{2}}{4}}\\
3650 return &=& max(0,mode - \Delta)\ or\ min(1,mode + \Delta)
3651\f}
3652
3653*/
3654
3656{
3657 Double_t alpha = (1.0 - level)/2;
3658 if (total == 0) return (bUpper) ? 1 : 0;
3659 Double_t average = ((Double_t)passed) / total;
3660 Double_t kappa = ROOT::Math::normal_quantile(1 - alpha,1);
3661
3662 Double_t mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa);
3663 Double_t delta = kappa / (total + kappa*kappa) * std::sqrt(total * average
3664 * (1 - average) + kappa * kappa / 4);
3665 if(bUpper)
3666 return ((mode + delta) > 1) ? 1.0 : (mode + delta);
3667 else
3668 return ((mode - delta) < 0) ? 0.0 : (mode - delta);
3669}
3670
3671////////////////////////////////////////////////////////////////////////////////
3672/// Addition operator
3673///
3674/// adds the corresponding histograms:
3675/// ~~~ {.cpp}
3676/// lhs.GetTotalHistogram() + rhs.GetTotalHistogram()
3677/// lhs.GetPassedHistogram() + rhs.GetPassedHistogram()
3678/// ~~~
3679/// the statistic option and the confidence level are taken from lhs
3680
3681const TEfficiency operator+(const TEfficiency& lhs,const TEfficiency& rhs)
3682{
3683 TEfficiency tmp(lhs);
3684 tmp += rhs;
3685 return tmp;
3686}
3687
3688#endif
void Class()
Definition: Class.C:29
double
Definition: Converters.cxx:939
#define b(i)
Definition: RSha256.hxx:100
#define e(i)
Definition: RSha256.hxx:103
const Ssiz_t kNPOS
Definition: RtypesCore.h:124
int Int_t
Definition: RtypesCore.h:45
int Ssiz_t
Definition: RtypesCore.h:67
unsigned int UInt_t
Definition: RtypesCore.h:46
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
long long Long64_t
Definition: RtypesCore.h:80
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:364
static void indent(ostringstream &buf, int indent_level)
#define gDirectory
Definition: TDirectory.h:348
const TEfficiency operator+(const TEfficiency &lhs, const TEfficiency &rhs)
Addition operator.
const Double_t kDefBetaAlpha
Definition: TEfficiency.cxx:36
const Double_t kDefWeight
Definition: TEfficiency.cxx:40
const Double_t kDefBetaBeta
Definition: TEfficiency.cxx:37
const TEfficiency::EStatOption kDefStatOpt
Definition: TEfficiency.cxx:39
const Double_t kDefConfLevel
Definition: TEfficiency.cxx:38
static unsigned int total
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
#define gROOT
Definition: TROOT.h:404
#define gPad
Definition: TVirtualPad.h:288
static struct mg_connection * fc(struct mg_context *ctx)
Definition: civetweb.c:3728
User class for performing function minimization.
virtual bool Minimize(int maxIter, double absTol=1.E-8, double relTol=1.E-10)
Find minimum position iterating until convergence specified by the absolute and relative tolerance or...
virtual double XMinimum() const
Return current estimate of the position of the minimum.
void SetFunction(const ROOT::Math::IGenFunction &f, double xlow, double xup)
Sets function to be minimized.
void SetNpx(int npx)
Set the number of point used to bracket root using a grid.
virtual double FValMinimum() const
Return function value at current estimate of the minimum.
Template class to wrap any C++ callable object which takes one argument i.e.
Double_t At(Int_t i) const
Definition: TArrayD.h:79
Double_t * fArray
Definition: TArrayD.h:30
const Double_t * GetArray() const
Definition: TArrayD.h:43
Int_t fN
Definition: TArray.h:38
Fill Area Attributes class.
Definition: TAttFill.h:19
void Copy(TAttFill &attfill) const
Copy this fill attributes to a new TAttFill.
Definition: TAttFill.cxx:204
virtual void SaveFillAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1001)
Save fill attributes as C++ statement(s) on output stream out.
Definition: TAttFill.cxx:236
Line Attributes class.
Definition: TAttLine.h:18
void Copy(TAttLine &attline) const
Copy this line attributes to a new TAttLine.
Definition: TAttLine.cxx:175
virtual void SaveLineAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1, Int_t widdef=1)
Save line attributes as C++ statement(s) on output stream out.
Definition: TAttLine.cxx:273
Marker Attributes class.
Definition: TAttMarker.h:19
virtual void SaveMarkerAttributes(std::ostream &out, const char *name, Int_t coldef=1, Int_t stydef=1, Int_t sizdef=1)
Save line attributes as C++ statement(s) on output stream out.
Definition: TAttMarker.cxx:345
void Copy(TAttMarker &attmarker) const
Copy this marker attributes to a new TAttMarker.
Definition: TAttMarker.cxx:241
Class to manage histogram axis.
Definition: TAxis.h:30
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
Definition: TAxis.cxx:823
Bool_t IsVariableBinSize() const
Definition: TAxis.h:136
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:478
const TArrayD * GetXbins() const
Definition: TAxis.h:130
Double_t GetXmax() const
Definition: TAxis.h:134
const char * GetBinLabel(Int_t bin) const
Return label for bin.
Definition: TAxis.cxx:440
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:518
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:419
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
const char * GetTitle() const
Returns title of object.
Definition: TAxis.h:129
THashList * GetLabels() const
Definition: TAxis.h:117
Binomial fitter for the division of two histograms.
TFitResultPtr Fit(TF1 *f1, Option_t *option="")
Carry out the fit of the given function to the given histograms.
Collection abstract base class.
Definition: TCollection.h:65
virtual Bool_t IsEmpty() const
Definition: TCollection.h:188
TDirectory::TContext keeps track and restore the current directory.
Definition: TDirectory.h:89
Describe directory structure in memory.
Definition: TDirectory.h:45
virtual void Append(TObject *obj, Bool_t replace=kFALSE)
Append object to this directory.
Definition: TDirectory.cxx:191
virtual TObject * Remove(TObject *)
Remove an object from the in-memory list.
Class to handle efficiency histograms.
Definition: TEfficiency.h:28
static Bool_t FeldmanCousinsInterval(Double_t total, Double_t passed, Double_t level, Double_t &lower, Double_t &upper)
Calculates the interval boundaries using the frequentist methods of Feldman-Cousins.
static Double_t BetaMode(Double_t alpha, Double_t beta)
Compute the mode of the beta distribution.
Bool_t SetPassedEvents(Int_t bin, Int_t events)
Sets the number of passed events in the given global bin.
TH2 * CreateHistogram(Option_t *opt="") const
Create the histogram used to be painted (for dim=2 TEfficiency) The return object is managed by the c...
static Bool_t BetaShortestInterval(Double_t level, Double_t alpha, Double_t beta, Double_t &lower, Double_t &upper)
Calculates the boundaries for a shortest confidence interval for a Beta distribution.
static Bool_t CheckWeights(const TH1 &pass, const TH1 &total)
Check if both histogram are weighted.
static Double_t BetaMean(Double_t alpha, Double_t beta)
Compute the mean (average) of the beta distribution.
TEfficiency()
Default constructor.
Double_t GetBetaAlpha(Int_t bin=-1) const
Definition: TEfficiency.h:106
void FillWeighted(Bool_t bPassed, Double_t weight, Double_t x, Double_t y=0, Double_t z=0)
This function is used for filling the two histograms with a weight.
~TEfficiency()
default destructor
TList * GetListOfFunctions()
static Double_t Bayesian(Double_t total, Double_t passed, Double_t level, Double_t alpha, Double_t beta, Bool_t bUpper, Bool_t bShortest=false)
Calculates the boundaries for a Bayesian confidence interval (shortest or central interval depending ...
static Double_t AgrestiCoull(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Agresti-Coull interval.
Long64_t Merge(TCollection *list)
Merges the TEfficiency objects in the given list to the given TEfficiency object using the operator+=...
std::vector< std::pair< Double_t, Double_t > > fBeta_bin_params
Parameter for prior beta distribution different bin by bin (default vector is empty)
Definition: TEfficiency.h:48
static Double_t FeldmanCousins(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Feldman-Cousins interval.
EStatOption fStatisticOption
Defines how the confidence intervals are determined.
Definition: TEfficiency.h:57
void SetStatisticOption(EStatOption option)
Sets the statistic option which affects the calculation of the confidence interval.
void SetWeight(Double_t weight)
Sets the global weight for this TEfficiency object.
TH1 * fTotalHistogram
Histogram for total number of events.
Definition: TEfficiency.h:58
Int_t GetDimension() const
returns the dimension of the current TEfficiency object
TEfficiency & operator+=(const TEfficiency &rhs)
Adds the histograms of another TEfficiency object to current histograms.
Bool_t SetBins(Int_t nx, Double_t xmin, Double_t xmax)
Set the bins for the underlined passed and total histograms If the class have been already filled the...
void Build(const char *name, const char *title)
Building standard data structure of a TEfficiency object.
TH1 * GetCopyPassedHisto() const
Returns a cloned version of fPassedHistogram.
Double_t GetEfficiencyErrorUp(Int_t bin) const
Returns the upper error on the efficiency in the given global bin.
void Draw(Option_t *opt="")
Draws the current TEfficiency object.
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a graph.
Double_t fBeta_alpha
Global parameter for prior beta distribution (default = 1)
Definition: TEfficiency.h:46
Bool_t UsesBayesianStat() const
Definition: TEfficiency.h:156
void SetBetaBeta(Double_t beta)
Sets the shape parameter β.
Double_t GetConfidenceLevel() const
Definition: TEfficiency.h:108
static Bool_t CheckBinning(const TH1 &pass, const TH1 &total)
Checks binning for each axis.
static Double_t BetaCentralInterval(Double_t level, Double_t alpha, Double_t beta, Bool_t bUpper)
Calculates the boundaries for a central confidence interval for a Beta distribution.
Int_t GetGlobalBin(Int_t binx, Int_t biny=0, Int_t binz=0) const
Returns the global bin number which can be used as argument for the following functions:
TH1 * fPassedHistogram
Histogram for events which passed certain criteria.
Definition: TEfficiency.h:56
static Double_t MidPInterval(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries using the mid-P binomial interval (Lancaster method) from B.
void SetBetaAlpha(Double_t alpha)
Sets the shape parameter α.
@ kIsBayesian
Bayesian statistics are used.
Definition: TEfficiency.h:62
@ kUseWeights
Use weights.
Definition: TEfficiency.h:66
@ kPosteriorMode
Use posterior mean for best estimate (Bayesian statistics)
Definition: TEfficiency.h:63
@ kUseBinPrior
Use a different prior for each bin.
Definition: TEfficiency.h:65
@ kShortestInterval
Use shortest interval.
Definition: TEfficiency.h:64
static Bool_t CheckEntries(const TH1 &pass, const TH1 &total, Option_t *opt="")
Checks whether bin contents are compatible with binomial statistics.
static Double_t Normal(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Returns the confidence limits for the efficiency supposing that the efficiency follows a normal distr...
Double_t fWeight
Weight for all events (default = 1)
Definition: TEfficiency.h:59
Bool_t SetPassedHistogram(const TH1 &rPassed, Option_t *opt)
Sets the histogram containing the passed events.
Bool_t SetTotalEvents(Int_t bin, Int_t events)
Sets the number of total events in the given global bin.
Double_t GetBetaBeta(Int_t bin=-1) const
Definition: TEfficiency.h:107
Double_t(* fBoundary)(Double_t, Double_t, Double_t, Bool_t)
! Pointer to a method calculating the boundaries of confidence intervals
Definition: TEfficiency.h:50
void FillGraph(TGraphAsymmErrors *graph, Option_t *opt) const
Fill the graph to be painted with information from TEfficiency Internal method called by TEfficiency:...
void SetName(const char *name)
Sets the name.
void FillHistogram(TH2 *h2) const
Fill the 2d histogram to be painted with information from TEfficiency 2D Internal method called by TE...
Int_t FindFixBin(Double_t x, Double_t y=0, Double_t z=0) const
Returns the global bin number containing the given values.
TDirectory * fDirectory
! Pointer to directory holding this TEfficiency object
Definition: TEfficiency.h:52
static Double_t Combine(Double_t &up, Double_t &low, Int_t n, const Int_t *pass, const Int_t *total, Double_t alpha, Double_t beta, Double_t level=0.683, const Double_t *w=0, Option_t *opt="")
void SetUseWeightedEvents(Bool_t on=kTRUE)
static Double_t Wilson(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Wilson interval.
TEfficiency & operator=(const TEfficiency &rhs)
Assignment operator.
Double_t fConfLevel
Confidence level (default = 0.683, 1 sigma)
Definition: TEfficiency.h:51
Double_t fBeta_beta
Global parameter for prior beta distribution (default = 1)
Definition: TEfficiency.h:47
void SavePrimitive(std::ostream &out, Option_t *opt="")
Have histograms fixed bins along each axis?
Double_t GetEfficiency(Int_t bin) const
Returns the efficiency in the given global bin.
Bool_t SetTotalHistogram(const TH1 &rTotal, Option_t *opt)
Sets the histogram containing all events.
void Fill(Bool_t bPassed, Double_t x, Double_t y=0, Double_t z=0)
This function is used for filling the two histograms.
void SetDirectory(TDirectory *dir)
Sets the directory holding this TEfficiency object.
TGraphAsymmErrors * fPaintGraph
! Temporary graph for painting
Definition: TEfficiency.h:54
TGraphAsymmErrors * CreateGraph(Option_t *opt="") const
Create the graph used be painted (for dim=1 TEfficiency) The return object is managed by the caller.
EStatOption GetStatisticOption() const
Definition: TEfficiency.h:121
TList * fFunctions
->Pointer to list of functions
Definition: TEfficiency.h:53
void Paint(Option_t *opt)
Paints this TEfficiency object.
void SetBetaBinParameters(Int_t bin, Double_t alpha, Double_t beta)
Sets different shape parameter α and β for the prior distribution for each bin.
static Bool_t CheckConsistency(const TH1 &pass, const TH1 &total, Option_t *opt="")
Checks the consistence of the given histograms.
Double_t GetWeight() const
Definition: TEfficiency.h:123
TH1 * GetCopyTotalHisto() const
Returns a cloned version of fTotalHistogram.
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
static Double_t ClopperPearson(Double_t total, Double_t passed, Double_t level, Bool_t bUpper)
Calculates the boundaries for the frequentist Clopper-Pearson interval.
void SetConfidenceLevel(Double_t level)
Sets the confidence level (0 < level < 1) The default value is 1-sigma :~ 0.683.
Double_t GetEfficiencyErrorLow(Int_t bin) const
Returns the lower error on the efficiency in the given global bin.
EStatOption
Enumeration type for different statistic options for calculating confidence intervals kF* ....
Definition: TEfficiency.h:32
@ kBJeffrey
Jeffrey interval (Prior ~ Beta(0.5,0.5)
Definition: TEfficiency.h:38
@ kFWilson
Wilson interval.
Definition: TEfficiency.h:35
@ kFAC
Agresti-Coull interval.
Definition: TEfficiency.h:36
@ kMidP
Mid-P Lancaster interval.
Definition: TEfficiency.h:41
@ kBUniform
Prior ~ Uniform = Beta(1,1)
Definition: TEfficiency.h:39
@ kFFC
Feldman-Cousins interval.
Definition: TEfficiency.h:37
@ kBBayesian
User specified Prior ~ Beta(fBeta_alpha,fBeta_beta)
Definition: TEfficiency.h:40
@ kFNormal
Normal approximation.
Definition: TEfficiency.h:34
@ kFCP
Clopper-Pearson interval (recommended by PDG)
Definition: TEfficiency.h:33
void SetTitle(const char *title)
Sets the title.
TFitResultPtr Fit(TF1 *f1, Option_t *opt="")
Fits the efficiency using the TBinomialEfficiencyFitter class.
TH2 * fPaintHisto
! Temporary histogram for painting
Definition: TEfficiency.h:55
1-Dim function class
Definition: TF1.h:213
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:32
TGraph with asymmetric error bars.
virtual void Paint(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:2038
virtual void PaintStats(TF1 *fit)
Draw the stats.
Definition: TGraph.cxx:2065
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a graph.
Definition: TGraph.cxx:813
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TGraph.cxx:990
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
virtual void SetDirectory(TDirectory *dir)
By default, when a histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8780
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6680
virtual void SetNormFactor(Double_t factor=1)
Definition: TH1.h:404
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8984
TAxis * GetZaxis()
Definition: TH1.h:322
virtual void GetStats(Double_t *stats) const
fill the array stats from the contents of this histogram The array stats must be correctly dimensione...
Definition: TH1.cxx:7721
virtual Int_t GetNbinsY() const
Definition: TH1.h:297
virtual Int_t GetNbinsZ() const
Definition: TH1.h:298
virtual Int_t GetDimension() const
Definition: TH1.h:282
@ kIsAverage
Bin contents are average (used by Add)
Definition: TH1.h:170
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:7071
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:320
virtual Int_t GetNcells() const
Definition: TH1.h:299
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2741
virtual Int_t GetBin(Int_t binx, Int_t biny=0, Int_t binz=0) const
Return Global bin number corresponding to binx,y,z.
Definition: TH1.cxx:4894
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition: TH1.cxx:823
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3351
TAxis * GetYaxis()
Definition: TH1.h:321
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9065
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8995
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4387
@ kNstat
Size of statistics data (up to TProfile3D)
Definition: TH1.h:183
virtual void SetName(const char *name)
Change the name of this histogram.
Definition: TH1.cxx:8803
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4994
virtual TArrayD * GetSumw2()
Definition: TH1.h:312
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TH1.cxx:3247
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:9006
virtual void Paint(Option_t *option="")
Control routine to paint any kind of histograms.
Definition: TH1.cxx:6157
virtual Int_t GetSumw2N() const
Definition: TH1.h:314
virtual void SetBins(Int_t nx, Double_t xmin, Double_t xmax)
Redefine x axis parameters.
Definition: TH1.cxx:8610
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8863
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a line.
Definition: TH1.cxx:2812
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8833
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:251
Service class for 2-D histogram classes.
Definition: TH2.h:30
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2571
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:305
The 3-D histogram classes derived from the 1-D histogram classes.
Definition: TH3.h:31
A doubly linked list.
Definition: TList.h:44
virtual void Add(TObject *obj)
Definition: TList.h:87
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:822
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:470
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:659
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:29
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:359
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
@ kNotDeleted
object has not been deleted
Definition: TObject.h:78
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
virtual void AppendPad(Option_t *option="")
Append graphics object to current pad.
Definition: TObject.cxx:107
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TObject.cxx:666
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:696
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:445
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:921
void ResetBit(UInt_t f)
Definition: TObject.h:186
@ kInvalidObject
if object ctor succeeded but object should not be used
Definition: TObject.h:68
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:867
Basic string class.
Definition: TString.h:136
Ssiz_t Length() const
Definition: TString.h:410
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1150
TString & Insert(Ssiz_t pos, const char *s)
Definition: TString.h:649
Ssiz_t First(char c) const
Find first occurrence of a character c.
Definition: TString.cxx:523
const char * Data() const
Definition: TString.h:369
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:692
Bool_t IsNull() const
Definition: TString.h:407
TString & Append(const char *cs)
Definition: TString.h:564
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:624
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:639
double beta_pdf(double x, double a, double b)
Probability density function of the beta distribution.
double beta_cdf(double x, double a, double b)
Cumulative distribution function of the beta distribution Upper tail of the integral of the beta_pdf.
double beta_cdf_c(double x, double a, double b)
Complement of the cumulative distribution function of the beta distribution.
double normal_quantile(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the lower tail of the normal (Gaussian) distri...
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
double beta_quantile_c(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the lower tail of the beta distribution (beta_...
double beta_quantile(double x, double a, double b)
Inverse ( ) of the cumulative distribution function of the upper tail of the beta distribution (beta_...
double beta(double x, double y)
Calculates the beta function.
const Double_t sigma
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TGraphErrors * gr
Definition: legend1.C:25
TH1F * h1
Definition: legend1.C:5
TF1 * f1
Definition: legend1.C:11
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:430
const double xbins[xbins_n]
Definition: graph.py:1
Beta_interval_length(Double_t level, Double_t alpha, Double_t beta)
Double_t operator()(double lower) const
auto * a
Definition: textangle.C:12