ROOT   6.14/05 Reference Guide
TRolke.cxx
Go to the documentation of this file.
1 // @(#)root/physics:$Id$
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
26 1: 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
35 2: 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
44 3: 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
52 4: 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
60 5: SetGaussBkgKnownEff(x,bm,sdb,e)
61 ~~~
62  Background: Gaussian
63  Efficiency: Known
64 ~~~
65  when background is Gaussian
66
67 6: SetKnownBkgBinomEff(x,z,b,m)
68 ~~~
69  Background: Known
70  Efficiency: Binomial
71 ~~~
72  when signal efficiency was determined from Monte Carlo
73
74 7: 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
446 bool TRolke::GetSensitivity(Double_t& low, Double_t& high, Double_t pPrecision)
447 {
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
481 bool TRolke::GetLimitsQuantile(Double_t& low, Double_t& high, Int_t& out_x, Double_t integral)
482 {
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
511 bool TRolke::GetLimitsML(Double_t& low, Double_t& high, Int_t& out_x)
512 {
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 {
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;
558  ComputeInterval(rolke_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
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
579 void 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
593 void TRolke::Print(Option_t*) const {
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
638  if (fNumWarningsDeprecated2<2 ) {
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
913 done:
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
1023 struct 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 ////////////////////////////////////////////////////////////////////////////////
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 }
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 &#39;Interval&#39;.
Definition: TRolke.cxx:754
Double_t PoissonI(Double_t x, Double_t par)
TRolke(Double_t CL=0.9, Option_t *option="")
Constructor with optional Confidence Level argument.
Definition: TRolke.cxx:175
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
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t mid
Definition: TRolke.cxx:630
Double_t Log(Double_t x)
Definition: TMath.h:759
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 f_em
Definition: TRolke.h:51
const char Option_t
Definition: RtypesCore.h:62
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
Int_t f_y
Definition: TRolke.h:48
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)
Double_t f_sdb
Definition: TRolke.h:55
#define N
Definition: test.py:1
void SetModelParameters()
Definition: TRolke.cxx:691
Int_t f_x
Definition: TRolke.h:47
Int_t fNumWarningsDeprecated2
Definition: TRolke.h:43
#define f(i)
Definition: RSha256.hxx:104
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
int Int_t
Definition: RtypesCore.h:41
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
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
you should not use this method at all Int_t y
Definition: TRolke.cxx:630
static constexpr double mu0
Int_t f_z
Definition: TRolke.h:49
double sqrt(double)
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
static Double_t EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate mononomial.
Definition: TRolke.cxx:1429
Double_t f_bm
Definition: TRolke.h:50
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
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
you should not use this method at all Int_t Int_t Double_t Double_t em
Definition: TRolke.cxx:630
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 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
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
Double_t f_b
Definition: TRolke.h:57
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 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
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 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
you should not use this method at all Int_t Int_t Double_t bm
Definition: TRolke.cxx:630
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
you should not use this method at all * sde
Definition: TRolke.cxx:630
TCanvas * roots()
Definition: roots.C:1
Double_t f_sde
Definition: TRolke.h:54
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t Int_t m
Definition: TRolke.cxx:637
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
void Print(Option_t *) const
Dump internals. Print members.
Definition: TRolke.cxx:593
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
SVector< double, 2 > v
Definition: Dict.h:5
void SetSwitch(bool bnd)
Deprecated name for SetBounding.
Definition: TRolke.cxx:579
Double_t fUpperLimit
Definition: TRolke.h:38
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 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 f_tau
Definition: TRolke.h:56
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
virtual ~TRolke()
Destructor.
Definition: TRolke.cxx:189
Int_t f_m
Definition: TRolke.h:58
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
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 LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Definition: TRolke.cxx:1062
Double_t GetLowerLimit()
Calculate and get lower limit for the pre-specified model.
Definition: TRolke.cxx:407
#define ClassImp(name)
Definition: Rtypes.h:359
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
double Double_t
Definition: RtypesCore.h:55
Double_t GetUpperLimit()
Calculate and get upper limit for the pre-specified model.
Definition: TRolke.cxx:397
Double_t f_e
Definition: TRolke.h:52
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
static constexpr double s
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
void SetBounding(const bool bnd)
Definition: TRolke.h:184
Int_t f_mid
Definition: TRolke.h:53
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
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
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
static Double_t EvalPolynomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate polynomial.
Definition: TRolke.cxx:1412
Double_t GetBackground()
Return a simple background value (estimate/truth) given the pre-specified model.
Definition: TRolke.cxx:417
auto * l
Definition: textangle.C:4
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Double_t fCL
Definition: TRolke.h:37
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:486
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
This class computes confidence intervals for the rate of a Poisson process in the presence of uncerta...
Definition: TRolke.h:33
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
Bool_t RootsCubic(const Double_t coef[4], Double_t &a, Double_t &b, Double_t &c)
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
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
Int_t fNumWarningsDeprecated1
Definition: TRolke.h:42
Double_t fLowerLimit
Definition: TRolke.h:39
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t sdb
Definition: TRolke.cxx:630
const Int_t n
Definition: legend1.C:16
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
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