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