Logo ROOT   6.18/05
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////////////////////////////////////////////////////////////////////////////////
615///default constructor
616///
617///should not be used explicitly
618
620fBeta_alpha(kDefBetaAlpha),
621fBeta_beta(kDefBetaBeta),
622fBoundary(0),
623fConfLevel(kDefConfLevel),
624fDirectory(0),
625fFunctions(0),
626fPaintGraph(0),
627fPaintHisto(0),
628fPassedHistogram(0),
629fTotalHistogram(0),
630fWeight(kDefWeight)
631{
633
634 // create 2 dummy histograms
635 fPassedHistogram = new TH1F("h_passed","passed",10,0,10);
636 fTotalHistogram = new TH1F("h_total","total",10,0,10);
637}
638
639////////////////////////////////////////////////////////////////////////////////
640///constructor using two existing histograms as input
641///
642///Input: passed - contains the events fulfilling some criteria
643/// total - contains all investigated events
644///
645///Notes: - both histograms have to fulfill the conditions of CheckConsistency
646/// - dimension of the resulting efficiency object depends
647/// on the dimension of the given histograms
648/// - Clones of both histograms are stored internally
649/// - The function SetName(total.GetName() + "_clone") is called to set
650/// the names of the new object and the internal histograms..
651/// - The created TEfficiency object is NOT appended to a directory. It
652/// will not be written to disk during the next TFile::Write() command
653/// in order to prevent duplication of data. If you want to save this
654/// TEfficiency object anyway, you can either append it to a
655/// directory by calling SetDirectory(TDirectory*) or write it
656/// explicitly to disk by calling Write().
657
658TEfficiency::TEfficiency(const TH1& passed,const TH1& total):
659fBeta_alpha(kDefBetaAlpha),
660fBeta_beta(kDefBetaBeta),
661fConfLevel(kDefConfLevel),
662fDirectory(0),
663fFunctions(0),
664fPaintGraph(0),
665fPaintHisto(0),
666fWeight(kDefWeight)
667{
668 //check consistency of histograms
669 if(CheckConsistency(passed,total)) {
672 fTotalHistogram = (TH1*)total.Clone();
673 fPassedHistogram = (TH1*)passed.Clone();
674 TH1::AddDirectory(bStatus);
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
693 fTotalHistogram = new TH1D("h1_total","h1 (total)",10,0,10);
694 fPassedHistogram = new TH1D("h1_passed","h1 (passed)",10,0,10);
695 TH1::AddDirectory(bStatus);
696 }
697
698 SetBit(kPosteriorMode,false);
700
702 SetDirectory(0);
703}
704
705////////////////////////////////////////////////////////////////////////////////
706/// Create 1-dimensional TEfficiency object with variable bin size.
707///
708/// Constructor creates two new and empty histograms with a given binning
709///
710/// Input:
711///
712/// - `name`: the common part of the name for both histograms (no blanks)
713/// fTotalHistogram has name: name + "_total"
714/// fPassedHistogram has name: name + "_passed"
715/// - `title`: the common part of the title for both histogram
716/// fTotalHistogram has title: title + " (total)"
717/// fPassedHistogram has title: title + " (passed)"
718/// It is possible to label the axis by passing a title with
719/// the following format: "title;xlabel;ylabel".
720/// - `nbins`: number of bins on the x-axis
721/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
722/// xbins[nbinsx] ... lower edge for overflow bin
723
724TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbins,
725 const Double_t* xbins):
726fBeta_alpha(kDefBetaAlpha),
727fBeta_beta(kDefBetaBeta),
728fConfLevel(kDefConfLevel),
729fDirectory(0),
730fFunctions(0),
731fPaintGraph(0),
732fPaintHisto(0),
733fWeight(kDefWeight)
734{
737 fTotalHistogram = new TH1D("total","total",nbins,xbins);
738 fPassedHistogram = new TH1D("passed","passed",nbins,xbins);
739 TH1::AddDirectory(bStatus);
740
741 Build(name,title);
742}
743
744////////////////////////////////////////////////////////////////////////////////
745/// Create 1-dimensional TEfficiency object with fixed bins size.
746///
747/// Constructor creates two new and empty histograms with a fixed binning.
748///
749/// Input:
750///
751/// - `name`: the common part of the name for both histograms(no blanks)
752/// fTotalHistogram has name: name + "_total"
753/// fPassedHistogram has name: name + "_passed"
754/// - `title`: the common part of the title for both histogram
755/// fTotalHistogram has title: title + " (total)"
756/// fPassedHistogram has title: title + " (passed)"
757/// It is possible to label the axis by passing a title with
758/// the following format: "title;xlabel;ylabel".
759/// - `nbinsx`: number of bins on the x-axis
760/// - `xlow`: lower edge of first bin
761/// - `xup`: upper edge of last bin
762
763TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
764 Double_t xlow,Double_t xup):
765fBeta_alpha(kDefBetaAlpha),
766fBeta_beta(kDefBetaBeta),
767fConfLevel(kDefConfLevel),
768fDirectory(0),
769fFunctions(0),
770fPaintGraph(0),
771fPaintHisto(0),
772fWeight(kDefWeight)
773{
776 fTotalHistogram = new TH1D("total","total",nbinsx,xlow,xup);
777 fPassedHistogram = new TH1D("passed","passed",nbinsx,xlow,xup);
778 TH1::AddDirectory(bStatus);
779
780 Build(name,title);
781}
782
783////////////////////////////////////////////////////////////////////////////////
784/// Create 2-dimensional TEfficiency object with fixed bin size.
785///
786/// Constructor creates two new and empty histograms with a fixed binning.
787///
788/// Input:
789///
790/// - `name`: the common part of the name for both histograms(no blanks)
791/// fTotalHistogram has name: name + "_total"
792/// fPassedHistogram has name: name + "_passed"
793/// - `title`: the common part of the title for both histogram
794/// fTotalHistogram has title: title + " (total)"
795/// fPassedHistogram has title: title + " (passed)"
796/// It is possible to label the axis by passing a title with
797/// the following format: "title;xlabel;ylabel;zlabel".
798/// - `nbinsx`: number of bins on the x-axis
799/// - `xlow`: lower edge of first x-bin
800/// - `xup`: upper edge of last x-bin
801/// - `nbinsy`: number of bins on the y-axis
802/// - `ylow`: lower edge of first y-bin
803/// - `yup`: upper edge of last y-bin
804
805TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
806 Double_t xlow,Double_t xup,Int_t nbinsy,
807 Double_t ylow,Double_t yup):
808fBeta_alpha(kDefBetaAlpha),
809fBeta_beta(kDefBetaBeta),
810fConfLevel(kDefConfLevel),
811fDirectory(0),
812fFunctions(0),
813fPaintGraph(0),
814fPaintHisto(0),
815fWeight(kDefWeight)
816{
819 fTotalHistogram = new TH2D("total","total",nbinsx,xlow,xup,nbinsy,ylow,yup);
820 fPassedHistogram = new TH2D("passed","passed",nbinsx,xlow,xup,nbinsy,ylow,yup);
821 TH1::AddDirectory(bStatus);
822
823 Build(name,title);
824}
825
826////////////////////////////////////////////////////////////////////////////////
827/// Create 2-dimensional TEfficiency object with variable bin size.
828///
829/// Constructor creates two new and empty histograms with a given binning.
830///
831/// Input:
832///
833/// - `name`: the common part of the name for both histograms(no blanks)
834/// fTotalHistogram has name: name + "_total"
835/// fPassedHistogram has name: name + "_passed"
836/// - `title`: the common part of the title for both histogram
837/// fTotalHistogram has title: title + " (total)"
838/// fPassedHistogram has title: title + " (passed)"
839/// It is possible to label the axis by passing a title with
840/// the following format: "title;xlabel;ylabel;zlabel".
841/// - `nbinsx`: number of bins on the x-axis
842/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
843/// xbins[nbinsx] ... lower edge for overflow x-bin
844/// - `nbinsy`: number of bins on the y-axis
845/// - `ybins`: array of length (nbins + 1) with low-edges for each bin
846/// ybins[nbinsy] ... lower edge for overflow y-bin
847
848TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
849 const Double_t* xbins,Int_t nbinsy,
850 const Double_t* ybins):
851fBeta_alpha(kDefBetaAlpha),
852fBeta_beta(kDefBetaBeta),
853fConfLevel(kDefConfLevel),
854fDirectory(0),
855fFunctions(0),
856fPaintGraph(0),
857fPaintHisto(0),
858fWeight(kDefWeight)
859{
862 fTotalHistogram = new TH2D("total","total",nbinsx,xbins,nbinsy,ybins);
863 fPassedHistogram = new TH2D("passed","passed",nbinsx,xbins,nbinsy,ybins);
864 TH1::AddDirectory(bStatus);
865
866 Build(name,title);
867}
868
869////////////////////////////////////////////////////////////////////////////////
870/// Create 3-dimensional TEfficiency object with fixed bin size.
871///
872/// Constructor creates two new and empty histograms with a fixed binning.
873///
874/// Input:
875///
876/// - `name`: the common part of the name for both histograms(no blanks)
877/// fTotalHistogram has name: name + "_total"
878/// fPassedHistogram has name: name + "_passed"
879/// - `title`: the common part of the title for both histogram
880/// fTotalHistogram has title: title + " (total)"
881/// fPassedHistogram has title: title + " (passed)"
882/// It is possible to label the axis by passing a title with
883/// the following format: "title;xlabel;ylabel;zlabel".
884/// - `nbinsx`: number of bins on the x-axis
885/// - `xlow`: lower edge of first x-bin
886/// - `xup`: upper edge of last x-bin
887/// - `nbinsy`: number of bins on the y-axis
888/// - `ylow`: lower edge of first y-bin
889/// - `yup`: upper edge of last y-bin
890/// - `nbinsz`: number of bins on the z-axis
891/// - `zlow`: lower edge of first z-bin
892/// - `zup`: upper edge of last z-bin
893
894TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
895 Double_t xlow,Double_t xup,Int_t nbinsy,
896 Double_t ylow,Double_t yup,Int_t nbinsz,
897 Double_t zlow,Double_t zup):
898fBeta_alpha(kDefBetaAlpha),
899fBeta_beta(kDefBetaBeta),
900fConfLevel(kDefConfLevel),
901fDirectory(0),
902fFunctions(0),
903fPaintGraph(0),
904fPaintHisto(0),
905fWeight(kDefWeight)
906{
909 fTotalHistogram = new TH3D("total","total",nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup);
910 fPassedHistogram = new TH3D("passed","passed",nbinsx,xlow,xup,nbinsy,ylow,yup,nbinsz,zlow,zup);
911 TH1::AddDirectory(bStatus);
912
913 Build(name,title);
914}
915
916////////////////////////////////////////////////////////////////////////////////
917/// Create 3-dimensional TEfficiency object with variable bin size.
918///
919/// Constructor creates two new and empty histograms with a given binning.
920///
921/// Input:
922///
923/// - `name`: the common part of the name for both histograms(no blanks)
924/// fTotalHistogram has name: name + "_total"
925/// fPassedHistogram has name: name + "_passed"
926/// - `title`: the common part of the title for both histogram
927/// fTotalHistogram has title: title + " (total)"
928/// fPassedHistogram has title: title + " (passed)"
929/// It is possible to label the axis by passing a title with
930/// the following format: "title;xlabel;ylabel;zlabel".
931/// - `nbinsx`: number of bins on the x-axis
932/// - `xbins`: array of length (nbins + 1) with low-edges for each bin
933/// xbins[nbinsx] ... lower edge for overflow x-bin
934/// - `nbinsy`: number of bins on the y-axis
935/// - `ybins`: array of length (nbins + 1) with low-edges for each bin
936/// xbins[nbinsx] ... lower edge for overflow y-bin
937/// - `nbinsz`: number of bins on the z-axis
938/// - `zbins`: array of length (nbins + 1) with low-edges for each bin
939/// xbins[nbinsx] ... lower edge for overflow z-bin
940
941TEfficiency::TEfficiency(const char* name,const char* title,Int_t nbinsx,
942 const Double_t* xbins,Int_t nbinsy,
943 const Double_t* ybins,Int_t nbinsz,
944 const Double_t* zbins):
945fBeta_alpha(kDefBetaAlpha),
946fBeta_beta(kDefBetaBeta),
947fConfLevel(kDefConfLevel),
948fDirectory(0),
949fFunctions(0),
950fPaintGraph(0),
951fPaintHisto(0),
952fWeight(kDefWeight)
953{
956 fTotalHistogram = new TH3D("total","total",nbinsx,xbins,nbinsy,ybins,nbinsz,zbins);
957 fPassedHistogram = new TH3D("passed","passed",nbinsx,xbins,nbinsy,ybins,nbinsz,zbins);
958 TH1::AddDirectory(bStatus);
959
960 Build(name,title);
961}
962
963////////////////////////////////////////////////////////////////////////////////
964/// Copy constructor.
965///
966///The list of associated objects (e.g. fitted functions) is not copied.
967///
968///Note:
969///
970/// - SetName(rEff.GetName() + "_copy") is called to set the names of the
971/// object and the histograms.
972/// - The titles are set by calling SetTitle("[copy] " + rEff.GetTitle()).
973/// - The copied TEfficiency object is NOT appended to a directory. It
974/// will not be written to disk during the next TFile::Write() command
975/// in order to prevent duplication of data. If you want to save this
976/// TEfficiency object anyway, you can either append it to a directory
977/// by calling SetDirectory(TDirectory*) or write it explicitly to disk
978/// by calling Write().
979
981TNamed(),
982TAttLine(),
983TAttFill(),
984TAttMarker(),
985fBeta_alpha(rEff.fBeta_alpha),
986fBeta_beta(rEff.fBeta_beta),
987fBeta_bin_params(rEff.fBeta_bin_params),
988fConfLevel(rEff.fConfLevel),
989fDirectory(0),
990fFunctions(0),
991fPaintGraph(0),
992fPaintHisto(0),
993fWeight(rEff.fWeight)
994{
995 // copy TObject bits
996 ((TObject&)rEff).Copy(*this);
997
1000 fTotalHistogram = (TH1*)((rEff.fTotalHistogram)->Clone());
1001 fPassedHistogram = (TH1*)((rEff.fPassedHistogram)->Clone());
1002 TH1::AddDirectory(bStatus);
1003
1004 TString name = rEff.GetName();
1005 name += "_copy";
1006 SetName(name);
1007 TString title = "[copy] ";
1008 title += rEff.GetTitle();
1009 SetTitle(title);
1010
1012
1013 SetDirectory(0);
1014
1015 //copy style
1016 rEff.TAttLine::Copy(*this);
1017 rEff.TAttFill::Copy(*this);
1018 rEff.TAttMarker::Copy(*this);
1019}
1020
1021////////////////////////////////////////////////////////////////////////////////
1022///default destructor
1023
1025{
1026 //delete all function in fFunctions
1027 // use same logic as in TH1 destructor
1028 // (see TH1::~TH1 code in TH1.cxx)
1029 if(fFunctions) {
1031 TObject* obj = 0;
1032 while ((obj = fFunctions->First())) {
1033 while(fFunctions->Remove(obj)) { }
1034 if (!obj->TestBit(kNotDeleted)) {
1035 break;
1036 }
1037 delete obj;
1038 obj = 0;
1039 }
1040 delete fFunctions;
1041 fFunctions = 0;
1042 }
1043
1044 if(fDirectory)
1045 fDirectory->Remove(this);
1046
1047 delete fTotalHistogram;
1048 delete fPassedHistogram;
1049 delete fPaintGraph;
1050 delete fPaintHisto;
1051}
1052
1053////////////////////////////////////////////////////////////////////////////////
1054/**
1055 Calculates the boundaries for the frequentist Agresti-Coull interval
1056
1057 \param total number of total events
1058 \param passed 0 <= number of passed events <= total
1059 \param level confidence level
1060 \param bUpper true - upper boundary is returned
1061 false - lower boundary is returned
1062
1063
1064 \f{eqnarray*}{
1065 \alpha &=& 1 - \frac{level}{2} \\
1066 \kappa &=& \Phi^{-1}(1 - \alpha,1)\ ... normal\ quantile\ function\\
1067 mode &=& \frac{passed + \frac{\kappa^{2}}{2}}{total + \kappa^{2}}\\
1068 \Delta &=& \kappa * \sqrt{\frac{mode * (1 - mode)}{total + \kappa^{2}}}\\
1069 return &=& max(0,mode - \Delta)\ or\ min(1,mode + \Delta)
1070 \f}
1071
1072*/
1073
1075{
1076 Double_t alpha = (1.0 - level)/2;
1077 Double_t kappa = ROOT::Math::normal_quantile(1 - alpha,1);
1078
1079 Double_t mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa);
1080 Double_t delta = kappa * std::sqrt(mode * (1 - mode) / (total + kappa * kappa));
1081
1082 if(bUpper)
1083 return ((mode + delta) > 1) ? 1.0 : (mode + delta);
1084 else
1085 return ((mode - delta) < 0) ? 0.0 : (mode - delta);
1086}
1087
1088////////////////////////////////////////////////////////////////////////////////
1089/// Calculates the boundaries for the frequentist Feldman-Cousins interval
1090///
1091/// \param total number of total events
1092/// \param passed 0 <= number of passed events <= total
1093/// \param level confidence level
1094/// \param bUpper: true - upper boundary is returned
1095/// false - lower boundary is returned
1096
1098{
1099 Double_t lower = 0;
1100 Double_t upper = 1;
1101 if (!FeldmanCousinsInterval(total,passed,level, lower, upper)) {
1102 ::Error("FeldmanCousins","Error running FC method - return 0 or 1");
1103 }
1104 return (bUpper) ? upper : lower;
1105}
1106
1107////////////////////////////////////////////////////////////////////////////////
1108/// Calculates the interval boundaries using the frequentist methods of Feldman-Cousins
1109///
1110/// \param[in] total number of total events
1111/// \param[in] passed 0 <= number of passed events <= total
1112/// \param[in] level confidence level
1113/// \param[out] lower lower boundary returned on exit
1114/// \param[out] upper lower boundary returned on exit
1115/// \return a flag with the status of the calculation
1116///
1117/// Calculation:
1118///
1119/// The Feldman-Cousins is a frequentist method where the interval is estimated using a Neyman construction where the ordering
1120/// is based on the likelihood ratio:
1121/// \f[
1122/// LR = \frac{Binomial(k | N, \epsilon)}{Binomial(k | N, \hat{\epsilon} ) }
1123/// \f]
1124/// See G. J. Feldman and R. D. Cousins, Phys. Rev. D57 (1998) 3873
1125/// and R. D. Cousins, K. E. Hymes, J. Tucker, Nuclear Instruments and Methods in Physics Research A 612 (2010) 388
1126///
1127/// Implemented using classes developed by Jordan Tucker and Luca Lista
1128/// See File hist/hist/src/TEfficiencyHelper.h
1129
1131{
1133 double alpha = 1.-level;
1134 fc.Init(alpha);
1135 fc.Calculate(passed, total);
1136 lower = fc.Lower();
1137 upper = fc.Upper();
1138 return true;
1139}
1140
1141////////////////////////////////////////////////////////////////////////////////
1142/// Calculates the boundaries using the mid-P binomial
1143/// interval (Lancaster method) from B. Cousing and J. Tucker.
1144/// See http://arxiv.org/abs/0905.3831 for a description and references for the method
1145///
1146/// Modify equal_tailed to get the kind of interval you want.
1147/// Can also be converted to interval on ratio of poisson means X/Y by the substitutions
1148/// ~~~ {.cpp}
1149/// X = passed
1150/// total = X + Y
1151/// lower_poisson = lower/(1 - lower)
1152/// upper_poisson = upper/(1 - upper)
1153/// ~~~
1154
1156{
1157 const double alpha = 1. - level;
1158 const bool equal_tailed = true; // change if you don;t want equal tailed interval
1159 const double alpha_min = equal_tailed ? alpha/2 : alpha;
1160 const double tol = 1e-9; // tolerance
1161 double pmin = 0;
1162 double pmax = 0;
1163 double p = 0;
1164
1165 pmin = 0; pmax = 1;
1166
1167
1168 // treat special case for 0<passed<1
1169 // do a linear interpolation of the upper limit values
1170 if ( passed > 0 && passed < 1) {
1171 double p0 = MidPInterval(total,0.0,level,bUpper);
1172 double p1 = MidPInterval(total,1.0,level,bUpper);
1173 p = (p1 - p0) * passed + p0;
1174 return p;
1175 }
1176
1177 while (std::abs(pmax - pmin) > tol) {
1178 p = (pmin + pmax)/2;
1179 //double v = 0.5 * ROOT::Math::binomial_pdf(int(passed), p, int(total));
1180 // make it work for non integer using the binomial - beta relationship
1181 double v = 0.5 * ROOT::Math::beta_pdf(p, passed+1., total-passed+1)/(total+1);
1182 //if (passed > 0) v += ROOT::Math::binomial_cdf(int(passed - 1), p, int(total));
1183 // compute the binomial cdf at passed -1
1184 if ( (passed-1) >= 0) v += ROOT::Math::beta_cdf_c(p, passed, total-passed+1);
1185
1186 double vmin = (bUpper) ? alpha_min : 1.- alpha_min;
1187 if (v > vmin)
1188 pmin = p;
1189 else
1190 pmax = p;
1191 }
1192
1193 return p;
1194}
1195
1196
1197////////////////////////////////////////////////////////////////////////////////
1198/**
1199Calculates the boundaries for a Bayesian confidence interval (shortest or central interval depending on the option)
1200
1201\param[in] total number of total events
1202\param[in] passed 0 <= number of passed events <= total
1203\param[in] level confidence level
1204\param[in] alpha shape parameter > 0 for the prior distribution (fBeta_alpha)
1205\param[in] beta shape parameter > 0 for the prior distribution (fBeta_beta)
1206\param[in] bUpper
1207 - true - upper boundary is returned
1208 - false - lower boundary is returned
1209\param[in] bShortest ??
1210
1211Note: In the case central confidence interval is calculated.
1212 when passed = 0 (or passed = total) the lower (or upper)
1213 interval values will be larger than 0 (or smaller than 1).
1214
1215Calculation:
1216
1217The posterior probability in bayesian statistics is given by:
1218\f[
1219 P(\varepsilon |k,N) \propto L(\varepsilon|k,N) \times Prior(\varepsilon)
1220\f]
1221As an efficiency can be interpreted as probability of a positive outcome of
1222a Bernoullli trial the likelihood function is given by the binomial
1223distribution:
1224\f[
1225 L(\varepsilon|k,N) = Binomial(N,k) \varepsilon ^{k} (1 - \varepsilon)^{N-k}
1226\f]
1227At the moment only beta distributions are supported as prior probabilities
1228of the efficiency (\f$ B(\alpha,\beta)\f$ is the beta function):
1229\f[
1230 Prior(\varepsilon) = \frac{1}{B(\alpha,\beta)} \varepsilon ^{\alpha - 1} (1 - \varepsilon)^{\beta - 1}
1231\f]
1232The posterior probability is therefore again given by a beta distribution:
1233\f[
1234 P(\varepsilon |k,N) \propto \varepsilon ^{k + \alpha - 1} (1 - \varepsilon)^{N - k + \beta - 1}
1235\f]
1236In case of central intervals
1237the lower boundary for the equal-tailed confidence interval is given by the
1238inverse cumulative (= quantile) function for the quantile \f$ \frac{1 - level}{2} \f$.
1239The upper boundary for the equal-tailed confidence interval is given by the
1240inverse cumulative (= quantile) function for the quantile \f$ \frac{1 + level}{2} \f$.
1241Hence it is the solution \f$ \varepsilon \f$ of the following equation:
1242\f[
1243 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}
1244\f]
1245In the case of shortest interval the minimum interval around the mode is found by minimizing the length of all intervals width the
1246given probability content. See TEfficiency::BetaShortestInterval
1247*/
1248
1250{
1251 Double_t a = double(passed)+alpha;
1252 Double_t b = double(total-passed)+beta;
1253
1254 if (bShortest) {
1255 double lower = 0;
1256 double upper = 1;
1257 BetaShortestInterval(level,a,b,lower,upper);
1258 return (bUpper) ? upper : lower;
1259 }
1260 else
1261 return BetaCentralInterval(level, a, b, bUpper);
1262}
1263
1264////////////////////////////////////////////////////////////////////////////////
1265/// Calculates the boundaries for a central confidence interval for a Beta distribution
1266///
1267/// \param[in] level confidence level
1268/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1269/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1270/// \param[in] bUpper true - upper boundary is returned
1271/// false - lower boundary is returned
1272
1274{
1275 if(bUpper) {
1276 if((a > 0) && (b > 0))
1277 return ROOT::Math::beta_quantile((1+level)/2,a,b);
1278 else {
1279 gROOT->Error("TEfficiency::BayesianCentral","Invalid input parameters - return 1");
1280 return 1;
1281 }
1282 }
1283 else {
1284 if((a > 0) && (b > 0))
1285 return ROOT::Math::beta_quantile((1-level)/2,a,b);
1286 else {
1287 gROOT->Error("TEfficiency::BayesianCentral","Invalid input parameters - return 0");
1288 return 0;
1289 }
1290 }
1291}
1292
1293struct Beta_interval_length {
1294 Beta_interval_length(Double_t level,Double_t alpha,Double_t beta ) :
1295 fCL(level), fAlpha(alpha), fBeta(beta)
1296 {}
1297
1298 Double_t LowerMax() {
1299 // max allowed value of lower given the interval size
1300 return ROOT::Math::beta_quantile_c(fCL, fAlpha,fBeta);
1301 }
1302
1303 Double_t operator() (double lower) const {
1304 // return length of interval
1305 Double_t plow = ROOT::Math::beta_cdf(lower, fAlpha, fBeta);
1306 Double_t pup = plow + fCL;
1307 double upper = ROOT::Math::beta_quantile(pup, fAlpha,fBeta);
1308 return upper-lower;
1309 }
1310 Double_t fCL; // interval size (confidence level)
1311 Double_t fAlpha; // beta distribution alpha parameter
1312 Double_t fBeta; // beta distribution beta parameter
1313
1314};
1315
1316////////////////////////////////////////////////////////////////////////////////
1317/// Calculates the boundaries for a shortest confidence interval for a Beta distribution
1318///
1319/// \param[in] level confidence level
1320/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1321/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1322/// \param[out] upper upper boundary is returned
1323/// \param[out] lower lower boundary is returned
1324///
1325/// The lower/upper boundary are then obtained by finding the shortest interval of the beta distribution
1326/// contained the desired probability level.
1327/// The length of all possible intervals is minimized in order to find the shortest one
1328
1330{
1331 if (a <= 0 || b <= 0) {
1332 lower = 0; upper = 1;
1333 gROOT->Error("TEfficiency::BayesianShortest","Invalid input parameters - return [0,1]");
1334 return kFALSE;
1335 }
1336
1337 // treat here special cases when mode == 0 or 1
1338 double mode = BetaMode(a,b);
1339 if (mode == 0.0) {
1340 lower = 0;
1341 upper = ROOT::Math::beta_quantile(level, a, b);
1342 return kTRUE;
1343 }
1344 if (mode == 1.0) {
1345 lower = ROOT::Math::beta_quantile_c(level, a, b);
1346 upper = 1.0;
1347 return kTRUE;
1348 }
1349 // special case when the shortest interval is undefined return the central interval
1350 // can happen for a posterior when passed=total=0
1351 //
1352 if ( a==b && a<=1.0) {
1353 lower = BetaCentralInterval(level,a,b,kFALSE);
1354 upper = BetaCentralInterval(level,a,b,kTRUE);
1355 return kTRUE;
1356 }
1357
1358 // for the other case perform a minimization
1359 // make a function of the length of the posterior interval as a function of lower bound
1360 Beta_interval_length intervalLength(level,a,b);
1361 // minimize the interval length
1364 minim.SetFunction(func, 0, intervalLength.LowerMax() );
1365 minim.SetNpx(2); // no need to bracket with many iterations. Just do few times to estimate some better points
1366 bool ret = minim.Minimize(100, 1.E-10,1.E-10);
1367 if (!ret) {
1368 gROOT->Error("TEfficiency::BayesianShortes","Error finding the shortest interval");
1369 return kFALSE;
1370 }
1371 lower = minim.XMinimum();
1372 upper = lower + minim.FValMinimum();
1373 return kTRUE;
1374}
1375
1376////////////////////////////////////////////////////////////////////////////////
1377/// Compute the mean (average) of the beta distribution
1378///
1379/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1380/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1381///
1382
1384{
1385 if (a <= 0 || b <= 0 ) {
1386 gROOT->Error("TEfficiency::BayesianMean","Invalid input parameters - return 0");
1387 return 0;
1388 }
1389
1390 Double_t mean = a / (a + b);
1391 return mean;
1392}
1393
1394////////////////////////////////////////////////////////////////////////////////
1395/// Compute the mode of the beta distribution
1396///
1397/// \param[in] a parameter > 0 for the beta distribution (for a posterior is passed + prior_alpha
1398/// \param[in] b parameter > 0 for the beta distribution (for a posterior is (total-passed) + prior_beta
1399///
1400/// note the mode is defined for a Beta(a,b) only if (a,b)>1 (a = passed+alpha; b = total-passed+beta)
1401/// return then the following in case (a,b) < 1:
1402/// - if (a==b) return 0.5 (it is really undefined)
1403/// - if (a < b) return 0;
1404/// - if (a > b) return 1;
1405
1407{
1408 if (a <= 0 || b <= 0 ) {
1409 gROOT->Error("TEfficiency::BayesianMode","Invalid input parameters - return 0");
1410 return 0;
1411 }
1412 if ( a <= 1 || b <= 1) {
1413 if ( a < b) return 0;
1414 if ( a > b) return 1;
1415 if (a == b) return 0.5; // cannot do otherwise
1416 }
1417
1418 // since a and b are > 1 here denominator cannot be 0 or < 0
1419 Double_t mode = (a - 1.0) / (a + b -2.0);
1420 return mode;
1421}
1422////////////////////////////////////////////////////////////////////////////////
1423/// Building standard data structure of a TEfficiency object
1424///
1425/// Notes:
1426/// - calls: SetName(name), SetTitle(title)
1427/// - set the statistic option to the default (kFCP)
1428/// - appends this object to the current directory SetDirectory(gDirectory)
1429
1430void TEfficiency::Build(const char* name,const char* title)
1431{
1432 SetName(name);
1433 SetTitle(title);
1434
1437
1438 SetBit(kPosteriorMode,false);
1440 SetBit(kUseWeights,false);
1441
1442 //set normalisation factors to 0, otherwise the += may not work properly
1445}
1446
1447////////////////////////////////////////////////////////////////////////////////
1448/// Checks binning for each axis
1449///
1450/// It is assumed that the passed histograms have the same dimension.
1451
1453{
1454
1455 const TAxis* ax1 = 0;
1456 const TAxis* ax2 = 0;
1457
1458 //check binning along axis
1459 for(Int_t j = 0; j < pass.GetDimension(); ++j) {
1460 switch(j) {
1461 case 0:
1462 ax1 = pass.GetXaxis();
1463 ax2 = total.GetXaxis();
1464 break;
1465 case 1:
1466 ax1 = pass.GetYaxis();
1467 ax2 = total.GetYaxis();
1468 break;
1469 case 2:
1470 ax1 = pass.GetZaxis();
1471 ax2 = total.GetZaxis();
1472 break;
1473 }
1474
1475 if(ax1->GetNbins() != ax2->GetNbins()) {
1476 gROOT->Info("TEfficiency::CheckBinning","Histograms are not consistent: they have different number of bins");
1477 return false;
1478 }
1479
1480 for(Int_t i = 1; i <= ax1->GetNbins() + 1; ++i)
1481 if(!TMath::AreEqualRel(ax1->GetBinLowEdge(i), ax2->GetBinLowEdge(i), 1.E-15)) {
1482 gROOT->Info("TEfficiency::CheckBinning","Histograms are not consistent: they have different bin edges");
1483 return false;
1484 }
1485
1486
1487 }
1488
1489 return true;
1490}
1491
1492////////////////////////////////////////////////////////////////////////////////
1493/// Checks the consistence of the given histograms
1494///
1495/// The histograms are considered as consistent if:
1496/// - both have the same dimension
1497/// - both have the same binning
1498/// - pass.GetBinContent(i) <= total.GetBinContent(i) for each bin i
1499///
1500
1502{
1503 if(pass.GetDimension() != total.GetDimension()) {
1504 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects have different dimensions");
1505 return false;
1506 }
1507
1508 if(!CheckBinning(pass,total)) {
1509 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects have different binning");
1510 return false;
1511 }
1512
1513 if(!CheckEntries(pass,total)) {
1514 gROOT->Error("TEfficiency::CheckConsistency","passed TEfficiency objects do not have consistent bin contents");
1515 return false;
1516 }
1517
1518 return true;
1519}
1520
1521////////////////////////////////////////////////////////////////////////////////
1522/// Checks whether bin contents are compatible with binomial statistics
1523///
1524/// The following inequality has to be valid for each bin i:
1525/// total.GetBinContent(i) >= pass.GetBinContent(i)
1526///
1527///
1528///
1529/// Note:
1530///
1531/// - It is assumed that both histograms have the same dimension and binning.
1532
1534{
1535
1536 //check: pass <= total
1537 Int_t nbinsx, nbinsy, nbinsz, nbins;
1538
1539 nbinsx = pass.GetNbinsX();
1540 nbinsy = pass.GetNbinsY();
1541 nbinsz = pass.GetNbinsZ();
1542
1543 switch(pass.GetDimension()) {
1544 case 1: nbins = nbinsx + 2; break;
1545 case 2: nbins = (nbinsx + 2) * (nbinsy + 2); break;
1546 case 3: nbins = (nbinsx + 2) * (nbinsy + 2) * (nbinsz + 2); break;
1547 default: nbins = 0;
1548 }
1549
1550 for(Int_t i = 0; i < nbins; ++i) {
1551 if(pass.GetBinContent(i) > total.GetBinContent(i)) {
1552 gROOT->Info("TEfficiency::CheckEntries","Histograms are not consistent: passed bin content > total bin content");
1553 return false;
1554 }
1555 }
1556
1557 return true;
1558}
1559
1560////////////////////////////////////////////////////////////////////////////////
1561/// Check if both histogram are weighted. If they are weighted a true is returned
1562///
1564{
1565 if (pass.GetSumw2N() == 0 && total.GetSumw2N() == 0) return false;
1566
1567 // check also that the total sum of weight and weight squares are consistent
1568 Double_t statpass[TH1::kNstat];
1569 Double_t stattotal[TH1::kNstat];
1570
1571 pass.GetStats(statpass);
1572 total.GetStats(stattotal);
1573
1574 double tolerance = (total.IsA() == TH1F::Class() ) ? 1.E-5 : 1.E-12;
1575
1576 //require: sum of weights == sum of weights^2
1577 if(!TMath::AreEqualRel(statpass[0],statpass[1],tolerance) ||
1578 !TMath::AreEqualRel(stattotal[0],stattotal[1],tolerance) ) {
1579 return true;
1580 }
1581
1582 // histograms are not weighted
1583 return false;
1584
1585}
1586
1587
1588////////////////////////////////////////////////////////////////////////////////
1589/// Create the graph used be painted (for dim=1 TEfficiency)
1590/// The return object is managed by the caller
1591
1593{
1594 if (GetDimension() != 1) {
1595 Error("CreatePaintingGraph","Call this function only for dimension == 1");
1596 return 0;
1597 }
1598
1599
1600 Int_t npoints = fTotalHistogram->GetNbinsX();
1601 TGraphAsymmErrors * graph = new TGraphAsymmErrors(npoints);
1602 graph->SetName("eff_graph");
1603 FillGraph(graph,opt);
1604
1605 return graph;
1606}
1607
1608
1609////////////////////////////////////////////////////////////////////////////////
1610/// Fill the graph to be painted with information from TEfficiency
1611/// Internal method called by TEfficiency::Paint or TEfficiency::CreateGraph
1612
1614{
1615 TString option = opt;
1616 option.ToLower();
1617
1618 Bool_t plot0Bins = false;
1619 if (option.Contains("e0") ) plot0Bins = true;
1620
1621 Double_t x,y,xlow,xup,ylow,yup;
1622 //point i corresponds to bin i+1 in histogram
1623 // point j is point graph index
1624 // LM: cannot use TGraph::SetPoint because it deletes the underlying
1625 // histogram each time (see TGraph::SetPoint)
1626 // so use it only when extra points are added to the graph
1627 Int_t j = 0;
1628 double * px = graph->GetX();
1629 double * py = graph->GetY();
1630 double * exl = graph->GetEXlow();
1631 double * exh = graph->GetEXhigh();
1632 double * eyl = graph->GetEYlow();
1633 double * eyh = graph->GetEYhigh();
1634 Int_t npoints = fTotalHistogram->GetNbinsX();
1635 for (Int_t i = 0; i < npoints; ++i) {
1636 if (!plot0Bins && fTotalHistogram->GetBinContent(i+1) == 0 ) continue;
1638 y = GetEfficiency(i+1);
1640 xup = fTotalHistogram->GetBinWidth(i+1) - xlow;
1641 ylow = GetEfficiencyErrorLow(i+1);
1642 yup = GetEfficiencyErrorUp(i+1);
1643 // in the case the graph already existed and extra points have been added
1644 if (j >= graph->GetN() ) {
1645 graph->SetPoint(j,x,y);
1646 graph->SetPointError(j,xlow,xup,ylow,yup);
1647 }
1648 else {
1649 px[j] = x;
1650 py[j] = y;
1651 exl[j] = xlow;
1652 exh[j] = xup;
1653 eyl[j] = ylow;
1654 eyh[j] = yup;
1655 }
1656 j++;
1657 }
1658
1659 // tell the graph the effective number of points
1660 graph->Set(j);
1661 //refresh title before painting if changed
1662 TString oldTitle = graph->GetTitle();
1663 TString newTitle = GetTitle();
1664 if (oldTitle != newTitle ) {
1665 graph->SetTitle(newTitle);
1666 }
1667
1668 // set the axis labels
1671 if (xlabel) graph->GetXaxis()->SetTitle(xlabel);
1672 if (ylabel) graph->GetYaxis()->SetTitle(ylabel);
1673
1674 //copying style information
1678
1679 // this method forces the graph to compute correctly the axis
1680 // according to the given points
1681 graph->GetHistogram();
1682
1683}
1684
1685////////////////////////////////////////////////////////////////////////////////
1686/// Create the histogram used to be painted (for dim=2 TEfficiency)
1687/// The return object is managed by the caller
1688
1690{
1691 if (GetDimension() != 2) {
1692 Error("CreatePaintingistogram","Call this function only for dimension == 2");
1693 return 0;
1694 }
1695
1696 Int_t nbinsx = fTotalHistogram->GetNbinsX();
1697 Int_t nbinsy = fTotalHistogram->GetNbinsY();
1698 TAxis * xaxis = fTotalHistogram->GetXaxis();
1699 TAxis * yaxis = fTotalHistogram->GetYaxis();
1700 TH2 * hist = 0;
1701
1702 if (xaxis->IsVariableBinSize() && yaxis->IsVariableBinSize() )
1703 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXbins()->GetArray(),
1704 nbinsy,yaxis->GetXbins()->GetArray());
1705 else if (xaxis->IsVariableBinSize() && ! yaxis->IsVariableBinSize() )
1706 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXbins()->GetArray(),
1707 nbinsy,yaxis->GetXmin(), yaxis->GetXmax());
1708 else if (!xaxis->IsVariableBinSize() && yaxis->IsVariableBinSize() )
1709 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXmin(), xaxis->GetXmax(),
1710 nbinsy,yaxis->GetXbins()->GetArray());
1711 else
1712 hist = new TH2F("eff_histo",GetTitle(),nbinsx,xaxis->GetXmin(), xaxis->GetXmax(),
1713 nbinsy,yaxis->GetXmin(), yaxis->GetXmax());
1714
1715
1716 hist->SetDirectory(0);
1717
1718 FillHistogram(hist);
1719
1720 return hist;
1721}
1722
1723////////////////////////////////////////////////////////////////////////////////
1724/// Fill the 2d histogram to be painted with information from TEfficiency 2D
1725/// Internal method called by TEfficiency::Paint or TEfficiency::CreatePaintingGraph
1726
1728{
1729 //refresh title before each painting
1730 hist->SetTitle(GetTitle());
1731
1732 // set the axis labels
1736 if (xlabel) hist->GetXaxis()->SetTitle(xlabel);
1737 if (ylabel) hist->GetYaxis()->SetTitle(ylabel);
1738 if (zlabel) hist->GetZaxis()->SetTitle(zlabel);
1739
1740 Int_t bin;
1741 Int_t nbinsx = hist->GetNbinsX();
1742 Int_t nbinsy = hist->GetNbinsY();
1743 for(Int_t i = 0; i < nbinsx + 2; ++i) {
1744 for(Int_t j = 0; j < nbinsy + 2; ++j) {
1745 bin = GetGlobalBin(i,j);
1746 hist->SetBinContent(bin,GetEfficiency(bin));
1747 }
1748 }
1749
1750 //copying style information
1751 TAttLine::Copy(*hist);
1752 TAttFill::Copy(*hist);
1753 TAttMarker::Copy(*hist);
1754 hist->SetStats(0);
1755
1756 return;
1757
1758}
1759////////////////////////////////////////////////////////////////////////////////
1760/**
1761Calculates the boundaries for the frequentist Clopper-Pearson interval
1762
1763This interval is recommended by the PDG.
1764
1765\param[in] total number of total events
1766\param[in] passed 0 <= number of passed events <= total
1767\param[in] level confidence level
1768\param[in] bUpper true - upper boundary is returned
1769 ;false - lower boundary is returned
1770
1771Calculation:
1772
1773The lower boundary of the Clopper-Pearson interval is the "exact" inversion
1774of the test:
1775 \f{eqnarray*}{
1776 P(x \geq passed; total) &=& \frac{1 - level}{2}\\
1777 P(x \geq passed; total) &=& 1 - P(x \leq passed - 1; total)\\
1778 &=& 1 - \frac{1}{norm} * \int_{0}^{1 - \varepsilon} t^{total - passed} (1 - t)^{passed - 1} dt\\
1779 &=& 1 - \frac{1}{norm} * \int_{\varepsilon}^{1} t^{passed - 1} (1 - t)^{total - passed} dt\\
1780 &=& \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed - 1} (1 - t)^{total - passed} dt\\
1781 &=& I_{\varepsilon}(passed,total - passed + 1)
1782 \f}
1783The lower boundary is therefore given by the \f$ \frac{1 - level}{2}\f$ quantile
1784of the beta distribution.
1785
1786The upper boundary of the Clopper-Pearson interval is the "exact" inversion
1787of the test:
1788 \f{eqnarray*}{
1789 P(x \leq passed; total) &=& \frac{1 - level}{2}\\
1790 P(x \leq passed; total) &=& \frac{1}{norm} * \int_{0}^{1 - \varepsilon} t^{total - passed - 1} (1 - t)^{passed} dt\\
1791 &=& \frac{1}{norm} * \int_{\varepsilon}^{1} t^{passed} (1 - t)^{total - passed - 1} dt\\
1792 &=& 1 - \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed} (1 - t)^{total - passed - 1} dt\\
1793 \Rightarrow 1 - \frac{1 - level}{2} &=& \frac{1}{norm} * \int_{0}^{\varepsilon} t^{passed} (1 - t)^{total - passed -1} dt\\
1794 \frac{1 + level}{2} &=& I_{\varepsilon}(passed + 1,total - passed)
1795 \f}
1796The upper boundary is therefore given by the \f$\frac{1 + level}{2}\f$ quantile
1797of the beta distribution.
1798
1799Note: The connection between the binomial distribution and the regularized
1800 incomplete beta function \f$ I_{\varepsilon}(\alpha,\beta)\f$ has been used.
1801*/
1802
1804{
1805 Double_t alpha = (1.0 - level) / 2;
1806 if(bUpper)
1807 return ((passed == total) ? 1.0 : ROOT::Math::beta_quantile(1 - alpha,passed + 1,total-passed));
1808 else
1809 return ((passed == 0) ? 0.0 : ROOT::Math::beta_quantile(alpha,passed,total-passed+1.0));
1810}
1811////////////////////////////////////////////////////////////////////////////////
1812/**
1813 Calculates the combined efficiency and its uncertainties
1814
1815 This method does a bayesian combination of the given samples.
1816
1817 \param[in] up contains the upper limit of the confidence interval afterwards
1818 \param[in] low contains the lower limit of the confidence interval afterwards
1819 \param[in] n number of samples which are combined
1820 \param[in] pass array of length n containing the number of passed events
1821 \param[in] total array of length n containing the corresponding numbers of total events
1822 \param[in] alpha shape parameters for the beta distribution as prior
1823 \param[in] beta shape parameters for the beta distribution as prior
1824 \param[in] level desired confidence level
1825 \param[in] w weights for each sample; if not given, all samples get the weight 1
1826 The weights do not need to be normalized, since they are internally renormalized
1827 to the number of effective entries.
1828 \param[in] opt
1829 - mode : The mode is returned instead of the mean of the posterior as best value
1830 When using the mode the shortest interval is also computed instead of the central one
1831 - shortest: compute shortest interval (done by default if mode option is set)
1832 - central: compute central interval (done by default if mode option is NOT set)
1833
1834 Calculation:
1835
1836 The combined posterior distributions is calculated from the Bayes theorem assuming a common prior Beta distribution.
1837 It is easy to proof that the combined posterior is then:
1838 \f{eqnarray*}{
1839 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)\\
1840 w_{i} &=& weight\ for\ each\ sample\ renormalized\ to\ the\ effective\ entries\\
1841 w^{'}_{i} &=& w_{i} \frac{ \sum_{i} {w_{i} } } { \sum_{i} {w_{i}^{2} } }
1842 \f}
1843
1844 The estimated efficiency is the mode (or the mean) of the obtained posterior distribution
1845
1846 The boundaries of the confidence interval for a confidence level (1 - a)
1847 are given by the a/2 and 1-a/2 quantiles of the resulting cumulative
1848 distribution.
1849
1850 Example (uniform prior distribution):
1851
1852Begin_Macro(source)
1853{
1854 TCanvas* c1 = new TCanvas("c1","",600,800);
1855 c1->Divide(1,2);
1856 c1->SetFillStyle(1001);
1857 c1->SetFillColor(kWhite);
1858
1859 TF1* p1 = new TF1("p1","TMath::BetaDist(x,19,9)",0,1);
1860 TF1* p2 = new TF1("p2","TMath::BetaDist(x,4,8)",0,1);
1861 TF1* comb = new TF1("comb2","TMath::BetaDist(x,[0],[1])",0,1);
1862 double nrm = 1./(0.6*0.6+0.4*0.4); // weight normalization
1863 double a = 0.6*18.0 + 0.4*3.0 + 1.0; // new alpha parameter of combined beta dist.
1864 double b = 0.6*10+0.4*7+1.0; // new beta parameter of combined beta dist.
1865 comb->SetParameters(nrm*a ,nrm *b );
1866 TF1* const1 = new TF1("const1","0.05",0,1);
1867 TF1* const2 = new TF1("const2","0.95",0,1);
1868
1869 p1->SetLineColor(kRed);
1870 p1->SetTitle("combined posteriors;#epsilon;P(#epsilon|k,N)");
1871 p2->SetLineColor(kBlue);
1872 comb->SetLineColor(kGreen+2);
1873
1874 TLegend* leg1 = new TLegend(0.12,0.65,0.5,0.85);
1875 leg1->AddEntry(p1,"k1 = 18, N1 = 26","l");
1876 leg1->AddEntry(p2,"k2 = 3, N2 = 10","l");
1877 leg1->AddEntry(comb,"combined: p1 = 0.6, p2=0.4","l");
1878
1879 c1->cd(1);
1880 comb->Draw();
1881 p1->Draw("same");
1882 p2->Draw("same");
1883 leg1->Draw("same");
1884 c1->cd(2);
1885 const1->SetLineWidth(1);
1886 const2->SetLineWidth(1);
1887 TGraph* gr = (TGraph*)comb->DrawIntegral();
1888 gr->SetTitle("cumulative function of combined posterior with boundaries for cl = 95%;#epsilon;CDF");
1889 const1->Draw("same");
1890 const2->Draw("same");
1891
1892 c1->cd(0);
1893 return c1;
1894}
1895End_Macro
1896
1897**/
1898////////////////////////////////////////////////////////////////////
1900 const Int_t* pass,const Int_t* total,
1901 Double_t alpha, Double_t beta,
1902 Double_t level,const Double_t* w,Option_t* opt)
1903{
1904 TString option(opt);
1905 option.ToLower();
1906
1907 //LM: new formula for combination
1908 // works only if alpha beta are the same always
1909 // the weights are normalized to w(i) -> N_eff w(i)/ Sum w(i)
1910 // i.e. w(i) -> Sum (w(i) / Sum (w(i)^2) * w(i)
1911 // norm = Sum (w(i) / Sum (w(i)^2)
1912 double ntot = 0;
1913 double ktot = 0;
1914 double sumw = 0;
1915 double sumw2 = 0;
1916 for (int i = 0; i < n ; ++i) {
1917 if(pass[i] > total[i]) {
1918 ::Error("TEfficiency::Combine","total events = %i < passed events %i",total[i],pass[i]);
1919 ::Info("TEfficiency::Combine","stop combining");
1920 return -1;
1921 }
1922
1923 ntot += w[i] * total[i];
1924 ktot += w[i] * pass[i];
1925 sumw += w[i];
1926 sumw2 += w[i]*w[i];
1927 //mean += w[i] * (pass[i] + alpha[i])/(total[i] + alpha[i] + beta[i]);
1928 }
1929 double norm = sumw/sumw2;
1930 ntot *= norm;
1931 ktot *= norm;
1932 if(ktot > ntot) {
1933 ::Error("TEfficiency::Combine","total = %f < passed %f",ntot,ktot);
1934 ::Info("TEfficiency::Combine","stop combining");
1935 return -1;
1936 }
1937
1938 double a = ktot + alpha;
1939 double b = ntot - ktot + beta;
1940
1941 double mean = a/(a+b);
1942 double mode = BetaMode(a,b);
1943
1944
1945 Bool_t shortestInterval = option.Contains("sh") || ( option.Contains("mode") && !option.Contains("cent") );
1946
1947 if (shortestInterval)
1948 BetaShortestInterval(level, a, b, low, up);
1949 else {
1950 low = BetaCentralInterval(level, a, b, false);
1951 up = BetaCentralInterval(level, a, b, true);
1952 }
1953
1954 if (option.Contains("mode")) return mode;
1955 return mean;
1956
1957}
1958////////////////////////////////////////////////////////////////////////////////
1959/// Combines a list of 1-dimensional TEfficiency objects
1960///
1961/// A TGraphAsymmErrors object is returned which contains the estimated
1962/// efficiency and its uncertainty for each bin.
1963/// If the combination fails, a zero pointer is returned.
1964///
1965/// At the moment the combining is only implemented for bayesian statistics.
1966///
1967/// \param[in] pList list containing TEfficiency objects which should be combined
1968/// only one-dimensional efficiencies are taken into account
1969/// \param[in] option
1970/// - s : strict combining; only TEfficiency objects with the same beta
1971/// prior and the flag kIsBayesian == true are combined
1972/// If not specified the prior parameter of the first TEfficiency object is used
1973/// - v : verbose mode; print information about combining
1974/// - cl=x : set confidence level (0 < cl < 1). If not specified, the
1975/// confidence level of the first TEfficiency object is used.
1976/// - mode Use mode of combined posterior as estimated value for the efficiency
1977/// - shortest: compute shortest interval (done by default if mode option is set)
1978/// - central: compute central interval (done by default if mode option is NOT set)
1979/// \param[in] n number of weights (has to be the number of one-dimensional
1980/// TEfficiency objects in pList)
1981/// If no weights are passed, the internal weights GetWeight() of
1982/// the given TEfficiency objects are used.
1983/// \param[in] w array of length n with weights for each TEfficiency object in
1984/// pList (w[0] correspond to pList->First ... w[n-1] -> pList->Last)
1985/// The weights do not have to be normalised.
1986///
1987/// For each bin the calculation is done by the Combine(double&, double& ...) method.
1988
1990 Int_t n,const Double_t* w)
1991{
1992 TString opt = option;
1993 opt.ToLower();
1994
1995 //parameter of prior distribution, confidence level and normalisation factor
1996 Double_t alpha = -1;
1997 Double_t beta = -1;
1998 Double_t level = 0;
1999
2000 //flags for combining
2001 Bool_t bStrict = false;
2002 Bool_t bOutput = false;
2003 Bool_t bWeights = false;
2004 //list of all information needed to weight and combine efficiencies
2005 std::vector<TH1*> vTotal; vTotal.reserve(n);
2006 std::vector<TH1*> vPassed; vPassed.reserve(n);
2007 std::vector<Double_t> vWeights; vWeights.reserve(n);
2008 // std::vector<Double_t> vAlpha;
2009 // std::vector<Double_t> vBeta;
2010
2011 if(opt.Contains("s")) {
2012 opt.ReplaceAll("s","");
2013 bStrict = true;
2014 }
2015
2016 if(opt.Contains("v")) {
2017 opt.ReplaceAll("v","");
2018 bOutput = true;
2019 }
2020
2021 if(opt.Contains("cl=")) {
2022 Ssiz_t pos = opt.Index("cl=") + 3;
2023 level = atof( opt(pos,opt.Length() ).Data() );
2024 if((level <= 0) || (level >= 1))
2025 level = 0;
2026 opt.ReplaceAll("cl=","");
2027 }
2028
2029 //are weights explicitly given
2030 if(n && w) {
2031 bWeights = true;
2032 for(Int_t k = 0; k < n; ++k) {
2033 if(w[k] > 0)
2034 vWeights.push_back(w[k]);
2035 else {
2036 gROOT->Error("TEfficiency::Combine","invalid custom weight found w = %.2lf",w[k]);
2037 gROOT->Info("TEfficiency::Combine","stop combining");
2038 return 0;
2039 }
2040 }
2041 }
2042
2043 TIter next(pList);
2044 TObject* obj = 0;
2045 TEfficiency* pEff = 0;
2046 while((obj = next())) {
2047 pEff = dynamic_cast<TEfficiency*>(obj);
2048 //is object a TEfficiency object?
2049 if(pEff) {
2050 if(pEff->GetDimension() > 1)
2051 continue;
2052 if(!level) level = pEff->GetConfidenceLevel();
2053
2054 if(alpha<1) alpha = pEff->GetBetaAlpha();
2055 if(beta<1) beta = pEff->GetBetaBeta();
2056
2057 //if strict combining, check priors, confidence level and statistic
2058 if(bStrict) {
2059 if(alpha != pEff->GetBetaAlpha())
2060 continue;
2061 if(beta != pEff->GetBetaBeta())
2062 continue;
2063 if(!pEff->UsesBayesianStat())
2064 continue;
2065 }
2066
2067 vTotal.push_back(pEff->fTotalHistogram);
2068 vPassed.push_back(pEff->fPassedHistogram);
2069
2070 //no weights given -> use weights of TEfficiency objects
2071 if(!bWeights)
2072 vWeights.push_back(pEff->fWeight);
2073
2074 //strict combining -> using global prior
2075 // if(bStrict) {
2076 // vAlpha.push_back(alpha);
2077 // vBeta.push_back(beta);
2078 // }
2079 // else {
2080 // vAlpha.push_back(pEff->GetBetaAlpha());
2081 // vBeta.push_back(pEff->GetBetaBeta());
2082 // }
2083 }
2084 }
2085
2086 //no TEfficiency objects found
2087 if(vTotal.empty()) {
2088 gROOT->Error("TEfficiency::Combine","no TEfficiency objects in given list");
2089 gROOT->Info("TEfficiency::Combine","stop combining");
2090 return 0;
2091 }
2092
2093 //invalid number of custom weights
2094 if(bWeights && (n != (Int_t)vTotal.size())) {
2095 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());
2096 gROOT->Info("TEfficiency::Combine","stop combining");
2097 return 0;
2098 }
2099
2100 Int_t nbins_max = vTotal.at(0)->GetNbinsX();
2101 //check binning of all histograms
2102 for(UInt_t i=0; i<vTotal.size(); ++i) {
2103 if (!TEfficiency::CheckBinning(*vTotal.at(0),*vTotal.at(i)) )
2104 gROOT->Warning("TEfficiency::Combine","histograms have not the same binning -> results may be useless");
2105 if(vTotal.at(i)->GetNbinsX() < nbins_max) nbins_max = vTotal.at(i)->GetNbinsX();
2106 }
2107
2108 //display information about combining
2109 if(bOutput) {
2110 gROOT->Info("TEfficiency::Combine","combining %i TEfficiency objects",(Int_t)vTotal.size());
2111 if(bWeights)
2112 gROOT->Info("TEfficiency::Combine","using custom weights");
2113 if(bStrict) {
2114 gROOT->Info("TEfficiency::Combine","using the following prior probability for the efficiency: P(e) ~ Beta(e,%.3lf,%.3lf)",alpha,beta);
2115 }
2116 else
2117 gROOT->Info("TEfficiency::Combine","using individual priors of each TEfficiency object");
2118 gROOT->Info("TEfficiency::Combine","confidence level = %.2lf",level);
2119 }
2120
2121 //create TGraphAsymmErrors with efficiency
2122 std::vector<Double_t> x(nbins_max);
2123 std::vector<Double_t> xlow(nbins_max);
2124 std::vector<Double_t> xhigh(nbins_max);
2125 std::vector<Double_t> eff(nbins_max);
2126 std::vector<Double_t> efflow(nbins_max);
2127 std::vector<Double_t> effhigh(nbins_max);
2128
2129 //parameters for combining:
2130 //number of objects
2131 Int_t num = vTotal.size();
2132 std::vector<Int_t> pass(num);
2133 std::vector<Int_t> total(num);
2134
2135 //loop over all bins
2136 Double_t low = 0;
2137 Double_t up = 0;
2138 for(Int_t i=1; i <= nbins_max; ++i) {
2139 //the binning of the x-axis is taken from the first total histogram
2140 x[i-1] = vTotal.at(0)->GetBinCenter(i);
2141 xlow[i-1] = x[i-1] - vTotal.at(0)->GetBinLowEdge(i);
2142 xhigh[i-1] = vTotal.at(0)->GetBinWidth(i) - xlow[i-1];
2143
2144 for(Int_t j = 0; j < num; ++j) {
2145 pass[j] = (Int_t)(vPassed.at(j)->GetBinContent(i) + 0.5);
2146 total[j] = (Int_t)(vTotal.at(j)->GetBinContent(i) + 0.5);
2147 }
2148
2149 //fill efficiency and errors
2150 eff[i-1] = Combine(up,low,num,&pass[0],&total[0],alpha,beta,level,&vWeights[0],opt.Data());
2151 //did an error occurred ?
2152 if(eff[i-1] == -1) {
2153 gROOT->Error("TEfficiency::Combine","error occurred during combining");
2154 gROOT->Info("TEfficiency::Combine","stop combining");
2155 return 0;
2156 }
2157 efflow[i-1]= eff[i-1] - low;
2158 effhigh[i-1]= up - eff[i-1];
2159 }//loop over all bins
2160
2161 TGraphAsymmErrors* gr = new TGraphAsymmErrors(nbins_max,&x[0],&eff[0],&xlow[0],&xhigh[0],&efflow[0],&effhigh[0]);
2162
2163 return gr;
2164}
2165
2166////////////////////////////////////////////////////////////////////////////////
2167/// Compute distance from point px,py to a graph.
2168///
2169/// Compute the closest distance of approach from point px,py to this line.
2170/// The distance is computed in pixels units.
2171///
2172/// Forward the call to the painted graph
2173
2175{
2176 if (fPaintGraph) return fPaintGraph->DistancetoPrimitive(px,py);
2177 if (fPaintHisto) return fPaintHisto->DistancetoPrimitive(px,py);
2178 return 0;
2179}
2180
2181
2182////////////////////////////////////////////////////////////////////////////////
2183/// Draws the current TEfficiency object
2184///
2185/// \param[in] opt
2186/// - 1-dimensional case: same options as TGraphAsymmErrors::Draw()
2187/// but as default "AP" is used
2188/// - 2-dimensional case: same options as TH2::Draw()
2189/// - 3-dimensional case: not yet supported
2190///
2191/// Specific TEfficiency drawing options:
2192/// - E0 - plot bins where the total number of passed events is zero
2193/// (the error interval will be [0,1] )
2194
2196{
2197 //check options
2198 TString option = opt;
2199 option.ToLower();
2200 // use by default "AP"
2201 if (option.IsNull() ) option = "ap";
2202
2203 if(gPad && !option.Contains("same"))
2204 gPad->Clear();
2205 else {
2206 // add always "a" if not present
2207 if (!option.Contains("a") ) option += "a";
2208 }
2209
2210 // add always p to the option
2211 if (!option.Contains("p") ) option += "p";
2212
2213
2214 AppendPad(option.Data());
2215}
2216
2217////////////////////////////////////////////////////////////////////////////////
2218/// Execute action corresponding to one event.
2219///
2220/// This member function is called when the drawn class is clicked with the locator
2221/// If Left button clicked on one of the line end points, this point
2222/// follows the cursor until button is released.
2223///
2224/// if Middle button clicked, the line is moved parallel to itself
2225/// until the button is released.
2226/// Forward the call to the underlying graph
2227
2229{
2230 if (fPaintGraph) fPaintGraph->ExecuteEvent(event,px,py);
2231 else if (fPaintHisto) fPaintHisto->ExecuteEvent(event,px,py);
2232}
2233
2234////////////////////////////////////////////////////////////////////////////////
2235/// This function is used for filling the two histograms.
2236///
2237/// \param[in] bPassed flag whether the current event passed the selection
2238/// - true: both histograms are filled
2239/// - false: only the total histogram is filled
2240/// \param[in] x x-value
2241/// \param[in] y y-value (use default=0 for 1-D efficiencies)
2242/// \param[in] z z-value (use default=0 for 2-D or 1-D efficiencies)
2243
2245{
2246 switch(GetDimension()) {
2247 case 1:
2249 if(bPassed)
2251 break;
2252 case 2:
2253 ((TH2*)(fTotalHistogram))->Fill(x,y);
2254 if(bPassed)
2255 ((TH2*)(fPassedHistogram))->Fill(x,y);
2256 break;
2257 case 3:
2258 ((TH3*)(fTotalHistogram))->Fill(x,y,z);
2259 if(bPassed)
2260 ((TH3*)(fPassedHistogram))->Fill(x,y,z);
2261 break;
2262 }
2263}
2264
2265////////////////////////////////////////////////////////////////////////////////
2266///This function is used for filling the two histograms with a weight.
2267///
2268/// \param[in] bPassed flag whether the current event passed the selection
2269/// - true: both histograms are filled
2270/// - false: only the total histogram is filled
2271/// \param[in] weight weight for the event
2272/// \param[in] x x-value
2273/// \param[in] y y-value (use default=0 for 1-D efficiencies)
2274/// \param[in] z z-value (use default=0 for 2-D or 1-D efficiencies)
2275///
2276/// Note: - this function will call SetUseWeightedEvents if it was not called by the user before
2277
2279{
2280 if(!TestBit(kUseWeights))
2281 {
2282 // Info("FillWeighted","call SetUseWeightedEvents() manually to ensure correct storage of sum of weights squared");
2284 }
2285
2286 switch(GetDimension()) {
2287 case 1:
2288 fTotalHistogram->Fill(x,weight);
2289 if(bPassed)
2290 fPassedHistogram->Fill(x,weight);
2291 break;
2292 case 2:
2293 ((TH2*)(fTotalHistogram))->Fill(x,y,weight);
2294 if(bPassed)
2295 ((TH2*)(fPassedHistogram))->Fill(x,y,weight);
2296 break;
2297 case 3:
2298 ((TH3*)(fTotalHistogram))->Fill(x,y,z,weight);
2299 if(bPassed)
2300 ((TH3*)(fPassedHistogram))->Fill(x,y,z,weight);
2301 break;
2302 }
2303}
2304
2305////////////////////////////////////////////////////////////////////////////////
2306/// Returns the global bin number containing the given values
2307///
2308/// Note:
2309///
2310/// - values which belong to dimensions higher than the current dimension
2311/// of the TEfficiency object are ignored (i.e. for 1-dimensional
2312/// efficiencies only the x-value is considered)
2313
2315{
2317 Int_t ny = 0;
2318 Int_t nz = 0;
2319
2320 switch(GetDimension()) {
2321 case 3: nz = fTotalHistogram->GetZaxis()->FindFixBin(z);
2322 case 2: ny = fTotalHistogram->GetYaxis()->FindFixBin(y);break;
2323 }
2324
2325 return GetGlobalBin(nx,ny,nz);
2326}
2327
2328////////////////////////////////////////////////////////////////////////////////
2329/// Fits the efficiency using the TBinomialEfficiencyFitter class
2330///
2331/// The resulting fit function is added to the list of associated functions.
2332///
2333/// Options:
2334/// - "+": previous fitted functions in the list are kept, by default
2335/// all functions in the list are deleted
2336/// - for more fitting options see TBinomialEfficiencyFitter::Fit
2337
2339{
2340 TString option = opt;
2341 option.ToLower();
2342
2343 //replace existing functions in list with same name
2344 Bool_t bDeleteOld = true;
2345 if(option.Contains("+")) {
2346 option.ReplaceAll("+","");
2347 bDeleteOld = false;
2348 }
2349
2351
2352 TFitResultPtr result = Fitter.Fit(f1,option.Data());
2353
2354 //create copy which is appended to the list
2355 TF1* pFunc = new TF1(*f1);
2356
2357 if(bDeleteOld) {
2358 TIter next(fFunctions);
2359 TObject* obj = 0;
2360 while((obj = next())) {
2361 if(obj->InheritsFrom(TF1::Class())) {
2362 fFunctions->Remove(obj);
2363 delete obj;
2364 }
2365 }
2366 }
2367
2368 // create list if necessary
2369 if(!fFunctions)
2370 fFunctions = new TList();
2371
2372 fFunctions->Add(pFunc);
2373
2374 return result;
2375}
2376
2377////////////////////////////////////////////////////////////////////////////////
2378/// Returns a cloned version of fPassedHistogram
2379///
2380/// Notes:
2381/// - The histogram is filled with unit weights. You might want to scale
2382/// it with the global weight GetWeight().
2383/// - The returned object is owned by the user who has to care about the
2384/// deletion of the new TH1 object.
2385/// - This histogram is by default NOT attached to the current directory
2386/// to avoid duplication of data. If you want to store it automatically
2387/// during the next TFile::Write() command, you have to attach it to
2388/// the corresponding directory.
2389///
2390/// ~~~~~~~{.cpp}
2391/// TFile* pFile = new TFile("passed.root","update");
2392/// TEfficiency* pEff = (TEfficiency*)gDirectory->Get("my_eff");
2393/// TH1* copy = pEff->GetCopyPassedHisto();
2394/// copy->SetDirectory(gDirectory);
2395/// pFile->Write();
2396/// ~~~~~~~
2397
2399{
2400 Bool_t bStatus = TH1::AddDirectoryStatus();
2402 TH1* tmp = (TH1*)(fPassedHistogram->Clone());
2403 TH1::AddDirectory(bStatus);
2404
2405 return tmp;
2406}
2407
2408////////////////////////////////////////////////////////////////////////////////
2409/// Returns a cloned version of fTotalHistogram
2410///
2411/// Notes:
2412/// - The histogram is filled with unit weights. You might want to scale
2413/// it with the global weight GetWeight().
2414/// - The returned object is owned by the user who has to care about the
2415/// deletion of the new TH1 object.
2416/// - This histogram is by default NOT attached to the current directory
2417/// to avoid duplication of data. If you want to store it automatically
2418/// during the next TFile::Write() command, you have to attach it to
2419/// the corresponding directory.
2420///
2421/// ~~~~~~~{.cpp}
2422/// TFile* pFile = new TFile("total.root","update");
2423/// TEfficiency* pEff = (TEfficiency*)gDirectory->Get("my_eff");
2424/// TH1* copy = pEff->GetCopyTotalHisto();
2425/// copy->SetDirectory(gDirectory);
2426/// pFile->Write();
2427/// ~~~~~~~
2428
2430{
2431 Bool_t bStatus = TH1::AddDirectoryStatus();
2433 TH1* tmp = (TH1*)(fTotalHistogram->Clone());
2434 TH1::AddDirectory(bStatus);
2435
2436 return tmp;
2437}
2438
2439////////////////////////////////////////////////////////////////////////////////
2440///returns the dimension of the current TEfficiency object
2441
2443{
2444 return fTotalHistogram->GetDimension();
2445}
2446
2447////////////////////////////////////////////////////////////////////////////////
2448/// Returns the efficiency in the given global bin
2449///
2450/// Note:
2451/// - The estimated efficiency depends on the chosen statistic option:
2452/// for frequentist ones:
2453/// \f$ \hat{\varepsilon} = \frac{passed}{total} \f$
2454/// for bayesian ones the expectation value of the resulting posterior
2455/// distribution is returned:
2456/// \f$ \hat{\varepsilon} = \frac{passed + \alpha}{total + \alpha + \beta} \f$
2457/// If the bit kPosteriorMode is set (or the method TEfficiency::UsePosteriorMode() has been called ) the
2458/// mode (most probable value) of the posterior is returned:
2459/// \f$ \hat{\varepsilon} = \frac{passed + \alpha -1}{total + \alpha + \beta -2} \f$
2460/// - If the denominator is equal to 0, an efficiency of 0 is returned.
2461/// - When \f$ passed + \alpha < 1 \f$ or \f$ total - passed + \beta < 1 \f$ the above
2462/// formula for the mode is not valid. In these cases values the estimated efficiency is 0 or 1.
2463
2465{
2468
2469 if(TestBit(kIsBayesian)) {
2470
2471 // parameters for the beta prior distribution
2474
2475 Double_t aa,bb;
2476 if(TestBit(kUseWeights))
2477 {
2479 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2481
2482 if (tw2 <= 0 ) return pw/tw;
2483
2484 // tw/tw2 renormalize the weights
2485 double norm = tw/tw2;
2486 aa = pw * norm + alpha;
2487 bb = (tw - pw) * norm + beta;
2488 }
2489 else
2490 {
2491 aa = passed + alpha;
2492 bb = total - passed + beta;
2493 }
2494
2495 if (!TestBit(kPosteriorMode) )
2496 return BetaMean(aa,bb);
2497 else
2498 return BetaMode(aa,bb);
2499
2500 }
2501 else
2502 return (total)? ((Double_t)passed)/total : 0;
2503}
2504
2505////////////////////////////////////////////////////////////////////////////////
2506/// Returns the lower error on the efficiency in the given global bin
2507///
2508/// The result depends on the current confidence level fConfLevel and the
2509/// chosen statistic option fStatisticOption. See SetStatisticOption(Int_t) for
2510/// more details.
2511///
2512/// Note: If the histograms are filled with weights, only bayesian methods and the
2513/// normal approximation are supported.
2514
2516{
2519
2520 Double_t eff = GetEfficiency(bin);
2521
2522 // check whether weights have been used
2523 if(TestBit(kUseWeights))
2524 {
2526 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2528 Double_t pw2 = fPassedHistogram->GetSumw2()->At(bin);
2529
2530 if(TestBit(kIsBayesian))
2531 {
2534
2535 if (tw2 <= 0) return 0;
2536
2537 // tw/tw2 renormalize the weights
2538 Double_t norm = tw/tw2;
2539 Double_t aa = pw * norm + alpha;
2540 Double_t bb = (tw - pw) * norm + beta;
2541 Double_t low = 0;
2542 Double_t upper = 1;
2545 }
2546 else {
2548 }
2549
2550 return eff - low;
2551 }
2552 else
2553 {
2555 {
2556 Warning("GetEfficiencyErrorLow","frequentist confidence intervals for weights are only supported by the normal approximation");
2557 Info("GetEfficiencyErrorLow","setting statistic option to kFNormal");
2558 const_cast<TEfficiency*>(this)->SetStatisticOption(kFNormal);
2559 }
2560
2561 Double_t variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
2562 Double_t sigma = sqrt(variance);
2563
2564 Double_t prob = 0.5 * (1.- fConfLevel);
2566
2567 // avoid to return errors which makes eff-err < 0
2568 return (eff - delta < 0) ? eff : delta;
2569 }
2570 }
2571 else
2572 {
2573 if(TestBit(kIsBayesian))
2574 {
2575 // parameters for the beta prior distribution
2578 return (eff - Bayesian(total,passed,fConfLevel,alpha,beta,false,TestBit(kShortestInterval)));
2579 }
2580 else
2581 return (eff - fBoundary(total,passed,fConfLevel,false));
2582 }
2583}
2584
2585////////////////////////////////////////////////////////////////////////////////
2586/// Returns the upper error on the efficiency in the given global bin
2587///
2588/// The result depends on the current confidence level fConfLevel and the
2589/// chosen statistic option fStatisticOption. See SetStatisticOption(Int_t) for
2590/// more details.
2591///
2592/// Note: If the histograms are filled with weights, only bayesian methods and the
2593/// normal approximation are supported.
2594
2596{
2599
2600 Double_t eff = GetEfficiency(bin);
2601
2602 // check whether weights have been used
2603 if(TestBit(kUseWeights))
2604 {
2606 Double_t tw2 = fTotalHistogram->GetSumw2()->At(bin);
2608 Double_t pw2 = fPassedHistogram->GetSumw2()->At(bin);
2609
2610 if(TestBit(kIsBayesian))
2611 {
2614
2615 if (tw2 <= 0) return 0;
2616
2617 // tw/tw2 renormalize the weights
2618 Double_t norm = tw/tw2;
2619 Double_t aa = pw * norm + alpha;
2620 Double_t bb = (tw - pw) * norm + beta;
2621 Double_t low = 0;
2622 Double_t upper = 1;
2625 }
2626 else {
2627 upper = TEfficiency::BetaCentralInterval(fConfLevel,aa,bb,true);
2628 }
2629
2630 return upper - eff;
2631 }
2632 else
2633 {
2635 {
2636 Warning("GetEfficiencyErrorUp","frequentist confidence intervals for weights are only supported by the normal approximation");
2637 Info("GetEfficiencyErrorUp","setting statistic option to kFNormal");
2638 const_cast<TEfficiency*>(this)->SetStatisticOption(kFNormal);
2639 }
2640
2641 Double_t variance = ( pw2 * (1. - 2 * eff) + tw2 * eff *eff ) / ( tw * tw) ;
2642 Double_t sigma = sqrt(variance);
2643
2644 Double_t prob = 0.5 * (1.- fConfLevel);
2646
2647 return (eff + delta > 1) ? 1.-eff : delta;
2648 }
2649 }
2650 else
2651 {
2652 if(TestBit(kIsBayesian))
2653 {
2654 // parameters for the beta prior distribution
2657 return (Bayesian(total,passed,fConfLevel,alpha,beta,true,TestBit(kShortestInterval)) - eff);
2658 }
2659 else
2660 return fBoundary(total,passed,fConfLevel,true) - eff;
2661 }
2662}
2663
2664////////////////////////////////////////////////////////////////////////////////
2665/// Returns the global bin number which can be used as argument for the
2666/// following functions:
2667///
2668/// - GetEfficiency(bin), GetEfficiencyErrorLow(bin), GetEfficiencyErrorUp(bin)
2669/// - SetPassedEvents(bin), SetTotalEvents(bin)
2670///
2671/// see TH1::GetBin() for conventions on numbering bins
2672
2674{
2675 return fTotalHistogram->GetBin(binx,biny,binz);
2676}
2677
2678////////////////////////////////////////////////////////////////////////////////
2679
2681{
2682 return (fFunctions) ? fFunctions : fFunctions = new TList();
2683}
2684
2685////////////////////////////////////////////////////////////////////////////////
2686/// Merges the TEfficiency objects in the given list to the given
2687/// TEfficiency object using the operator+=(TEfficiency&)
2688///
2689/// The merged result is stored in the current object. The statistic options and
2690/// the confidence level are taken from the current object.
2691///
2692/// This function should be used when all TEfficiency objects correspond to
2693/// the same process.
2694///
2695/// The new weight is set according to:
2696/// \f$ \frac{1}{w_{new}} = \sum_{i} \frac{1}{w_{i}} \f$
2697
2699{
2700 if(!pList->IsEmpty()) {
2701 TIter next(pList);
2702 TObject* obj = 0;
2703 TEfficiency* pEff = 0;
2704 while((obj = next())) {
2705 pEff = dynamic_cast<TEfficiency*>(obj);
2706 if(pEff) {
2707 *this += *pEff;
2708 }
2709 }
2710 }
2712}
2713
2714////////////////////////////////////////////////////////////////////////////////
2715/**
2716Returns the confidence limits for the efficiency supposing that the
2717efficiency follows a normal distribution with the rms below
2718
2719\param[in] total number of total events
2720\param[in] passed 0 <= number of passed events <= total
2721\param[in] level confidence level
2722\param[in] bUpper
2723 - true - upper boundary is returned
2724 - false - lower boundary is returned
2725
2726Calculation:
2727
2728\f{eqnarray*}{
2729 \hat{\varepsilon} &=& \frac{passed}{total}\\
2730 \sigma_{\varepsilon} &=& \sqrt{\frac{\hat{\varepsilon} (1 - \hat{\varepsilon})}{total}}\\
2731 \varepsilon_{low} &=& \hat{\varepsilon} \pm \Phi^{-1}(\frac{level}{2},\sigma_{\varepsilon})
2732\f}
2733*/
2734
2736{
2737 Double_t alpha = (1.0 - level)/2;
2738 if (total == 0) return (bUpper) ? 1 : 0;
2739 Double_t average = passed / total;
2740 Double_t sigma = std::sqrt(average * (1 - average) / total);
2741 Double_t delta = ROOT::Math::normal_quantile(1 - alpha,sigma);
2742
2743 if(bUpper)
2744 return ((average + delta) > 1) ? 1.0 : (average + delta);
2745 else
2746 return ((average - delta) < 0) ? 0.0 : (average - delta);
2747}
2748
2749////////////////////////////////////////////////////////////////////////////////
2750/// Adds the histograms of another TEfficiency object to current histograms
2751///
2752/// The statistic options and the confidence level remain unchanged.
2753///
2754/// fTotalHistogram += rhs.fTotalHistogram;
2755/// fPassedHistogram += rhs.fPassedHistogram;
2756///
2757/// calculates a new weight:
2758/// current weight of this TEfficiency object = \f$ w_{1} \f$
2759/// weight of rhs = \f$ w_{2} \f$
2760/// \f$ w_{new} = \frac{w_{1} \times w_{2}}{w_{1} + w_{2}} \f$
2761
2763{
2764
2765 if (fTotalHistogram == 0 && fPassedHistogram == 0) {
2766 // efficiency is empty just copy it over
2767 *this = rhs;
2768 return *this;
2769 }
2770 else if (fTotalHistogram == 0 || fPassedHistogram == 0) {
2771 Fatal("operator+=","Adding to a non consistent TEfficiency object which has not a total or a passed histogram ");
2772 return *this;
2773 }
2774
2775 if (rhs.fTotalHistogram == 0 && rhs.fPassedHistogram == 0 ) {
2776 Warning("operator+=","no operation: adding an empty object");
2777 return *this;
2778 }
2779 else if (rhs.fTotalHistogram == 0 || rhs.fPassedHistogram == 0 ) {
2780 Fatal("operator+=","Adding a non consistent TEfficiency object which has not a total or a passed histogram ");
2781 return *this;
2782 }
2783
2786
2789
2790 SetWeight((fWeight * rhs.GetWeight())/(fWeight + rhs.GetWeight()));
2791
2792 return *this;
2793}
2794
2795////////////////////////////////////////////////////////////////////////////////
2796/// Assignment operator
2797///
2798/// The histograms, statistic option, confidence level, weight and paint styles
2799/// of rhs are copied to the this TEfficiency object.
2800///
2801/// Note: - The list of associated functions is not copied. After this
2802/// operation the list of associated functions is empty.
2803
2805{
2806 if(this != &rhs)
2807 {
2808 //statistic options
2812 SetBetaBeta(rhs.GetBetaBeta());
2813 SetWeight(rhs.GetWeight());
2814
2815 //associated list of functions
2816 if(fFunctions)
2817 fFunctions->Delete();
2818
2819 //copy histograms
2820 delete fTotalHistogram;
2821 delete fPassedHistogram;
2822
2823 Bool_t bStatus = TH1::AddDirectoryStatus();
2827 TH1::AddDirectory(bStatus);
2828
2829 //delete temporary paint objects
2830 delete fPaintHisto;
2831 delete fPaintGraph;
2832 fPaintHisto = 0;
2833 fPaintGraph = 0;
2834
2835 //copy style
2836 rhs.TAttLine::Copy(*this);
2837 rhs.TAttFill::Copy(*this);
2838 rhs.TAttMarker::Copy(*this);
2839 }
2840
2841 return *this;
2842}
2843
2844////////////////////////////////////////////////////////////////////////////////
2845/// Paints this TEfficiency object
2846///
2847/// For details on the possible option see Draw(Option_t*)
2848///
2849/// Note for 1D classes
2850/// In 1D the TEfficiency uses a TGraphAsymmErrors for drawing
2851/// The TGraph is created only the first time Paint is used. The user can manipulate the
2852/// TGraph via the method TEfficiency::GetPaintedGraph()
2853/// The TGraph creates behing an histogram for the axis. The histogram is created also only the first time.
2854/// If the axis needs to be updated because in the meantime the class changed use this trick
2855/// which will trigger a re-calculation of the axis of the graph
2856/// TEfficiency::GetPaintedGraph()->Set(0)
2857///
2858/// Note that in order to access the painted graph via GetPaintedGraph() you need either to call Paint or better
2859/// gPad->Update();
2860///
2861
2863{
2864
2865
2866 if(!gPad)
2867 return;
2868
2869
2870 //use TGraphAsymmErrors for painting
2871 if(GetDimension() == 1) {
2872 if(!fPaintGraph) {
2873 fPaintGraph = CreateGraph(opt);
2874 }
2875 else
2876 // update existing graph already created
2877 FillGraph(fPaintGraph, opt);
2878
2879 //paint graph
2880
2881 fPaintGraph->Paint(opt);
2882
2883 //paint all associated functions
2884 if(fFunctions) {
2885 //paint box with fit parameters
2886 //the fit statistics will be painted if gStyle->SetOptFit(1) has been
2887 // called by the user
2888 TIter next(fFunctions);
2889 TObject* obj = 0;
2890 while((obj = next())) {
2891 if(obj->InheritsFrom(TF1::Class())) {
2892 fPaintGraph->PaintStats((TF1*)obj);
2893 ((TF1*)obj)->Paint("sameC");
2894 }
2895 }
2896 }
2897
2898 return;
2899 }
2900
2901 //use TH2 for painting
2902 if(GetDimension() == 2) {
2903 if(!fPaintHisto) {
2905 }
2906 else
2908
2909 //paint histogram
2910 fPaintHisto->Paint(opt);
2911 return;
2912 }
2913 Warning("Paint","Painting 3D efficiency is not implemented");
2914}
2915
2916////////////////////////////////////////////////////////////////////////////////
2917/// Have histograms fixed bins along each axis?
2918
2919void TEfficiency::SavePrimitive(std::ostream& out,Option_t* opt)
2920{
2921 Bool_t equi_bins = true;
2922
2923 //indentation
2924 TString indent = " ";
2925 //names for arrays containing the bin edges
2926 //static counter needed if more objects are saved
2927 static Int_t naxis = 0;
2928 TString sxaxis="xAxis",syaxis="yAxis",szaxis="zAxis";
2929
2930 //note the missing break statements!
2931 switch(GetDimension()) {
2932 case 3:
2933 equi_bins = equi_bins && !fTotalHistogram->GetZaxis()->GetXbins()->fArray
2935 case 2:
2936 equi_bins = equi_bins && !fTotalHistogram->GetYaxis()->GetXbins()->fArray
2938 case 1:
2939 equi_bins = equi_bins && !fTotalHistogram->GetXaxis()->GetXbins()->fArray
2941 }
2942
2943 //create arrays containing the variable binning
2944 if(!equi_bins) {
2945 Int_t i;
2946 ++naxis;
2947 sxaxis += naxis;
2948 syaxis += naxis;
2949 szaxis += naxis;
2950 //x axis
2951 out << indent << "Double_t " << sxaxis << "["
2952 << fTotalHistogram->GetXaxis()->GetXbins()->fN << "] = {";
2953 for (i = 0; i < fTotalHistogram->GetXaxis()->GetXbins()->fN; ++i) {
2954 if (i != 0) out << ", ";
2955 out << fTotalHistogram->GetXaxis()->GetXbins()->fArray[i];
2956 }
2957 out << "}; " << std::endl;
2958 //y axis
2959 if(GetDimension() > 1) {
2960 out << indent << "Double_t " << syaxis << "["
2961 << fTotalHistogram->GetYaxis()->GetXbins()->fN << "] = {";
2962 for (i = 0; i < fTotalHistogram->GetYaxis()->GetXbins()->fN; ++i) {
2963 if (i != 0) out << ", ";
2964 out << fTotalHistogram->GetYaxis()->GetXbins()->fArray[i];
2965 }
2966 out << "}; " << std::endl;
2967 }
2968 //z axis
2969 if(GetDimension() > 2) {
2970 out << indent << "Double_t " << szaxis << "["
2971 << fTotalHistogram->GetZaxis()->GetXbins()->fN << "] = {";
2972 for (i = 0; i < fTotalHistogram->GetZaxis()->GetXbins()->fN; ++i) {
2973 if (i != 0) out << ", ";
2974 out << fTotalHistogram->GetZaxis()->GetXbins()->fArray[i];
2975 }
2976 out << "}; " << std::endl;
2977 }
2978 }//creating variable binning
2979
2980 //TEfficiency pointer has efficiency name + counter
2981 static Int_t eff_count = 0;
2982 ++eff_count;
2983 TString eff_name = GetName();
2984 eff_name += eff_count;
2985
2986 const char* name = eff_name.Data();
2987
2988 //construct TEfficiency object
2989 const char quote = '"';
2990 out << indent << std::endl;
2991 out << indent << ClassName() << " * " << name << " = new " << ClassName()
2992 << "(" << quote << GetName() << quote << "," << quote
2993 << GetTitle() << quote <<",";
2994 //fixed bin size -> use n,min,max constructor
2995 if(equi_bins) {
2996 out << fTotalHistogram->GetXaxis()->GetNbins() << ","
2997 << fTotalHistogram->GetXaxis()->GetXmin() << ","
2999 if(GetDimension() > 1) {
3000 out << "," << fTotalHistogram->GetYaxis()->GetNbins() << ","
3001 << fTotalHistogram->GetYaxis()->GetXmin() << ","
3003 }
3004 if(GetDimension() > 2) {
3005 out << "," << fTotalHistogram->GetZaxis()->GetNbins() << ","
3006 << fTotalHistogram->GetZaxis()->GetXmin() << ","
3008 }
3009 }
3010 //variable bin size -> use n,*bins constructor
3011 else {
3012 out << fTotalHistogram->GetXaxis()->GetNbins() << "," << sxaxis;
3013 if(GetDimension() > 1)
3014 out << "," << fTotalHistogram->GetYaxis()->GetNbins() << ","
3015 << syaxis;
3016 if(GetDimension() > 2)
3017 out << "," << fTotalHistogram->GetZaxis()->GetNbins() << ","
3018 << szaxis;
3019 }
3020 out << ");" << std::endl;
3021 out << indent << std::endl;
3022
3023 //set statistic options
3024 out << indent << name << "->SetConfidenceLevel(" << fConfLevel << ");"
3025 << std::endl;
3026 out << indent << name << "->SetBetaAlpha(" << fBeta_alpha << ");"
3027 << std::endl;
3028 out << indent << name << "->SetBetaBeta(" << fBeta_beta << ");" << std::endl;
3029 out << indent << name << "->SetWeight(" << fWeight << ");" << std::endl;
3030 out << indent << name << "->SetStatisticOption(" << fStatisticOption << ");"
3031 << std::endl;
3032 out << indent << name << "->SetPosteriorMode(" << TestBit(kPosteriorMode) << ");" << std::endl;
3033 out << indent << name << "->SetShortestInterval(" << TestBit(kShortestInterval) << ");" << std::endl;
3034 if(TestBit(kUseWeights))
3035 out << indent << name << "->SetUseWeightedEvents();" << std::endl;
3036
3037 // save bin-by-bin prior parameters
3038 for(unsigned int i = 0; i < fBeta_bin_params.size(); ++i)
3039 {
3040 out << indent << name << "->SetBetaBinParameters(" << i << "," << fBeta_bin_params.at(i).first
3041 << "," << fBeta_bin_params.at(i).second << ");" << std::endl;
3042 }
3043
3044 //set bin contents
3045 Int_t nbins = fTotalHistogram->GetNbinsX() + 2;
3046 if(GetDimension() > 1)
3047 nbins *= fTotalHistogram->GetNbinsY() + 2;
3048 if(GetDimension() > 2)
3049 nbins *= fTotalHistogram->GetNbinsZ() + 2;
3050
3051 //important: set first total number than passed number
3052 for(Int_t i = 0; i < nbins; ++i) {
3053 out << indent << name <<"->SetTotalEvents(" << i << "," <<
3054 fTotalHistogram->GetBinContent(i) << ");" << std::endl;
3055 out << indent << name <<"->SetPassedEvents(" << i << "," <<
3056 fPassedHistogram->GetBinContent(i) << ");" << std::endl;
3057 }
3058
3059 //save list of functions
3060 TIter next(fFunctions);
3061 TObject* obj = 0;
3062 while((obj = next())) {
3063 obj->SavePrimitive(out,"nodraw");
3064 if(obj->InheritsFrom(TF1::Class())) {
3065 out << indent << name << "->GetListOfFunctions()->Add("
3066 << obj->GetName() << ");" << std::endl;
3067 }
3068 }
3069
3070 //set style
3074
3075 //draw TEfficiency object
3076 TString option = opt;
3077 option.ToLower();
3078 if (!option.Contains("nodraw"))
3079 out<< indent << name<< "->Draw(" << quote << opt << quote << ");"
3080 << std::endl;
3081}
3082
3083////////////////////////////////////////////////////////////////////////////////
3084/// Sets the shape parameter &alpha;
3085///
3086/// The prior probability of the efficiency is given by the beta distribution:
3087/// \f[
3088/// f(\varepsilon;\alpha;\beta) = \frac{1}{B(\alpha,\beta)} \varepsilon^{\alpha-1} (1 - \varepsilon)^{\beta-1}
3089/// \f]
3090///
3091/// Note: - both shape parameters have to be positive (i.e. > 0)
3092
3094{
3095 if(alpha > 0)
3096 fBeta_alpha = alpha;
3097 else
3098 Warning("SetBetaAlpha(Double_t)","invalid shape parameter %.2lf",alpha);
3099}
3100
3101////////////////////////////////////////////////////////////////////////////////
3102/// Sets the shape parameter &beta;
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(beta > 0)
3114 fBeta_beta = beta;
3115 else
3116 Warning("SetBetaBeta(Double_t)","invalid shape parameter %.2lf",beta);
3117}
3118
3119////////////////////////////////////////////////////////////////////////////////
3120/// Sets different shape parameter &alpha; and &beta;
3121/// for the prior distribution for each bin. By default the global parameter are used if they are not set
3122/// for the specific bin
3123/// The prior probability of the efficiency is given by the beta distribution:
3124/// \f[
3125/// f(\varepsilon;\alpha;\beta) = \frac{1}{B(\alpha,\beta)} \varepsilon^{\alpha-1} (1 - \varepsilon)^{\beta-1}
3126/// \f]
3127///
3128/// Note:
3129/// - both shape parameters have to be positive (i.e. > 0)
3130/// - bin gives the global bin number (cf. GetGlobalBin)
3131
3133{
3134 if (!fPassedHistogram || !fTotalHistogram) return;
3136 // doing this I get h1->fN which is available only for a TH1D
3137 UInt_t n = h1->GetBin(h1->GetNbinsX()+1, h1->GetNbinsY()+1, h1->GetNbinsZ()+1 ) + 1;
3138
3139 // in case vector is not created do with default alpha, beta params
3140 if (fBeta_bin_params.size() != n )
3141 fBeta_bin_params = std::vector<std::pair<Double_t, Double_t> >(n, std::make_pair(fBeta_alpha, fBeta_beta) );
3142
3143 // vector contains also values for under/overflows
3144 fBeta_bin_params[bin] = std::make_pair(alpha,beta);
3145 SetBit(kUseBinPrior,true);
3146
3147}
3148
3149////////////////////////////////////////////////////////////////////////////////
3150/// Set the bins for the underlined passed and total histograms
3151/// If the class have been already filled the previous contents will be lost
3152
3154{
3155 if (GetDimension() != 1) {
3156 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3157 return kFALSE;
3158 }
3159 if (fTotalHistogram->GetEntries() != 0 ) {
3160 Warning("SetBins","Histogram entries will be lost after SetBins");
3163 }
3166 return kTRUE;
3167}
3168
3169////////////////////////////////////////////////////////////////////////////////
3170/// Set the bins for the underlined passed and total histograms
3171/// If the class have been already filled the previous contents will be lost
3172
3174{
3175 if (GetDimension() != 1) {
3176 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3177 return kFALSE;
3178 }
3179 if (fTotalHistogram->GetEntries() != 0 ) {
3180 Warning("SetBins","Histogram entries will be lost after SetBins");
3183 }
3184 fPassedHistogram->SetBins(nx,xBins);
3185 fTotalHistogram->SetBins(nx,xBins);
3186 return kTRUE;
3187}
3188
3189////////////////////////////////////////////////////////////////////////////////
3190/// Set the bins for the underlined passed and total histograms
3191/// If the class have been already filled the previous contents will be lost
3192
3194{
3195 if (GetDimension() != 2) {
3196 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3197 return kFALSE;
3198 }
3199 if (fTotalHistogram->GetEntries() != 0 ) {
3200 Warning("SetBins","Histogram entries will be lost after SetBins");
3203 }
3206 return kTRUE;
3207}
3208
3209////////////////////////////////////////////////////////////////////////////////
3210/// Set the bins for the underlined passed and total histograms
3211/// If the class have been already filled the previous contents will be lost
3212
3213Bool_t TEfficiency::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins)
3214{
3215 if (GetDimension() != 2) {
3216 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3217 return kFALSE;
3218 }
3219 if (fTotalHistogram->GetEntries() != 0 ) {
3220 Warning("SetBins","Histogram entries will be lost after SetBins");
3223 }
3224 fPassedHistogram->SetBins(nx,xBins,ny,yBins);
3225 fTotalHistogram->SetBins(nx,xBins,ny,yBins);
3226 return kTRUE;
3227}
3228
3229////////////////////////////////////////////////////////////////////////////////
3230/// Set the bins for the underlined passed and total histograms
3231/// If the class have been already filled the previous contents will be lost
3232
3234 Int_t nz, Double_t zmin, Double_t zmax)
3235{
3236 if (GetDimension() != 3) {
3237 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3238 return kFALSE;
3239 }
3240 if (fTotalHistogram->GetEntries() != 0 ) {
3241 Warning("SetBins","Histogram entries will be lost after SetBins");
3244 }
3245 fPassedHistogram->SetBins(nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
3246 fTotalHistogram->SetBins (nx,xmin,xmax,ny,ymin,ymax,nz,zmin,zmax);
3247 return kTRUE;
3248}
3249
3250////////////////////////////////////////////////////////////////////////////////
3251/// Set the bins for the underlined passed and total histograms
3252/// If the class have been already filled the previous contents will be lost
3253
3254Bool_t TEfficiency::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins, Int_t nz,
3255 const Double_t *zBins )
3256{
3257 if (GetDimension() != 3) {
3258 Error("SetBins","Using wrong SetBins function for a %d-d histogram",GetDimension());
3259 return kFALSE;
3260 }
3261 if (fTotalHistogram->GetEntries() != 0 ) {
3262 Warning("SetBins","Histogram entries will be lost after SetBins");
3265 }
3266 fPassedHistogram->SetBins(nx,xBins,ny,yBins,nz,zBins);
3267 fTotalHistogram->SetBins(nx,xBins,ny,yBins,nz,zBins);
3268 return kTRUE;
3269}
3270
3271////////////////////////////////////////////////////////////////////////////////
3272/// Sets the confidence level (0 < level < 1)
3273/// The default value is 1-sigma :~ 0.683
3274
3276{
3277 if((level > 0) && (level < 1))
3278 fConfLevel = level;
3279 else
3280 Warning("SetConfidenceLevel(Double_t)","invalid confidence level %.2lf",level);
3281}
3282
3283////////////////////////////////////////////////////////////////////////////////
3284/// Sets the directory holding this TEfficiency object
3285///
3286/// A reference to this TEfficiency object is removed from the current
3287/// directory (if it exists) and a new reference to this TEfficiency object is
3288/// added to the given directory.
3289///
3290/// Notes: - If the given directory is 0, the TEfficiency object does not
3291/// belong to any directory and will not be written to file during the
3292/// next TFile::Write() command.
3293
3295{
3296 if(fDirectory == dir)
3297 return;
3298 if(fDirectory)
3299 fDirectory->Remove(this);
3300 fDirectory = dir;
3301 if(fDirectory)
3302 fDirectory->Append(this);
3303}
3304
3305////////////////////////////////////////////////////////////////////////////////
3306/// Sets the name
3307///
3308/// Note: The names of the internal histograms are set to "name + _total" and
3309/// "name + _passed" respectively.
3310
3312{
3314
3315 //setting the names (appending the correct ending)
3316 TString name_total = name + TString("_total");
3317 TString name_passed = name + TString("_passed");
3318 fTotalHistogram->SetName(name_total);
3319 fPassedHistogram->SetName(name_passed);
3320}
3321
3322////////////////////////////////////////////////////////////////////////////////
3323/// Sets the number of passed events in the given global bin
3324///
3325/// returns "true" if the number of passed events has been updated
3326/// otherwise "false" ist returned
3327///
3328/// Note: - requires: 0 <= events <= fTotalHistogram->GetBinContent(bin)
3329
3331{
3332 if(events <= fTotalHistogram->GetBinContent(bin)) {
3333 fPassedHistogram->SetBinContent(bin,events);
3334 return true;
3335 }
3336 else {
3337 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);
3338 return false;
3339 }
3340}
3341
3342////////////////////////////////////////////////////////////////////////////////
3343/// Sets the histogram containing the passed events
3344///
3345/// The given histogram is cloned and stored internally as histogram containing
3346/// the passed events. The given histogram has to be consistent with the current
3347/// fTotalHistogram (see CheckConsistency(const TH1&,const TH1&)).
3348/// The method returns whether the fPassedHistogram has been replaced (true) or
3349/// not (false).
3350///
3351/// Note: The list of associated functions fFunctions is cleared.
3352///
3353/// Option:
3354/// - "f": force the replacement without checking the consistency
3355/// This can lead to inconsistent histograms and useless results
3356/// or unexpected behaviour. But sometimes it might be the only
3357/// way to change the histograms. If you use this option, you
3358/// should ensure that the fTotalHistogram is replaced by a
3359/// consistent one (with respect to rPassed) as well.
3360
3362{
3363 TString option = opt;
3364 option.ToLower();
3365
3366 Bool_t bReplace = option.Contains("f");
3367
3368 if(!bReplace)
3369 bReplace = CheckConsistency(rPassed,*fTotalHistogram);
3370
3371 if(bReplace) {
3372 delete fPassedHistogram;
3373 Bool_t bStatus = TH1::AddDirectoryStatus();
3375 fPassedHistogram = (TH1*)(rPassed.Clone());
3377 TH1::AddDirectory(bStatus);
3378
3379 if(fFunctions)
3380 fFunctions->Delete();
3381
3382 //check whether both histograms are filled with weights
3383 bool useWeights = CheckWeights(rPassed,*fTotalHistogram);
3384
3385 SetUseWeightedEvents(useWeights);
3386
3387 return true;
3388 }
3389 else
3390 return false;
3391}
3392
3393////////////////////////////////////////////////////////////////////////////////
3394/// Sets the statistic option which affects the calculation of the confidence interval
3395///
3396/// Options:
3397/// - kFCP (=0)(default): using the Clopper-Pearson interval (recommended by PDG)
3398/// sets kIsBayesian = false
3399/// see also ClopperPearson
3400/// - kFNormal (=1) : using the normal approximation
3401/// sets kIsBayesian = false
3402/// see also Normal
3403/// - kFWilson (=2) : using the Wilson interval
3404/// sets kIsBayesian = false
3405/// see also Wilson
3406/// - kFAC (=3) : using the Agresti-Coull interval
3407/// sets kIsBayesian = false
3408/// see also AgrestiCoull
3409/// - kFFC (=4) : using the Feldman-Cousins frequentist method
3410/// sets kIsBayesian = false
3411/// see also FeldmanCousins
3412/// - kBJeffrey (=5) : using the Jeffrey interval
3413/// sets kIsBayesian = true, fBeta_alpha = 0.5 and fBeta_beta = 0.5
3414/// see also Bayesian
3415/// - kBUniform (=6) : using a uniform prior
3416/// sets kIsBayesian = true, fBeta_alpha = 1 and fBeta_beta = 1
3417/// see also Bayesian
3418/// - kBBayesian (=7) : using a custom prior defined by fBeta_alpha and fBeta_beta
3419/// sets kIsBayesian = true
3420/// see also Bayesian
3421/// - kMidP (=8) : using the Lancaster Mid-P method
3422/// sets kIsBayesian = false
3423
3424
3426{
3427 fStatisticOption = option;
3428
3429 switch(option)
3430 {
3431 case kFCP:
3433 SetBit(kIsBayesian,false);
3434 break;
3435 case kFNormal:
3436 fBoundary = &Normal;
3437 SetBit(kIsBayesian,false);
3438 break;
3439 case kFWilson:
3440 fBoundary = &Wilson;
3441 SetBit(kIsBayesian,false);
3442 break;
3443 case kFAC:
3445 SetBit(kIsBayesian,false);
3446 break;
3447 case kFFC:
3449 SetBit(kIsBayesian,false);
3450 break;
3451 case kMidP:
3453 SetBit(kIsBayesian,false);
3454 break;
3455 case kBJeffrey:
3456 fBeta_alpha = 0.5;
3457 fBeta_beta = 0.5;
3458 SetBit(kIsBayesian,true);
3459 SetBit(kUseBinPrior,false);
3460 break;
3461 case kBUniform:
3462 fBeta_alpha = 1;
3463 fBeta_beta = 1;
3464 SetBit(kIsBayesian,true);
3465 SetBit(kUseBinPrior,false);
3466 break;
3467 case kBBayesian:
3468 SetBit(kIsBayesian,true);
3469 break;
3470 default:
3473 SetBit(kIsBayesian,false);
3474 }
3475}
3476
3477////////////////////////////////////////////////////////////////////////////////
3478/// Sets the title
3479///
3480/// Notes:
3481/// - The titles of the internal histograms are set to "title + (total)"
3482/// or "title + (passed)" respectively.
3483/// - It is possible to label the axis of the histograms as usual (see
3484/// TH1::SetTitle).
3485///
3486/// Example: Setting the title to "My Efficiency" and label the axis
3487/// pEff->SetTitle("My Efficiency;x label;eff");
3488
3489void TEfficiency::SetTitle(const char* title)
3490{
3491
3492 //setting the titles (looking for the first semicolon and insert the tokens there)
3493 TString title_passed = title;
3494 TString title_total = title;
3495 Ssiz_t pos = title_passed.First(";");
3496 if (pos != kNPOS) {
3497 title_passed.Insert(pos," (passed)");
3498 title_total.Insert(pos," (total)");
3499 }
3500 else {
3501 title_passed.Append(" (passed)");
3502 title_total.Append(" (total)");
3503 }
3504 fPassedHistogram->SetTitle(title_passed);
3505 fTotalHistogram->SetTitle(title_total);
3506
3507 // strip (total) for the TEfficiency title
3508 // HIstogram SetTitle has already stripped the axis
3509 TString teffTitle = fTotalHistogram->GetTitle();
3510 teffTitle.ReplaceAll(" (total)","");
3511 TNamed::SetTitle(teffTitle);
3512
3513}
3514
3515////////////////////////////////////////////////////////////////////////////////
3516/// Sets the number of total events in the given global bin
3517///
3518/// returns "true" if the number of total events has been updated
3519/// otherwise "false" ist returned
3520///
3521/// Note: - requires: fPassedHistogram->GetBinContent(bin) <= events
3522
3524{
3525 if(events >= fPassedHistogram->GetBinContent(bin)) {
3526 fTotalHistogram->SetBinContent(bin,events);
3527 return true;
3528 }
3529 else {
3530 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);
3531 return false;
3532 }
3533}
3534
3535////////////////////////////////////////////////////////////////////////////////
3536/// Sets the histogram containing all events
3537///
3538/// The given histogram is cloned and stored internally as histogram containing
3539/// all events. The given histogram has to be consistent with the current
3540/// fPassedHistogram (see CheckConsistency(const TH1&,const TH1&)).
3541/// The method returns whether the fTotalHistogram has been replaced (true) or
3542/// not (false).
3543///
3544/// Note: The list of associated functions fFunctions is cleared.
3545///
3546/// Option:
3547/// - "f": force the replacement without checking the consistency
3548/// This can lead to inconsistent histograms and useless results
3549/// or unexpected behaviour. But sometimes it might be the only
3550/// way to change the histograms. If you use this option, you
3551/// should ensure that the fPassedHistogram is replaced by a
3552/// consistent one (with respect to rTotal) as well.
3553
3555{
3556 TString option = opt;
3557 option.ToLower();
3558
3559 Bool_t bReplace = option.Contains("f");
3560
3561 if(!bReplace)
3562 bReplace = CheckConsistency(*fPassedHistogram,rTotal);
3563
3564 if(bReplace) {
3565 delete fTotalHistogram;
3566 Bool_t bStatus = TH1::AddDirectoryStatus();
3568 fTotalHistogram = (TH1*)(rTotal.Clone());
3570 TH1::AddDirectory(bStatus);
3571
3572 if(fFunctions)
3573 fFunctions->Delete();
3574
3575 //check whether both histograms are filled with weights
3576 bool useWeights = CheckWeights(*fPassedHistogram,rTotal);
3577 SetUseWeightedEvents(useWeights);
3578
3579 return true;
3580 }
3581 else
3582 return false;
3583}
3584
3585////////////////////////////////////////////////////////////////////////////////
3586
3588{
3589 if (on && !TestBit(kUseWeights) )
3590 gROOT->Info("TEfficiency::SetUseWeightedEvents","Handle weighted events for computing efficiency");
3591
3592 SetBit(kUseWeights,on);
3593
3598}
3599
3600////////////////////////////////////////////////////////////////////////////////
3601/// Sets the global weight for this TEfficiency object
3602///
3603/// Note: - weight has to be positive ( > 0)
3604
3606{
3607 if(weight > 0)
3608 fWeight = weight;
3609 else
3610 Warning("SetWeight","invalid weight %.2lf",weight);
3611}
3612
3613////////////////////////////////////////////////////////////////////////////////
3614/**
3615Calculates the boundaries for the frequentist Wilson interval
3616
3617\param[in] total number of total events
3618\param[in] passed 0 <= number of passed events <= total
3619\param[in] level confidence level
3620\param[in] bUpper
3621 - true - upper boundary is returned
3622 - false - lower boundary is returned
3623
3624Calculation:
3625\f{eqnarray*}{
3626 \alpha &=& 1 - \frac{level}{2}\\
3627 \kappa &=& \Phi^{-1}(1 - \alpha,1) ...\ normal\ quantile\ function\\
3628 mode &=& \frac{passed + \frac{\kappa^{2}}{2}}{total + \kappa^{2}}\\
3629 \Delta &=& \frac{\kappa}{total + \kappa^{2}} * \sqrt{passed (1 - \frac{passed}{total}) + \frac{\kappa^{2}}{4}}\\
3630 return &=& max(0,mode - \Delta)\ or\ min(1,mode + \Delta)
3631\f}
3632
3633*/
3634
3636{
3637 Double_t alpha = (1.0 - level)/2;
3638 if (total == 0) return (bUpper) ? 1 : 0;
3639 Double_t average = ((Double_t)passed) / total;
3640 Double_t kappa = ROOT::Math::normal_quantile(1 - alpha,1);
3641
3642 Double_t mode = (passed + 0.5 * kappa * kappa) / (total + kappa * kappa);
3643 Double_t delta = kappa / (total + kappa*kappa) * std::sqrt(total * average
3644 * (1 - average) + kappa * kappa / 4);
3645 if(bUpper)
3646 return ((mode + delta) > 1) ? 1.0 : (mode + delta);
3647 else
3648 return ((mode - delta) < 0) ? 0.0 : (mode - delta);
3649}
3650
3651////////////////////////////////////////////////////////////////////////////////
3652/// Addition operator
3653///
3654/// adds the corresponding histograms:
3655/// ~~~ {.cpp}
3656/// lhs.GetTotalHistogram() + rhs.GetTotalHistogram()
3657/// lhs.GetPassedHistogram() + rhs.GetPassedHistogram()
3658/// ~~~
3659/// the statistic option and the confidence level are taken from lhs
3660
3661const TEfficiency operator+(const TEfficiency& lhs,const TEfficiency& rhs)
3662{
3663 TEfficiency tmp(lhs);
3664 tmp += rhs;
3665 return tmp;
3666}
3667
3668#endif
void Class()
Definition: Class.C:29
SVector< double, 2 > v
Definition: Dict.h:5
#define b(i)
Definition: RSha256.hxx:100
#define e(i)
Definition: RSha256.hxx:103
const Ssiz_t kNPOS
Definition: RtypesCore.h:111
int Int_t
Definition: RtypesCore.h:41
int Ssiz_t
Definition: RtypesCore.h:63
unsigned int UInt_t
Definition: RtypesCore.h:42
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
long long Long64_t
Definition: RtypesCore.h:69
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
static void indent(ostringstream &buf, int indent_level)
#define gDirectory
Definition: TDirectory.h:218
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:109
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
double sqrt(double)
TRObject operator()(const T1 &t1) const
#define gROOT
Definition: TROOT.h:414
#define gPad
Definition: TVirtualPad.h:286
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:202
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:234
Line Attributes class.
Definition: TAttLine.h:18
void Copy(TAttLine &attline) const
Copy this line attributes to a new TAttLine.
Definition: TAttLine.cxx:164
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:262
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:245
void Copy(TAttMarker &attmarker) const
Copy this marker attributes to a new TAttMarker.
Definition: TAttMarker.cxx:210
Class to manage histogram axis.
Definition: TAxis.h:30
Bool_t IsVariableBinSize() const
Definition: TAxis.h:136
const TArrayD * GetXbins() const
Definition: TAxis.h:130
Double_t GetXmax() const
Definition: TAxis.h:134
virtual Double_t GetBinLowEdge(Int_t bin) const
Return low edge of bin.
Definition: TAxis.cxx:504
virtual Int_t FindFixBin(Double_t x) const
Find bin number corresponding to abscissa x.
Definition: TAxis.cxx:405
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
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:63
virtual Bool_t IsEmpty() const
Definition: TCollection.h:186
Describe directory structure in memory.
Definition: TDirectory.h:34
virtual void Append(TObject *obj, Bool_t replace=kFALSE)
Append object to this directory.
Definition: TDirectory.cxx:190
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
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
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
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
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
temporary histogram for painting
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 α.
@ kShortestInterval
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
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)
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
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
pointer to a method calculating the boundaries of confidence intervals
Definition: TEfficiency.h:51
Double_t fBeta_beta
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
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 directory holding this TEfficiency object
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.
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 graph for painting
Definition: TEfficiency.h:55
1-Dim function class
Definition: TF1.h:211
Provides an indirection to the TFitResult class and with a semantics identical to a TFitResult pointe...
Definition: TFitResultPtr.h:31
TGraph with asymmetric error bars.
virtual void Paint(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:1959
virtual void PaintStats(TF1 *fit)
Draw the stats.
Definition: TGraph.cxx:1986
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a graph.
Definition: TGraph.cxx:791
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TGraph.cxx:968
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
The TH1 histogram class.
Definition: TH1.h:56
virtual void SetDirectory(TDirectory *dir)
By default when an histogram is created, it is added to the list of histogram objects in the current ...
Definition: TH1.cxx:8351
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6309
virtual void SetNormFactor(Double_t factor=1)
Definition: TH1.h:400
virtual Double_t GetBinCenter(Int_t bin) const
Return bin center for 1D histogram.
Definition: TH1.cxx:8554
TAxis * GetZaxis()
Definition: TH1.h:318
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:7308
virtual Int_t GetNbinsY() const
Definition: TH1.h:293
virtual Int_t GetNbinsZ() const
Definition: TH1.h:294
virtual Int_t GetDimension() const
Definition: TH1.h:278
static void AddDirectory(Bool_t add=kTRUE)
Sets the flag controlling the automatic add of histograms in memory.
Definition: TH1.cxx:1225
@ kIsAverage
Bin contents are average (used by Add)
Definition: TH1.h:166
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH1.cxx:6700
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
virtual Int_t GetNcells() const
Definition: TH1.h:295
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2664
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:4784
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
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:777
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3258
TAxis * GetYaxis()
Definition: TH1.h:317
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:8635
@ kNstat
Definition: TH1.h:179
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8565
virtual Double_t GetEntries() const
Return the current number of entries.
Definition: TH1.cxx:4277
virtual void SetName(const char *name)
Change the name of this histogram.
Definition: TH1.cxx:8374
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4882
virtual TArrayD * GetSumw2()
Definition: TH1.h:308
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Execute action corresponding to one event.
Definition: TH1.cxx:3154
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8576
virtual void Paint(Option_t *option="")
Control routine to paint any kind of histograms.
Definition: TH1.cxx:5809
virtual Int_t GetSumw2N() const
Definition: TH1.h:310
virtual void SetBins(Int_t nx, Double_t xmin, Double_t xmax)
Redefine x axis parameters.
Definition: TH1.cxx:8181
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8433
static Bool_t AddDirectoryStatus()
Static function: cannot be inlined on Windows/NT.
Definition: TH1.cxx:705
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance from point px,py to a line.
Definition: TH1.cxx:2719
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8404
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:289
2-D histogram with a float per channel (see TH1 documentation)}
Definition: TH2.h:248
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2440
3-D histogram with a double per channel (see TH1 documentation)}
Definition: TH3.h:302
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:819
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:467
virtual TObject * First() const
Return the first object in the list. Returns 0 when list is empty.
Definition: TList.cxx:656
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
@ kNotDeleted
object has not been deleted
Definition: TObject.h:78
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:357
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual void AppendPad(Option_t *option="")
Append graphics object to current pad.
Definition: TObject.cxx:105
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TObject.cxx:664
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
void ResetBit(UInt_t f)
Definition: TObject.h:171
@ 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:854
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
TString & Insert(Ssiz_t pos, const char *s)
Definition: TString.h:644
Ssiz_t First(char c) const
Find first occurrence of a character c.
Definition: TString.cxx:499
const char * Data() const
Definition: TString.h:364
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:687
Bool_t IsNull() const
Definition: TString.h:402
TString & Append(const char *cs)
Definition: TString.h:559
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
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
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:416
Definition: graph.py:1
auto * a
Definition: textangle.C:12