Logo ROOT   6.18/05
Reference Guide
TRolke.cxx
Go to the documentation of this file.
1// @(#)root/physics:$Id$
2// Author: Jan Conrad
3
4/** \class TRolke
5 \ingroup Physics
6 This class computes confidence intervals for the rate of a Poisson
7 process in the presence of uncertain background and/or efficiency.
8
9 The treatment and the resulting limits are fully frequentist. The
10 limit calculations make use of the profile likelihood method.
11
12\author Jan Conrad (CERN) 2004, Updated: Johan Lundberg (CERN) 2009
13
14 For a full list of methods and their syntax, and build instructions,
15 consult the header file TRolke.h.
16
17 Examples/tutorials are found in the separate file Rolke.C
18
19### TRolke implements the following Models
20
21 The signal is always assumed to be Poisson, with the following
22 combinations of models of background and detection efficiency:
23
24 If unsure, first consider model 3, 4 or 5.
25
261: SetPoissonBkgBinomEff(x,y,z,tau,m)
27~~~
28 Background: Poisson
29 Efficiency: Binomial
30~~~
31 when the background is simultaneously measured
32 from sidebands (or MC), and
33 the signal efficiency was determined from Monte Carlo
34
352: SetPoissonBkgGaussEff(x,y,em,sde,tau)
36~~~
37 Background: Poisson
38 Efficiency: Gaussian
39~~~
40 when the background is simultaneously measured
41 from sidebands (or MC), and
42 the efficiency is modeled as Gaussian
43
443: SetGaussBkgGaussEff(x,bm,em,sde,sdb)
45~~~
46 Background: Gaussian
47 Efficiency: Gaussian
48~~~
49 when background and efficiency can both be
50 modeled as Gaussian.
51
524: SetPoissonBkgKnownEff(x,y,tau,e)
53~~~
54 Background: Poisson
55 Efficiency: Known
56~~~
57 when the background is simultaneously measured
58 from sidebands (or MC).
59
605: SetGaussBkgKnownEff(x,bm,sdb,e)
61~~~
62 Background: Gaussian
63 Efficiency: Known
64~~~
65 when background is Gaussian
66
676: SetKnownBkgBinomEff(x,z,b,m)
68~~~
69 Background: Known
70 Efficiency: Binomial
71~~~
72 when signal efficiency was determined from Monte Carlo
73
747: SetKnownBkgGaussEff(x,em,sde,b)
75~~~
76 Background: Known
77 Efficiency: Gaussian
78~~~
79 when background is known and efficiency Gaussian
80
81### Parameters and further explanation
82
83#### For all models:
84~~~
85 x = number of observed events in the experiment
86~~~
87 Efficiency (e or em) is the detection probability for signal.
88 A low efficiency hence generally means weaker limits.
89 If the efficiency of an experiment (with analysis cuts) is
90 dealt with elsewhere, em or e can be set to one.
91
92#### For Poisson background measurements (sideband or MC):
93~~~
94 y = number of observed events in background region
95 tau =
96 Either: the ratio between signal and background region
97 in case background is observed.
98 Or: the ratio between observed and simulated live-time
99 in case background is determined from MC.
100~~~
101#### For Gaussian efficiency or background:
102~~~
103 bm = estimate of the background
104 sdb = corresponding standard deviation
105
106 em = estimate of the efficiency
107 sde = corresponding standard deviation
108~~~
109 If the efficiency scale of dealt with elsewhere,
110 set em to 1 and sde to the relative uncertainty.
111
112#### For Binomial signal efficiency:
113~~~
114 m = number of MC events generated
115 z = number of MC events observed
116~~~
117#### For the case of known background expectation or known efficiency:
118~~~
119 e = true efficiency (considered known)
120 b = background expectation value (considered known)
121~~~
122
123 The confidence level (CL) is set either at construction
124 time or with either of SetCL or SetCLSigmas
125
126 The TRolke method is very similar to the one used in MINUIT (MINOS).
127
128 Two options are offered to deal with cases where the maximum likelihood
129 estimate (MLE) is not in the physical region. Version "bounded likelihood"
130 is the one used by MINOS if bounds for the physical region are chosen.
131 Unbounded likelihood (the default) allows the MLE to be in the
132 unphysical region. It has however better coverage.
133 For more details consult the reference (see below).
134
135 For a description of the method and its properties:
136
137 W.Rolke, A. Lopez, J. Conrad and Fred James
138 "Limits and Confidence Intervals in presence of nuisance parameters"
139 http://lanl.arxiv.org/abs/physics/0403059
140 Nucl.Instrum.Meth.A551:493-503,2005
141
142#### Should I use TRolke, TFeldmanCousins, TLimit?
143
144 1. Does TRolke make TFeldmanCousins obsolete?
145 Certainly not. TFeldmanCousins is the fully frequentist construction and
146 should be used in case of no (or negligible) uncertainties. It is however
147 not capable of treating uncertainties in nuisance parameters. In other
148 words, it does not handle background expectations or signal efficiencies
149 which are known only with some limited accuracy.
150 TRolke is designed for this case and it is shown in the reference above
151 that it has good coverage properties for most cases, and can be used
152 where FeldmannCousins can't.
153
154 2. What are the advantages of TRolke over TLimit?
155 TRolke is fully frequentist. TLimit treats nuisance parameters Bayesian.
156 For a coverage study of a Bayesian method refer to
157 physics/0408039 (Tegenfeldt & J.C). However, this note studies
158 the coverage of Feldman&Cousins with Bayesian treatment of nuisance
159 parameters. To make a long story short: using the Bayesian method you
160 might introduce a small amount of over-coverage (though I haven't shown it
161 for TLimit). On the other hand, coverage of course is a not so interesting
162 when you consider yourself a Bayesian.
163*/
164
165#include "TRolke.h"
166#include "TMath.h"
167#include <iostream>
168
170
171////////////////////////////////////////////////////////////////////////////////
172/// Constructor with optional Confidence Level argument.
173/// 'option' is not used.
174
176: fCL(CL),
177 fUpperLimit(0.0),
178 fLowerLimit(0.0),
179 fBounding(false), // true gives bounded likelihood
180 fNumWarningsDeprecated1(0),
181 fNumWarningsDeprecated2(0)
182{
184}
185
186////////////////////////////////////////////////////////////////////////////////
187/// Destructor.
188
190{
191}
192
193////////////////////////////////////////////////////////////////////////////////
194/// Model 1: Background - Poisson, Efficiency - Binomial
195/// - x : number of observed events in the experiment
196/// - y : number of observed events in background region
197/// - z : number of MC events observed
198/// - tau : ratio parameter (read TRolke.cxx for details)
199/// - m : number of MC events generated
200
202{
204 x , // Int_t x,
205 y , // Int_t y,
206 z , // Int_t z,
207 0 , // Double_t bm,
208 0 , // Double_t em,
209 0 , // Double_t e,
210 1 , // Int_t mid,
211 0 , // Double_t sde,
212 0 , // Double_t sdb,
213 tau, // Double_t tau,
214 0 , // Double_t b,
215 m); // Int_t m
216}
217
218////////////////////////////////////////////////////////////////////////////////
219/// Model 2: Background - Poisson, Efficiency - Gaussian
220/// - x : number of observed events in the experiment
221/// - y : number of observed events in background region
222/// - em : estimate of the efficiency
223/// - tau : ratio parameter (read TRolke.cxx for details)
224/// - sde : efficiency estimate's standard deviation
225
227{
229 x , // Int_t x,
230 y , // Int_t y,
231 0 , // Int_t z,
232 0 , // Double_t bm,
233 em , // Double_t em,
234 0 , // Double_t e,
235 2 , // Int_t mid,
236 sde, // Double_t sde,
237 0 , // Double_t sdb,
238 tau, // Double_t tau,
239 0 , // Double_t b,
240 0); // Int_t m
241
242}
243
244////////////////////////////////////////////////////////////////////////////////
245/// Model 3: Background - Gaussian, Efficiency - Gaussian (x,bm,em,sde,sdb)
246/// - x : number of observed events in the experiment
247/// - bm : estimate of the background
248/// - em : estimate of the efficiency
249/// - sde : efficiency estimate's standard deviation
250/// - sdb : background estimate's standard deviation
251
253{
255 x , // Int_t x,
256 0 , // Int_t y,
257 0 , // Int_t z,
258 bm , // Double_t bm,
259 em , // Double_t em,
260 0 , // Double_t e,
261 3 , // Int_t mid,
262 sde, // Double_t sde,
263 sdb, // Double_t sdb,
264 0 , // Double_t tau,
265 0 , // Double_t b,
266 0); // Int_t m
267
268}
269
270////////////////////////////////////////////////////////////////////////////////
271/// Model 4: Background - Poisson, Efficiency - known (x,y,tau,e)
272/// - x : number of observed events in the experiment
273/// - y : number of observed events in background region
274/// - tau : ratio parameter (read TRolke.cxx for details)
275/// - e : true efficiency (considered known)
276
278{
280 x , // Int_t x,
281 y , // Int_t y,
282 0 , // Int_t z,
283 0 , // Double_t bm,
284 0 , // Double_t em,
285 e , // Double_t e,
286 4 , // Int_t mid,
287 0 , // Double_t sde,
288 0 , // Double_t sdb,
289 tau, // Double_t tau,
290 0 , // Double_t b,
291 0); // Int_t m
292
293}
294
295////////////////////////////////////////////////////////////////////////////////
296/// Model 5: Background - Gaussian, Efficiency - known (x,bm,sdb,e
297/// - x : number of observed events in the experiment
298/// - bm : estimate of the background
299/// - sdb : background estimate's standard deviation
300/// - e : true efficiency (considered known)
301
303{
305 x , // Int_t x,
306 0 , // Int_t y,
307 0 , // Int_t z,
308 bm , // Double_t bm,
309 0 , // Double_t em,
310 e , // Double_t e,
311 5 , // Int_t mid,
312 0 , // Double_t sde,
313 sdb, // Double_t sdb,
314 0 , // Double_t tau,
315 0 , // Double_t b,
316 0); // Int_t m
317
318}
319
320////////////////////////////////////////////////////////////////////////////////
321/// Model 6: Background - known, Efficiency - Binomial (x,z,m,b)
322/// - x : number of observed events in the experiment
323/// - z : number of MC events observed
324/// - m : number of MC events generated
325/// - b : background expectation value (considered known)
326
328{
330 x , // Int_t x,
331 0 , // Int_t y
332 z , // Int_t z,
333 0 , // Double_t bm,
334 0 , // Double_t em,
335 0 , // Double_t e,
336 6 , // Int_t mid,
337 0 , // Double_t sde,
338 0 , // Double_t sdb,
339 0 , // Double_t tau,
340 b , // Double_t b,
341 m); // Int_t m
342
343}
344
345////////////////////////////////////////////////////////////////////////////////
346/// Model 7: Background - known, Efficiency - Gaussian (x,em,sde,b)
347/// - x : number of observed events in the experiment
348/// - em : estimate of the efficiency
349/// - sde : efficiency estimate's standard deviation
350/// - b : background expectation value (considered known)
351
353{
355 x , // Int_t x,
356 0 , // Int_t y
357 0 , // Int_t z,
358 0 , // Double_t bm,
359 em , // Double_t em,
360 0 , // Double_t e,
361 7 , // Int_t mid,
362 sde, // Double_t sde,
363 0 , // Double_t sdb,
364 0 , // Double_t tau,
365 b , // Double_t b,
366 0); // Int_t m
367
368}
369
370////////////////////////////////////////////////////////////////////////////////
371/// Calculate and get the upper and lower limits for the pre-specified model.
372
374{
375 if ((f_mid<1)||(f_mid>7)) {
376 std::cerr << "TRolke - Error: Model id "<< f_mid<<std::endl;
377 if (f_mid<1) {
378 std::cerr << "TRolke - Please specify a model with e.g. 'SetGaussBkgGaussEff' (read the docs in Rolke.cxx )"<<std::endl;
379 }
380 return false;
381 }
382
384 low = fLowerLimit;
385 high = fUpperLimit;
386 if (low < high) {
387 return true;
388 }else{
389 std::cerr << "TRolke - Warning: no limits found" <<std::endl;
390 return false;
391 }
392}
393
394////////////////////////////////////////////////////////////////////////////////
395/// Calculate and get upper limit for the pre-specified model.
396
398{
399 Double_t low(0), high(0);
400 GetLimits(low,high);
401 return fUpperLimit;
402}
403
404////////////////////////////////////////////////////////////////////////////////
405/// Calculate and get lower limit for the pre-specified model.
406
408{
409 Double_t low(0), high(0);
410 GetLimits(low,high);
411 return fLowerLimit;
412}
413
414////////////////////////////////////////////////////////////////////////////////
415/// Return a simple background value (estimate/truth) given the pre-specified model.
416
418{
419 Double_t background = 0;
420 switch (f_mid) {
421 case 1:
422 case 2:
423 case 4:
424 background = f_y / f_tau;
425 break;
426 case 3:
427 case 5:
428 background = f_bm;
429 break;
430 case 6:
431 case 7:
432 background = f_b;
433 break;
434 default:
435 std::cerr << "TRolke::GetBackground(): Model NR: " <<
436 f_mid << " unknown"<<std::endl;
437 return 0;
438 }
439 return background;
440}
441
442////////////////////////////////////////////////////////////////////////////////
443/// get the upper and lower average limits based on the specified model.
444/// No uncertainties are considered for the Poisson weights in the averaging sum
445
447{
448 Double_t background = GetBackground();
449
450 Double_t weight = 0;
451 Double_t weightSum = 0;
452
453 int loop_x = 0;
454
455 while (1) {
457 weight = TMath::PoissonI(loop_x, background);
458 low += fLowerLimit * weight;
459 high += fUpperLimit * weight;
460 weightSum += weight;
461 if (loop_x > (background + 1)) { // don't stop too early
462 if ((weightSum > (1 - pPrecision)) || (weight < 1e-12)) break;
463 }
464 loop_x++;
465 }
466 low /= weightSum;
467 high /= weightSum;
468
469 return (low < high); // could also add more detailed test
470}
471
472////////////////////////////////////////////////////////////////////////////////
473/// get the upper and lower limits for the outcome corresponding to
474/// a given quantile.
475/// For integral=0.5 this gives the median limits
476/// in repeated experiments. The returned out_x is the corresponding
477/// (e.g. median) value of x.
478/// No uncertainties are considered for the Poisson weights when calculating
479/// the Poisson integral.
480
481bool TRolke::GetLimitsQuantile(Double_t& low, Double_t& high, Int_t& out_x, Double_t integral)
482{
483 Double_t background = GetBackground();
484 Double_t weight = 0;
485 Double_t weightSum = 0;
486 Int_t loop_x = 0;
487
488 while (1) {
489 weight = TMath::PoissonI(loop_x, background);
490 weightSum += weight;
491 if (weightSum >= integral) {
492 break;
493 }
494 loop_x++;
495 }
496
497 out_x = loop_x;
498
500 low = fLowerLimit;
501 high = fUpperLimit;
502 return (low < high); // could also add more detailed test
503
504}
505
506////////////////////////////////////////////////////////////////////////////////
507/// get the upper and lower limits for the most likely outcome.
508/// The returned out_x is the corresponding value of x
509/// No uncertainties are considered for the Poisson weights when finding ML.
510
512{
513 Double_t background = GetBackground();
514
515 Int_t loop_x = 0; // this can be optimized if needed.
516 Int_t loop_max = 1000 + (Int_t)background; // --||--
517
518 Double_t max = TMath::PoissonI(loop_x, background);
519 while (loop_x <= loop_max) {
520 if (TMath::PoissonI(loop_x + 1, background) < max) {
521 break;
522 }
523 loop_x++;
524 max = TMath::PoissonI(loop_x, background);
525 }
526 if (loop_x >= loop_max) {
527 std::cout << "internal error finding maximum of distribution" << std::endl;
528 return false;
529 }
530
531 out_x = loop_x;
532
534 low = fLowerLimit;
535 high = fUpperLimit;
536 return (low < high); // could also add more detailed test
537
538}
539
540////////////////////////////////////////////////////////////////////////////////
541/// get the value of x corresponding to rejection of the null hypothesis.
542/// This means a lower limit >0 with the pre-specified Confidence Level.
543/// Optionally give maxtry; the maximum value of x to try. Of not, or if
544/// maxtry<0 an automatic mode is used.
545
547{
548 Double_t background = GetBackground();
549
550 int j = 0;
551 int rolke_ncrit = -1;
552 int maxj =maxtry ;
553 if(maxtry<1){
554 maxj = 1000 + (Int_t)background; // max value to try
555 }
556 for (j = 0;j < maxj;j++) {
557 Int_t rolke_x = j;
559 double rolke_ll = fLowerLimit;
560 if (rolke_ll > 0) {
561 rolke_ncrit = j;
562 break;
563 }
564 }
565
566 if (rolke_ncrit == -1) {
567 std::cerr << "TRolke GetCriticalNumber : Error: problem finding rolke inverse. Specify a larger maxtry value. maxtry was: " << maxj << ". highest x considered was j "<< j<< std::endl;
568 ncrit = -1;
569 return false;
570 } else {
571 ncrit = rolke_ncrit;
572 return true;
573 }
574}
575
576////////////////////////////////////////////////////////////////////////////////
577/// Deprecated name for SetBounding.
578
579void TRolke::SetSwitch(bool bnd) {
581 std::cerr << "*******************************************" <<std::endl;
582 std::cerr << "TRolke - Warning: 'SetSwitch' is depricated and may be removed from future releases:" <<std::endl;
583 std::cerr << " - Use 'SetBounding' instead "<<std::endl;
584 std::cerr << "*******************************************" <<std::endl;
586 }
587 SetBounding(bnd);
588}
589
590////////////////////////////////////////////////////////////////////////////////
591/// Dump internals. Print members.
592
594 std::cout << "*******************************************" <<std::endl;
595 std::cout << "* TRolke::Print() - dump of internals: " <<std::endl;
596 std::cout << "*"<<std::endl;
597 std::cout << "* model id, mid = "<<f_mid <<std::endl;
598 std::cout << "*"<<std::endl;
599 std::cout << "* x = "<<f_x <<std::endl;
600 std::cout << "* bm = "<<f_bm <<std::endl;
601 std::cout << "* em = "<<f_em <<std::endl;
602 std::cout << "* sde = "<<f_sde <<std::endl;
603 std::cout << "* sdb = "<<f_sdb <<std::endl;
604 std::cout << "* y = "<<f_y <<std::endl;
605 std::cout << "* tau = "<<f_tau <<std::endl;
606 std::cout << "* e = "<<f_e <<std::endl;
607 std::cout << "* b = "<<f_b <<std::endl;
608 std::cout << "* m = "<<f_m <<std::endl;
609 std::cout << "* z = "<<f_z <<std::endl;
610 std::cout << "*"<<std::endl;
611 std::cout << "* CL = "<<fCL <<std::endl;
612 std::cout << "* Bounding = "<<fBounding <<std::endl;
613 std::cout << "*"<<std::endl;
614 std::cout << "* calculated on demand only:"<<std::endl;
615 std::cout << "* fUpperLimit = "<<fUpperLimit<<std::endl;
616 std::cout << "* fLowerLimit = "<<fLowerLimit<<std::endl;
617 std::cout << "*******************************************" <<std::endl;
618}
619
620////////////////////////////////////////////////////////////////////////////////
621/// Deprecated and error prone model selection interface.
622/// It's use is trongly discouraged. 'mid' is the model ID (1 to 7).
623/// This method is provided for backwards compatibility/developer use only. */
624/// - x : number of observed events in the experiment
625/// - y : number of observed events in background region
626/// - z : number of MC events observed
627/// - bm : estimate of the background
628/// - em : estimate of the efficiency
629/// - e : true efficiency (considered known)
630/// - mid : internal model id (really, you should not use this method at all)
631/// - sde : efficiency estimate's standard deviation
632/// - sdb : background estimate's standard deviation
633/// - tau : ratio parameter (read TRolke.cxx for details)
634/// - b : background expectation value (considered known)
635/// - m : number of MC events generated
636
639 std::cerr << "*******************************************" <<std::endl;
640 std::cerr << "TRolke - Warning: 'CalculateInterval' is depricated and may be removed from future releases:" <<std::endl;
641 std::cerr << " - Use e.g. 'SetGaussBkgGaussEff' and 'GetLimits' instead (read the docs in Rolke.cxx )"<<std::endl;
642 std::cerr << "*******************************************" <<std::endl;
644 }
646 x,
647 y,
648 z,
649 bm,
650 em,
651 e,
652 mid,
653 sde,
654 sdb,
655 tau,
656 b,
657 m);
659}
660
661////////////////////////////////////////////////////////////////////////////////
662/// - x : number of observed events in the experiment
663/// - y : number of observed events in background region
664/// - z : number of MC events observed
665/// - bm : estimate of the background
666/// - em : estimate of the efficiency
667/// - e : true efficiency (considered known)
668/// - mid : internal model id
669/// - sde : efficiency estimate's standard deviation
670/// - sdb : background estimate's standard deviation
671/// - tau : ratio parameter (read TRolke.cxx for details)
672/// - b : background expectation value (considered known)
673/// - m : number of MC events generated
674
676{
677 f_x = x ;
678 f_y = y ;
679 f_z = z ;
680 f_bm = bm ;
681 f_em = em ;
682 f_e = e ;
683 f_mid = mid ;
684 f_sde = sde ;
685 f_sdb = sdb ;
686 f_tau = tau ;
687 f_b = b ;
688 f_m = m ;
689}
690
692{
693/* Clear internal model */
694 SetModelParameters(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
695 f_mid=0;
696}
697
698////////////////////////////////////////////////////////////////////////////////
699/// ComputeInterval, the internals.
700/// - x : number of observed events in the experiment
701/// - y : number of observed events in background region
702/// - z : number of MC events observed
703/// - bm : estimate of the background
704/// - em : estimate of the efficiency
705/// - e : true efficiency (considered known)
706/// - mid : internal model id (really, you should not use this method at all)
707/// - sde : efficiency estimate's standard deviation
708/// - sdb : background estimate's standard deviation
709/// - tau : ratio parameter (read TRolke.cxx for details)
710/// - b : background expectation value (considered known)
711/// - m : number of MC events generated
712
714{
715 //calculate interval
716 Int_t done = 0;
717 Double_t limit[2];
718
719 limit[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
720
721 if (limit[1] > 0) {
722 done = 1;
723 }
724
725 if (! fBounding) {
726
727 Int_t trial_x = x;
728
729 while (done == 0) {
730 trial_x++;
731 limit[1] = Interval(trial_x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
732 if (limit[1] > 0) done = 1;
733 }
734 }
735
736 return limit[1];
737}
738
739////////////////////////////////////////////////////////////////////////////////
740/// Internal helper function 'Interval'
741/// - x : number of observed events in the experiment
742/// - y : number of observed events in background region
743/// - z : number of MC events observed
744/// - bm : estimate of the background
745/// - em : estimate of the efficiency
746/// - e : true efficiency (considered known)
747/// - mid : internal model id (really, you should not use this method at all)
748/// - sde : efficiency estimate's standard deviation
749/// - sdb : background estimate's standard deviation
750/// - tau : ratio parameter (read TRolke.cxx for details)
751/// - b : background expectation value (considered known)
752/// - m : number of MC events generated
753
755{
757 Double_t tempxy[2], limits[2] = {0, 0};
758 Double_t slope, fmid, low, flow, high, fhigh, test, ftest, mu0, maximum, target, l, f0;
759 Double_t med = 0;
760 Double_t maxiter = 1000, acc = 0.00001;
761 Int_t i;
762 Int_t bp = 0;
763
764 if ((mid != 3) && (mid != 5)) bm = y;
765 if ((mid == 3) || (mid == 5)) {
766 if (bm == 0) bm = 0.00001;
767 }
768
769 if ((mid == 6) || (mid == 7)) {
770 if (bm == 0) bm = 0.00001;
771 }
772
773 if ((mid <= 2) || (mid == 4)) bp = 1;
774
775
776 if (bp == 1 && x == 0 && bm > 0) {
777 for (i = 0; i < 2; i++) {
778 x++;
779 tempxy[i] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
780 }
781
782 slope = tempxy[1] - tempxy[0];
783 limits[1] = tempxy[0] - slope;
784 limits[0] = 0.0;
785 if (limits[1] < 0) limits[1] = 0.0;
786 goto done;
787 }
788
789 if (bp != 1 && x == 0) {
790
791 for (i = 0; i < 2; i++) {
792 x++;
793 tempxy[i] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
794 }
795 slope = tempxy[1] - tempxy[0];
796 limits[1] = tempxy[0] - slope;
797 limits[0] = 0.0;
798 if (limits[1] < 0) limits[1] = 0.0;
799 goto done;
800 }
801
802 if (bp != 1 && bm == 0) {
803 for (i = 0; i < 2; i++) {
804 bm++;
805 limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
806 tempxy[i] = limits[1];
807 }
808 slope = tempxy[1] - tempxy[0];
809 limits[1] = tempxy[0] - slope;
810 if (limits[1] < 0) limits[1] = 0;
811 goto done;
812 }
813
814 if (x == 0 && bm == 0) {
815 x++;
816 bm++;
817 limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
818 tempxy[0] = limits[1];
819 x = 1;
820 bm = 2;
821 limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
822 tempxy[1] = limits[1];
823 x = 2;
824 bm = 1;
825 limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
826 limits[1] = 3 * tempxy[0] - tempxy[1] - limits[1];
827 if (limits[1] < 0) limits[1] = 0;
828 goto done;
829 }
830
831 mu0 = Likelihood(0, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 1);
832 maximum = Likelihood(0, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 2);
833 test = 0;
834 f0 = Likelihood(test, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
835 if (fBounding) {
836 if (mu0 < 0) maximum = f0;
837 }
838
839 target = maximum - dchi2;
840 if (f0 > target) {
841 limits[0] = 0;
842 } else {
843 if (mu0 < 0) {
844 limits[0] = 0;
845 limits[1] = 0;
846 }
847
848 low = 0;
849 flow = f0;
850 high = mu0;
851 fhigh = maximum;
852 for (i = 0; i < maxiter; i++) {
853 l = (target - fhigh) / (flow - fhigh);
854 if (l < 0.2) l = 0.2;
855 if (l > 0.8) l = 0.8;
856 med = l * low + (1 - l) * high;
857 if (med < 0.01) {
858 limits[1] = 0.0;
859 goto done;
860 }
861 fmid = Likelihood(med, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
862 if (fmid > target) {
863 high = med;
864 fhigh = fmid;
865 } else {
866 low = med;
867 flow = fmid;
868 }
869 if ((high - low) < acc*high) break;
870 }
871 limits[0] = med;
872 }
873
874 if (mu0 > 0) {
875 low = mu0;
876 flow = maximum;
877 } else {
878 low = 0;
879 flow = f0;
880 }
881
882 test = low + 1 ;
883 ftest = Likelihood(test, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
884 if (ftest < target) {
885 high = test;
886 fhigh = ftest;
887 } else {
888 slope = (ftest - flow) / (test - low);
889 high = test + (target - ftest) / slope;
890 fhigh = Likelihood(high, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
891 }
892
893 for (i = 0; i < maxiter; i++) {
894 l = (target - fhigh) / (flow - fhigh);
895 if (l < 0.2) l = 0.2;
896 if (l > 0.8) l = 0.8;
897 med = l * low + (1. - l) * high;
898 fmid = Likelihood(med, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
899
900 if (fmid < target) {
901 high = med;
902 fhigh = fmid;
903 } else {
904 low = med;
905 flow = fmid;
906 }
907
908 if (high - low < acc*high) break;
909 }
910
911 limits[1] = med;
912
913done:
914
915 // apply known efficiency
916 if ((mid == 4) || (mid == 5)) {
917 limits[0] /= e;
918 limits[1] /= e;
919 }
920
921 fUpperLimit = limits[1];
922 fLowerLimit = TMath::Max(limits[0], 0.0);
923
924 return limits[1];
925}
926
927////////////////////////////////////////////////////////////////////////////////
928/// Internal helper function
929/// Chooses between the different profile likelihood functions to use for the
930/// different models.
931/// Returns evaluation of the profile likelihood functions.
932
934{
935 switch (mid) {
936 case 1:
937 return EvalLikeMod1(mu, x, y, z, tau, m, what);
938 case 2:
939 return EvalLikeMod2(mu, x, y, em, sde, tau, what);
940 case 3:
941 return EvalLikeMod3(mu, x, bm, em, sde, sdb, what);
942 case 4:
943 return EvalLikeMod4(mu, x, y, tau, what);
944 case 5:
945 return EvalLikeMod5(mu, x, bm, sdb, what);
946 case 6:
947 return EvalLikeMod6(mu, x, z, b, m, what);
948 case 7:
949 return EvalLikeMod7(mu, x, em, sde, b, what);
950 default:
951 std::cerr << "TRolke::Likelihood(...): Model NR: " <<
952 f_mid << " unknown"<<std::endl;
953 return 0;
954 }
955
956 return 0;
957}
958
959////////////////////////////////////////////////////////////////////////////////
960/// Calculates the Profile Likelihood for MODEL 1:
961/// Poisson background/ Binomial Efficiency
962/// - what = 1: Maximum likelihood estimate is returned
963/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
964/// - what = 3: Profile Likelihood of Test hypothesis is returned
965/// otherwise parameters as described in the beginning of the class)
966
968{
969 Double_t f = 0;
970 Double_t zm = Double_t(z) / m;
971
972 if (what == 1) {
973 f = (x - y / tau) / zm;
974 }
975
976 if (what == 2) {
977 mu = (x - y / tau) / zm;
978 Double_t b = y / tau;
979 Double_t e = zm;
980 f = LikeMod1(mu, b, e, x, y, z, tau, m);
981 }
982
983 if (what == 3) {
984 if (mu == 0) {
985 Double_t b = (x + y) / (1.0 + tau);
986 Double_t e = zm;
987 f = LikeMod1(mu, b, e, x, y, z, tau, m);
988 } else {
989 Double_t e = 0;
990 Double_t b = 0;
991 ProfLikeMod1(mu, b, e, x, y, z, tau, m);
992 f = LikeMod1(mu, b, e, x, y, z, tau, m);
993 }
994 }
995
996 return f;
997}
998
999////////////////////////////////////////////////////////////////////////////////
1000/// Profile Likelihood function for MODEL 1:
1001/// Poisson background/ Binomial Efficiency
1002
1004{
1005 double s = e*mu+b;
1006 double lls = - s;
1007 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1008 double bg = tau*b;
1009 double llb = -bg;
1010 if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y);
1011
1012 double lle = 0; // binomial log-like
1013 if (z == 0) lle = m * TMath::Log(1-e);
1014 else if ( z == m) lle = m * TMath::Log(e);
1015 else lle = z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z);
1016
1017 double f = 2*( lls + llb + lle);
1018 return f;
1019}
1020
1021
1022// this code is non-sense - // need to solve using Minuit
1023struct LikeFunction1 {
1024};
1025
1026////////////////////////////////////////////////////////////////////////////////
1027/// Helper for calculation of estimates of efficiency and background for model 1
1028
1030{
1031 Double_t med = 0.0, fmid;
1032 Int_t maxiter = 1000;
1033 Double_t acc = 0.00001;
1034 Double_t emin = ((m + mu * tau) - TMath::Sqrt((m + mu * tau) * (m + mu * tau) - 4 * mu * tau * z)) / 2 / mu / tau;
1035
1036 Double_t low = TMath::Max(1e-10, emin + 1e-10);
1037 Double_t high = 1 - 1e-10;
1038
1039 for (Int_t i = 0; i < maxiter; i++) {
1040 med = (low + high) / 2.;
1041
1042 fmid = LikeGradMod1(med, mu, x, y, z, tau, m);
1043
1044 if (high < 0.5) acc = 0.00001 * high;
1045 else acc = 0.00001 * (1 - high);
1046
1047 if ((high - low) < acc*high) break;
1048
1049 if (fmid > 0) low = med;
1050 else high = med;
1051 }
1052
1053 e = med;
1054 Double_t eta = Double_t(z) / e - Double_t(m - z) / (1 - e);
1055
1056 b = Double_t(y) / (tau - eta / mu);
1057}
1058
1059////////////////////////////////////////////////////////////////////////////////
1060/// Gradient model likelihood
1061
1063{
1064 Double_t eta, etaprime, bprime, f;
1065 eta = static_cast<double>(z) / e - static_cast<double>(m - z) / (1.0 - e);
1066 etaprime = (-1) * (static_cast<double>(m - z) / ((1.0 - e) * (1.0 - e)) + static_cast<double>(z) / (e * e));
1067 Double_t b = y / (tau - eta / mu);
1068 bprime = (b * b * etaprime) / mu / y;
1069 f = (mu + bprime) * (x / (e * mu + b) - 1) + (y / b - tau) * bprime + eta;
1070 return f;
1071}
1072
1073////////////////////////////////////////////////////////////////////////////////
1074/// Calculates the Profile Likelihood for MODEL 2:
1075/// Poisson background/ Gauss Efficiency
1076/// - what = 1: Maximum likelihood estimate is returned
1077/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1078/// - what = 3: Profile Likelihood of Test hypothesis is returned
1079/// otherwise parameters as described in the beginning of the class)
1080
1082{
1083 Double_t v = sde * sde;
1084 Double_t coef[4], roots[3];
1085 Double_t f = 0;
1086
1087 if (what == 1) {
1088 f = (x - y / tau) / em;
1089 }
1090
1091 if (what == 2) {
1092 mu = (x - y / tau) / em;
1093 Double_t b = y / tau;
1094 Double_t e = em;
1095 f = LikeMod2(mu, b, e, x, y, em, tau, v);
1096 }
1097
1098 if (what == 3) {
1099 if (mu == 0) {
1100 Double_t b = (x + y) / (1 + tau);
1101 Double_t e = em ;
1102 f = LikeMod2(mu, b, e, x, y, em, tau, v);
1103 } else {
1104 coef[3] = mu;
1105 coef[2] = mu * mu * v - 2 * em * mu - mu * mu * v * tau;
1106 coef[1] = (- x) * mu * v - mu * mu * mu * v * v * tau - mu * mu * v * em + em * mu * mu * v * tau + em * em * mu - y * mu * v;
1107 coef[0] = x * mu * mu * v * v * tau + x * em * mu * v - y * mu * mu * v * v + y * em * mu * v;
1108
1109 TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);
1110
1111 Double_t e = roots[1];
1112 Double_t b;
1113 if ( v > 0) b = y / (tau + (em - e) / mu / v);
1114 else b = y/tau;
1115 f = LikeMod2(mu, b, e, x, y, em, tau, v);
1116 }
1117 }
1118
1119 return f;
1120}
1121
1122////////////////////////////////////////////////////////////////////////////////
1123/// Profile Likelihood function for MODEL 2:
1124/// Poisson background/Gauss Efficiency
1125
1127{
1128 double s = e*mu+b;
1129 double lls = - s;
1130 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1131 double bg = tau*b;
1132 double llb = -bg;
1133 if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y);
1134 double lle = 0;
1135 if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1136
1137 return 2*( lls + llb + lle);
1138}
1139
1140////////////////////////////////////////////////////////////////////////////////
1141/// Calculates the Profile Likelihood for MODEL 3:
1142/// Gauss background/ Gauss Efficiency
1143/// - what = 1: Maximum likelihood estimate is returned
1144/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1145/// - what = 3: Profile Likelihood of Test hypothesis is returned
1146/// otherwise parameters as described in the beginning of the class)
1147
1149{
1150 Double_t f = 0.;
1151 Double_t v = sde * sde;
1152 Double_t u = sdb * sdb;
1153
1154 if (what == 1) {
1155 f = (x - bm) / em;
1156 }
1157
1158
1159 if (what == 2) {
1160 mu = (x - bm) / em;
1161 Double_t b = bm;
1162 Double_t e = em;
1163 f = LikeMod3(mu, b, e, x, bm, em, u, v);
1164 }
1165
1166
1167 if (what == 3) {
1168 if (mu == 0.0) {
1169 Double_t b = ((bm - u) + TMath::Sqrt((bm - u) * (bm - u) + 4 * x * u)) / 2.;
1170 Double_t e = em;
1171 f = LikeMod3(mu, b, e, x, bm, em, u, v);
1172 } else {
1173 Double_t e = em;
1174 Double_t b = bm;
1175 if ( v > 0) {
1176 Double_t temp[3];
1177 temp[0] = mu * mu * v + u;
1178 temp[1] = mu * mu * mu * v * v + mu * v * u - mu * mu * v * em + mu * v * bm - 2 * u * em;
1179 temp[2] = mu * mu * v * v * bm - mu * v * u * em - mu * v * bm * em + u * em * em - mu * mu * v * v * x;
1180 e = (-temp[1] + TMath::Sqrt(temp[1] * temp[1] - 4 * temp[0] * temp[2])) / 2 / temp[0];
1181 b = bm - (u * (em - e)) / v / mu;
1182 }
1183 f = LikeMod3(mu, b, e, x, bm, em, u, v);
1184 }
1185 }
1186
1187 return f;
1188}
1189
1190////////////////////////////////////////////////////////////////////////////////
1191/// Profile Likelihood function for MODEL 3:
1192/// Gauss background/Gauss Efficiency
1193
1195{
1196 double s = e*mu+b;
1197 double lls = - s;
1198 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1199 double llb = 0;
1200 if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;
1201 double lle = 0;
1202 if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1203
1204 return 2*( lls + llb + lle);
1205
1206}
1207
1208////////////////////////////////////////////////////////////////////////////////
1209/// Calculates the Profile Likelihood for MODEL 4:
1210/// Poiss background/Efficiency known
1211/// - what = 1: Maximum likelihood estimate is returned
1212/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1213/// - what = 3: Profile Likelihood of Test hypothesis is returned
1214/// otherwise parameters as described in the beginning of the class)
1215
1217{
1218 Double_t f = 0.0;
1219
1220 if (what == 1) f = x - y / tau;
1221 if (what == 2) {
1222 mu = x - y / tau;
1223 Double_t b = y / tau;
1224 f = LikeMod4(mu, b, x, y, tau);
1225 }
1226 if (what == 3) {
1227 if (mu == 0.0) {
1228 Double_t b = Double_t(x + y) / (1 + tau);
1229 f = LikeMod4(mu, b, x, y, tau);
1230 } else {
1231 Double_t b = (x + y - (1 + tau) * mu + sqrt((x + y - (1 + tau) * mu) * (x + y - (1 + tau) * mu) + 4 * (1 + tau) * y * mu)) / 2 / (1 + tau);
1232 f = LikeMod4(mu, b, x, y, tau);
1233 }
1234 }
1235 return f;
1236}
1237
1238////////////////////////////////////////////////////////////////////////////////
1239/// Profile Likelihood function for MODEL 4:
1240/// Poiss background/Efficiency known
1241
1243{
1244 double s = mu+b;
1245 double lls = - s;
1246 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1247 double bg = tau*b;
1248 double llb = -bg;
1249 if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y);
1250
1251 return 2*( lls + llb);
1252}
1253
1254////////////////////////////////////////////////////////////////////////////////
1255/// Calculates the Profile Likelihood for MODEL 5:
1256/// Gauss background/Efficiency known
1257/// - what = 1: Maximum likelihood estimate is returned
1258/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1259/// - what = 3: Profile Likelihood of Test hypothesis is returned
1260/// otherwise parameters as described in the beginning of the class)
1261
1263{
1264 Double_t u = sdb * sdb;
1265 Double_t f = 0;
1266
1267 if (what == 1) {
1268 f = x - bm;
1269 }
1270 if (what == 2) {
1271 mu = x - bm;
1272 Double_t b = bm;
1273 f = LikeMod5(mu, b, x, bm, u);
1274 }
1275
1276 if (what == 3) {
1277 Double_t b = ((bm - u - mu) + TMath::Sqrt((bm - u - mu) * (bm - u - mu) - 4 * (mu * u - mu * bm - u * x))) / 2;
1278 f = LikeMod5(mu, b, x, bm, u);
1279 }
1280 return f;
1281}
1282
1283////////////////////////////////////////////////////////////////////////////////
1284/// Profile Likelihood function for MODEL 5:
1285/// Gauss background/Efficiency known
1286
1288{
1289 double s = mu+b;
1290 double lls = - s;
1291 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1292 double llb = 0;
1293 if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;
1294
1295 return 2*( lls + llb);
1296}
1297
1298////////////////////////////////////////////////////////////////////////////////
1299/// Calculates the Profile Likelihood for MODEL 6:
1300/// Background known/Efficiency binomial
1301/// - what = 1: Maximum likelihood estimate is returned
1302/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1303/// - what = 3: Profile Likelihood of Test hypothesis is returned
1304/// otherwise parameters as described in the beginning of the class)
1305
1307{
1308 Double_t coef[4], roots[3];
1309 Double_t f = 0.;
1310 Double_t zm = Double_t(z) / m;
1311
1312 if (what == 1) {
1313 f = (x - b) / zm;
1314 }
1315
1316 if (what == 2) {
1317 mu = (x - b) / zm;
1318 Double_t e = zm;
1319 f = LikeMod6(mu, b, e, x, z, m);
1320 }
1321 if (what == 3) {
1322 Double_t e;
1323 if (mu == 0) {
1324 e = zm;
1325 } else {
1326 coef[3] = mu * mu;
1327 coef[2] = mu * b - mu * x - mu * mu - mu * m;
1328 coef[1] = mu * x - mu * b + mu * z - m * b;
1329 coef[0] = b * z;
1330 TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);
1331 e = roots[1];
1332 }
1333 f = LikeMod6(mu, b, e, x, z, m);
1334 }
1335 return f;
1336}
1337
1338////////////////////////////////////////////////////////////////////////////////
1339/// Profile Likelihood function for MODEL 6:
1340/// background known/ Efficiency binomial
1341
1343{
1344 double s = e*mu+b;
1345 double lls = - s;
1346 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1347
1348 double lle = 0;
1349 if (z == 0) lle = m * TMath::Log(1-e);
1350 else if ( z == m) lle = m * TMath::Log(e);
1351 else lle = z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z);
1352
1353 return 2* (lls + lle);
1354}
1355
1356
1357////////////////////////////////////////////////////////////////////////////////
1358/// Calculates the Profile Likelihood for MODEL 7:
1359/// background known/Efficiency Gauss
1360/// - what = 1: Maximum likelihood estimate is returned
1361/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1362/// - what = 3: Profile Likelihood of Test hypothesis is returned
1363/// otherwise parameters as described in the beginning of the class)
1364
1366{
1367 Double_t v = sde * sde;
1368 Double_t f = 0.;
1369
1370 if (what == 1) {
1371 f = (x - b) / em;
1372 }
1373
1374 if (what == 2) {
1375 mu = (x - b) / em;
1376 Double_t e = em;
1377 f = LikeMod7(mu, b, e, x, em, v);
1378 }
1379
1380 if (what == 3) {
1381 Double_t e;
1382 if (mu == 0) {
1383 e = em;
1384 } else {
1385 e = (-(mu * em - b - mu * mu * v) - TMath::Sqrt((mu * em - b - mu * mu * v) * (mu * em - b - mu * mu * v) + 4 * mu * (x * mu * v - mu * b * v + b * em))) / (- mu) / 2;
1386 }
1387 f = LikeMod7(mu, b, e, x, em, v);
1388 }
1389
1390 return f;
1391}
1392
1393////////////////////////////////////////////////////////////////////////////////
1394/// Profile Likelihood function for MODEL 6:
1395/// background known/ Efficiency gaussian
1396
1398{
1399 double s = e*mu+b;
1400 double lls = - s;
1401 if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1402
1403 double lle = 0;
1404 if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1405
1406 return 2*( lls + lle);
1407}
1408
1409////////////////////////////////////////////////////////////////////////////////
1410/// Evaluate polynomial
1411
1413{
1414 const Int_t *p;
1415 p = coef;
1416 Double_t ans = *p++;
1417 Int_t i = N;
1418
1419 do
1420 ans = ans * x + *p++;
1421 while (--i);
1422
1423 return ans;
1424}
1425
1426////////////////////////////////////////////////////////////////////////////////
1427/// Evaluate mononomial
1428
1430{
1431 Double_t ans;
1432 const Int_t *p;
1433
1434 p = coef;
1435 ans = x + *p++;
1436 Int_t i = N - 1;
1437
1438 do
1439 ans = ans * x + *p++;
1440 while (--i);
1441
1442 return ans;
1443}
1444
1445////////////////////////////////////////////////////////////////////////////////
1446/// LogFactorial function (use the logGamma function via the relation Gamma(n+1) = n!
1447
1449{
1450 return TMath::LnGamma(n+1);
1451}
SVector< double, 2 > v
Definition: Dict.h:5
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
#define N
double sqrt(double)
This class computes confidence intervals for the rate of a Poisson process in the presence of uncerta...
Definition: TRolke.h:34
Double_t f_sdb
Definition: TRolke.h:55
Double_t EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t sde, Double_t sdb, Int_t what)
Calculates the Profile Likelihood for MODEL 3: Gauss background/ Gauss Efficiency.
Definition: TRolke.cxx:1148
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
Double_t Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
Internal helper function 'Interval'.
Definition: TRolke.cxx:754
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
Double_t LikeMod7(Double_t mu, Double_t b, Double_t e, Int_t x, Double_t em, Double_t v)
Profile Likelihood function for MODEL 6: background known/ Efficiency gaussian.
Definition: TRolke.cxx:1397
Double_t EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t b, Int_t m, Int_t what)
Calculates the Profile Likelihood for MODEL 6: Background known/Efficiency binomial.
Definition: TRolke.cxx:1306
Double_t f_sde
Definition: TRolke.h:54
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
Double_t f_em
Definition: TRolke.h:51
static Double_t EvalPolynomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate polynomial.
Definition: TRolke.cxx:1412
Double_t f_b
Definition: TRolke.h:57
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
bool fBounding
Definition: TRolke.h:40
Double_t Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what)
Internal helper function Chooses between the different profile likelihood functions to use for the di...
Definition: TRolke.cxx:933
Int_t fNumWarningsDeprecated1
Definition: TRolke.h:42
Double_t LikeMod1(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Profile Likelihood function for MODEL 1: Poisson background/ Binomial Efficiency.
Definition: TRolke.cxx:1003
Double_t EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Int_t what)
Calculates the Profile Likelihood for MODEL 4: Poiss background/Efficiency known.
Definition: TRolke.cxx:1216
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
Double_t LikeMod5(Double_t mu, Double_t b, Int_t x, Double_t bm, Double_t u)
Profile Likelihood function for MODEL 5: Gauss background/Efficiency known.
Definition: TRolke.cxx:1287
Double_t GetUpperLimit()
Calculate and get upper limit for the pre-specified model.
Definition: TRolke.cxx:397
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
static Double_t EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate mononomial.
Definition: TRolke.cxx:1429
Double_t EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Int_t what)
Calculates the Profile Likelihood for MODEL 5: Gauss background/Efficiency known.
Definition: TRolke.cxx:1262
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
Double_t LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Gradient model likelihood.
Definition: TRolke.cxx:1062
TRolke(Double_t CL=0.9, Option_t *option="")
Constructor with optional Confidence Level argument.
Definition: TRolke.cxx:175
void ProfLikeMod1(Double_t mu, Double_t &b, Double_t &e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Helper for calculation of estimates of efficiency and background for model 1.
Definition: TRolke.cxx:1029
Double_t LikeMod3(Double_t mu, Double_t b, Double_t e, Int_t x, Double_t bm, Double_t em, Double_t u, Double_t v)
Profile Likelihood function for MODEL 3: Gauss background/Gauss Efficiency.
Definition: TRolke.cxx:1194
Double_t EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t sde, Double_t b, Int_t what)
Calculates the Profile Likelihood for MODEL 7: background known/Efficiency Gauss.
Definition: TRolke.cxx:1365
Int_t f_x
Definition: TRolke.h:47
Double_t fUpperLimit
Definition: TRolke.h:38
Int_t f_m
Definition: TRolke.h:58
void SetSwitch(bool bnd)
Deprecated name for SetBounding.
Definition: TRolke.cxx:579
Int_t f_z
Definition: TRolke.h:49
Double_t fCL
Definition: TRolke.h:37
void SetModelParameters()
Definition: TRolke.cxx:691
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
Double_t CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:637
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
Double_t LogFactorial(Int_t n)
LogFactorial function (use the logGamma function via the relation Gamma(n+1) = n!
Definition: TRolke.cxx:1448
Double_t ComputeInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
ComputeInterval, the internals.
Definition: TRolke.cxx:713
virtual ~TRolke()
Destructor.
Definition: TRolke.cxx:189
Int_t f_mid
Definition: TRolke.h:53
Double_t EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what)
Calculates the Profile Likelihood for MODEL 1: Poisson background/ Binomial Efficiency.
Definition: TRolke.cxx:967
Double_t GetLowerLimit()
Calculate and get lower limit for the pre-specified model.
Definition: TRolke.cxx:407
Int_t fNumWarningsDeprecated2
Definition: TRolke.h:43
Double_t LikeMod2(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Double_t em, Double_t tau, Double_t v)
Profile Likelihood function for MODEL 2: Poisson background/Gauss Efficiency.
Definition: TRolke.cxx:1126
Int_t f_y
Definition: TRolke.h:48
Double_t GetBackground()
Return a simple background value (estimate/truth) given the pre-specified model.
Definition: TRolke.cxx:417
Double_t f_tau
Definition: TRolke.h:56
Double_t fLowerLimit
Definition: TRolke.h:39
Double_t LikeMod6(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t z, Int_t m)
Profile Likelihood function for MODEL 6: background known/ Efficiency binomial.
Definition: TRolke.cxx:1342
Double_t EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t sde, Double_t tau, Int_t what)
Calculates the Profile Likelihood for MODEL 2: Poisson background/ Gauss Efficiency.
Definition: TRolke.cxx:1081
Double_t f_e
Definition: TRolke.h:52
Double_t LikeMod4(Double_t mu, Double_t b, Int_t x, Int_t y, Double_t tau)
Profile Likelihood function for MODEL 4: Poiss background/Efficiency known.
Definition: TRolke.cxx:1242
void Print(Option_t *) const
Dump internals. Print members.
Definition: TRolke.cxx:593
Double_t f_bm
Definition: TRolke.h:50
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
static constexpr double s
static constexpr double mu0
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
Double_t PoissonI(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par) This is a non-smooth function.
Definition: TMath.cxx:599
Double_t Log(Double_t x)
Definition: TMath.h:748
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:486
Bool_t RootsCubic(const Double_t coef[4], Double_t &a, Double_t &b, Double_t &c)
Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where.
Definition: TMath.cxx:1091
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition: TMath.cxx:2164
Definition: test.py:1
TCanvas * roots()
Definition: roots.C:1
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4