 ROOT   Reference Guide Rolke.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook -nodraw
4/// Example of the usage of the TRolke class
5/// The TRolke class computes the profile likelihood
6/// confidence limits for 7 different model assumptions
7/// on systematic/statistical uncertainties
8///
10///
11/// \macro_output
12/// \macro_code
13///
14/// \authors Jan Conrad. Johan Lundberg
15
16#include "TROOT.h"
17#include "TSystem.h"
18#include "TRolke.h"
19#include "Riostream.h"
20
21void Rolke()
22{
23 // variables used throughout the example
24 Double_t bm;
25 Double_t tau;
26 Int_t mid;
27 Int_t m;
28 Int_t z;
29 Int_t y;
30 Int_t x;
31 Double_t e;
32 Double_t em;
33 Double_t sde;
34 Double_t sdb;
35 Double_t b;
36
37 Double_t alpha; //Confidence Level
38
39 // make TRolke objects
40 TRolke tr; //
41
42 Double_t ul ; // upper limit
43 Double_t ll ; // lower limit
44
45//-----------------------------------------------
46// Model 1 assumes:
47//
48// Poisson uncertainty in the background estimate
49// Binomial uncertainty in the efficiency estimate
50//
51 cout << endl<<" ======================================================== " <<endl;
52 mid =1;
53 x = 5; // events in the signal region
54 y = 10; // events observed in the background region
55 tau = 2.5; // ratio between size of signal/background region
56 m = 100; // MC events have been produced (signal)
57 z = 50; // MC events have been observed (signal)
58
59 alpha=0.9; //Confidence Level
60
61 tr.SetCL(alpha);
62
63 tr.SetPoissonBkgBinomEff(x,y,z,tau,m);
64 tr.GetLimits(ll,ul);
65
66 cout << "For model 1: Poisson / Binomial" << endl;
67 cout << "the Profile Likelihood interval is :" << endl;
68 cout << "[" << ll << "," << ul << "]" << endl;
69
70//-----------------------------------------------
71// Model 2 assumes:
72//
73// Poisson uncertainty in the background estimate
74// Gaussian uncertainty in the efficiency estimate
75//
76 cout << endl<<" ======================================================== " <<endl;
77 mid =2;
78 y = 3 ; // events observed in the background region
79 x = 10 ; // events in the signal region
80 tau = 2.5; // ratio between size of signal/background region
81 em = 0.9; // measured efficiency
82 sde = 0.05; // standard deviation of efficiency
83 alpha =0.95; // Confidence L evel
84
85 tr.SetCL(alpha);
86
87 tr.SetPoissonBkgGaussEff(x,y,em,tau,sde);
88 tr.GetLimits(ll,ul);
89
90 cout << "For model 2 : Poisson / Gaussian" << endl;
91 cout << "the Profile Likelihood interval is :" << endl;
92 cout << "[" << ll << "," << ul << "]" << endl;
93
94//-----------------------------------------------
95// Model 3 assumes:
96//
97// Gaussian uncertainty in the background estimate
98// Gaussian uncertainty in the efficiency estimate
99//
100 cout << endl<<" ======================================================== " <<endl;
101 mid =3;
102 bm = 5; // expected background
103 x = 10; // events in the signal region
104 sdb = 0.5; // standard deviation in background estimate
105 em = 0.9; // measured efficiency
106 sde = 0.05; // standard deviation of efficiency
107 alpha =0.99; // Confidence Level
108
109 tr.SetCL(alpha);
110
111 tr.SetGaussBkgGaussEff(x,bm,em,sde,sdb);
112 tr.GetLimits(ll,ul);
113 cout << "For model 3 : Gaussian / Gaussian" << endl;
114 cout << "the Profile Likelihood interval is :" << endl;
115 cout << "[" << ll << "," << ul << "]" << endl;
116
117 cout << "***************************************" << endl;
118 cout << "* some more example's for gauss/gauss *" << endl;
119 cout << "* *" << endl;
120 Double_t slow,shigh;
121 tr.GetSensitivity(slow,shigh);
122 cout << "sensitivity:" << endl;
123 cout << "[" << slow << "," << shigh << "]" << endl;
124
125 int outx;
126 tr.GetLimitsQuantile(slow,shigh,outx,0.5);
127 cout << "median limit:" << endl;
128 cout << "[" << slow << "," << shigh << "] @ x =" << outx <<endl;
129
130 tr.GetLimitsML(slow,shigh,outx);
131 cout << "ML limit:" << endl;
132 cout << "[" << slow << "," << shigh << "] @ x =" << outx <<endl;
133
134 tr.GetSensitivity(slow,shigh);
135 cout << "sensitivity:" << endl;
136 cout << "[" << slow << "," << shigh << "]" << endl;
137
138 tr.GetLimits(ll,ul);
139 cout << "the Profile Likelihood interval is :" << endl;
140 cout << "[" << ll << "," << ul << "]" << endl;
141
142 Int_t ncrt;
143
144 tr.GetCriticalNumber(ncrt);
145 cout << "critical number: " << ncrt << endl;
146
147 tr.SetCLSigmas(5);
148 tr.GetCriticalNumber(ncrt);
149 cout << "critical number for 5 sigma: " << ncrt << endl;
150
151 cout << "***************************************" << endl;
152
153//-----------------------------------------------
154// Model 4 assumes:
155//
156// Poisson uncertainty in the background estimate
157// known efficiency
158//
159 cout << endl<<" ======================================================== " <<endl;
160 mid =4;
161 y = 7; // events observed in the background region
162 x = 1; // events in the signal region
163 tau = 5; // ratio between size of signal/background region
164 e = 0.25; // efficiency
165
166 alpha =0.68; // Confidence L evel
167
168 tr.SetCL(alpha);
169
170 tr.SetPoissonBkgKnownEff(x,y,tau,e);
171 tr.GetLimits(ll,ul);
172
173 cout << "For model 4 : Poissonian / Known" << endl;
174 cout << "the Profile Likelihood interval is :" << endl;
175 cout << "[" << ll << "," << ul << "]" << endl;
176
177//-----------------------------------------------
178// Model 5 assumes:
179//
180// Gaussian uncertainty in the background estimate
181// Known efficiency
182//
183 cout << endl<<" ======================================================== " <<endl;
184 mid =5;
185 bm = 0; // measured background expectation
186 x = 1 ; // events in the signal region
187 e = 0.65; // known eff
188 sdb = 1.0; // standard deviation of background estimate
189 alpha =0.799999; // Confidence Level
190
191 tr.SetCL(alpha);
192
193 tr.SetGaussBkgKnownEff(x,bm,sdb,e);
194 tr.GetLimits(ll,ul);
195
196 cout << "For model 5 : Gaussian / Known" << endl;
197 cout << "the Profile Likelihood interval is :" << endl;
198 cout << "[" << ll << "," << ul << "]" << endl;
199
200//-----------------------------------------------
201// Model 6 assumes:
202//
203// Known background
204// Binomial uncertainty in the efficiency estimate
205//
206 cout << endl<<" ======================================================== " <<endl;
207 mid =6;
208 b = 10; // known background
209 x = 25; // events in the signal region
210 z = 500; // Number of observed signal MC events
211 m = 750; // Number of produced MC signal events
212 alpha =0.9; // Confidence L evel
213
214 tr.SetCL(alpha);
215
216 tr.SetKnownBkgBinomEff(x, z,m,b);
217 tr.GetLimits(ll,ul);
218
219 cout << "For model 6 : Known / Binomial" << endl;
220 cout << "the Profile Likelihood interval is :" << endl;
221 cout << "[" << ll << "," << ul << "]" << endl;
222
223//-----------------------------------------------
224// Model 7 assumes:
225//
226// Known Background
227// Gaussian uncertainty in the efficiency estimate
228//
229 cout << endl<<" ======================================================== " <<endl;
230 mid =7;
231 x = 15; // events in the signal region
232 em = 0.77; // measured efficiency
233 sde = 0.15; // standard deviation of efficiency estimate
234 b = 10; // known background
235 alpha =0.95; // Confidence L evel
236
237 y = 1;
238
239 tr.SetCL(alpha);
240
241 tr.SetKnownBkgGaussEff(x,em,sde,b);
242 tr.GetLimits(ll,ul);
243
244 cout << "For model 7 : Known / Gaussian " << endl;
245 cout << "the Profile Likelihood interval is :" << endl;
246 cout << "[" << ll << "," << ul << "]" << endl;
247
248//-----------------------------------------------
249// Example of bounded and unbounded likelihood
250// Example for Model 1
251
252 bm = 0.0;
253 tau = 5;
254 mid = 1;
255 m = 100;
256 z = 90;
257 y = 15;
258 x = 0;
259 alpha = 0.90;
260
261 tr.SetCL(alpha);
262 tr.SetPoissonBkgBinomEff(x,y,z,tau,m);
263 tr.SetBounding(true); //bounded
264 tr.GetLimits(ll,ul);
265
266 cout << "Example of the effect of bounded vs unbounded, For model 1" << endl;
267 cout << "the BOUNDED Profile Likelihood interval is :" << endl;
268 cout << "[" << ll << "," << ul << "]" << endl;
269
270
271 tr.SetBounding(false); //unbounded
272 tr.GetLimits(ll,ul);
273
274 cout << "the UNBOUNDED Profile Likelihood interval is :" << endl;
275 cout << "[" << ll << "," << ul << "]" << endl;
276}
#define b(i)
Definition: RSha256.hxx:100
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
double Double_t
Definition: RtypesCore.h:59
This class computes confidence intervals for the rate of a Poisson process in the presence of uncerta...
Definition: TRolke.h:34
void SetGaussBkgKnownEff(Int_t x, Double_t bm, Double_t sdb, Double_t e)
Model 5: Background - Gaussian, Efficiency - known (x,bm,sdb,e.
Definition: TRolke.cxx:302
bool GetLimitsML(Double_t &low, Double_t &high, Int_t &out_x)
get the upper and lower limits for the most likely outcome.
Definition: TRolke.cxx:511
bool GetCriticalNumber(Int_t &ncrit, Int_t maxtry=-1)
get the value of x corresponding to rejection of the null hypothesis.
Definition: TRolke.cxx:546
void SetKnownBkgBinomEff(Int_t x, Int_t z, Int_t m, Double_t b)
Model 6: Background - known, Efficiency - Binomial (x,z,m,b)
Definition: TRolke.cxx:327
void SetPoissonBkgKnownEff(Int_t x, Int_t y, Double_t tau, Double_t e)
Model 4: Background - Poisson, Efficiency - known (x,y,tau,e)
Definition: TRolke.cxx:277
void SetBounding(const bool bnd)
Definition: TRolke.h:184
void SetKnownBkgGaussEff(Int_t x, Double_t em, Double_t sde, Double_t b)
Model 7: Background - known, Efficiency - Gaussian (x,em,sde,b)
Definition: TRolke.cxx:352
void SetGaussBkgGaussEff(Int_t x, Double_t bm, Double_t em, Double_t sde, Double_t sdb)
Model 3: Background - Gaussian, Efficiency - Gaussian (x,bm,em,sde,sdb)
Definition: TRolke.cxx:252
void SetPoissonBkgGaussEff(Int_t x, Int_t y, Double_t em, Double_t tau, Double_t sde)
Model 2: Background - Poisson, Efficiency - Gaussian.
Definition: TRolke.cxx:226
bool GetSensitivity(Double_t &low, Double_t &high, Double_t pPrecision=0.00001)
get the upper and lower average limits based on the specified model.
Definition: TRolke.cxx:446
bool GetLimitsQuantile(Double_t &low, Double_t &high, Int_t &out_x, Double_t integral=0.5)
get the upper and lower limits for the outcome corresponding to a given quantile.
Definition: TRolke.cxx:481
bool GetLimits(Double_t &low, Double_t &high)
Calculate and get the upper and lower limits for the pre-specified model.
Definition: TRolke.cxx:373
void SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Model 1: Background - Poisson, Efficiency - Binomial.
Definition: TRolke.cxx:201
void SetCL(Double_t CL)
Definition: TRolke.h:124
void SetCLSigmas(Double_t CLsigmas)
Definition: TRolke.h:129
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
auto * m
Definition: textangle.C:8