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