Logo ROOT   6.12/07
Reference Guide
TRolke.cxx
Go to the documentation of this file.
1 // @(#)root/physics:$Id$
2 // Author: Jan Conrad
3 
4 /** \class TRolke
5  \ingroup Physics
6  This class computes confidence intervals for the rate of a Poisson
7  process in the presence of uncertain background and/or efficiency.
8 
9  The treatment and the resulting limits are fully frequentist. The
10  limit calculations make use of the profile likelihood method.
11 
12 \author Jan Conrad (CERN) 2004, Updated: Johan Lundberg (CERN) 2009
13 
14  For a full list of methods and their syntax, and build instructions,
15  consult the header file TRolke.h.
16 
17  Examples/tutorials are found in the separate file Rolke.C
18 
19 ### TRolke implements the following Models
20 
21  The signal is always assumed to be Poisson, with the following
22  combinations of models of background and detection efficiency:
23 
24  If unsure, first consider model 3, 4 or 5.
25 
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 ////////////////////////////////////////////////////////////////////////////////
1060 /// Gradient model likelihood
1061 
1063 {
1064  Double_t eta, etaprime, bprime, f;
1065  eta = static_cast<double>(z) / e - static_cast<double>(m - z) / (1.0 - e);
1066  etaprime = (-1) * (static_cast<double>(m - z) / ((1.0 - e) * (1.0 - e)) + static_cast<double>(z) / (e * e));
1067  Double_t b = y / (tau - eta / mu);
1068  bprime = (b * b * etaprime) / mu / y;
1069  f = (mu + bprime) * (x / (e * mu + b) - 1) + (y / b - tau) * bprime + eta;
1070  return f;
1071 }
1072 
1073 ////////////////////////////////////////////////////////////////////////////////
1074 /// Calculates the Profile Likelihood for MODEL 2:
1075 /// Poisson background/ Gauss Efficiency
1076 /// - what = 1: Maximum likelihood estimate is returned
1077 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1078 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1079 /// otherwise parameters as described in the beginning of the class)
1080 
1082 {
1083  Double_t v = sde * sde;
1084  Double_t coef[4], roots[3];
1085  Double_t f = 0;
1086 
1087  if (what == 1) {
1088  f = (x - y / tau) / em;
1089  }
1090 
1091  if (what == 2) {
1092  mu = (x - y / tau) / em;
1093  Double_t b = y / tau;
1094  Double_t e = em;
1095  f = LikeMod2(mu, b, e, x, y, em, tau, v);
1096  }
1097 
1098  if (what == 3) {
1099  if (mu == 0) {
1100  Double_t b = (x + y) / (1 + tau);
1101  Double_t e = em ;
1102  f = LikeMod2(mu, b, e, x, y, em, tau, v);
1103  } else {
1104  coef[3] = mu;
1105  coef[2] = mu * mu * v - 2 * em * mu - mu * mu * v * tau;
1106  coef[1] = (- x) * mu * v - mu * mu * mu * v * v * tau - mu * mu * v * em + em * mu * mu * v * tau + em * em * mu - y * mu * v;
1107  coef[0] = x * mu * mu * v * v * tau + x * em * mu * v - y * mu * mu * v * v + y * em * mu * v;
1108 
1109  TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);
1110 
1111  Double_t e = roots[1];
1112  Double_t b;
1113  if ( v > 0) b = y / (tau + (em - e) / mu / v);
1114  else b = y/tau;
1115  f = LikeMod2(mu, b, e, x, y, em, tau, v);
1116  }
1117  }
1118 
1119  return f;
1120 }
1121 
1122 ////////////////////////////////////////////////////////////////////////////////
1123 /// Profile Likelihood function for MODEL 2:
1124 /// Poisson background/Gauss Efficiency
1125 
1127 {
1128  double s = e*mu+b;
1129  double lls = - s;
1130  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1131  double bg = tau*b;
1132  double llb = -bg;
1133  if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y);
1134  double lle = 0;
1135  if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1136 
1137  return 2*( lls + llb + lle);
1138 }
1139 
1140 ////////////////////////////////////////////////////////////////////////////////
1141 /// Calculates the Profile Likelihood for MODEL 3:
1142 /// Gauss background/ Gauss Efficiency
1143 /// - what = 1: Maximum likelihood estimate is returned
1144 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1145 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1146 /// otherwise parameters as described in the beginning of the class)
1147 
1149 {
1150  Double_t f = 0.;
1151  Double_t v = sde * sde;
1152  Double_t u = sdb * sdb;
1153 
1154  if (what == 1) {
1155  f = (x - bm) / em;
1156  }
1157 
1158 
1159  if (what == 2) {
1160  mu = (x - bm) / em;
1161  Double_t b = bm;
1162  Double_t e = em;
1163  f = LikeMod3(mu, b, e, x, bm, em, u, v);
1164  }
1165 
1166 
1167  if (what == 3) {
1168  if (mu == 0.0) {
1169  Double_t b = ((bm - u) + TMath::Sqrt((bm - u) * (bm - u) + 4 * x * u)) / 2.;
1170  Double_t e = em;
1171  f = LikeMod3(mu, b, e, x, bm, em, u, v);
1172  } else {
1173  Double_t e = em;
1174  Double_t b = bm;
1175  if ( v > 0) {
1176  Double_t temp[3];
1177  temp[0] = mu * mu * v + u;
1178  temp[1] = mu * mu * mu * v * v + mu * v * u - mu * mu * v * em + mu * v * bm - 2 * u * em;
1179  temp[2] = mu * mu * v * v * bm - mu * v * u * em - mu * v * bm * em + u * em * em - mu * mu * v * v * x;
1180  e = (-temp[1] + TMath::Sqrt(temp[1] * temp[1] - 4 * temp[0] * temp[2])) / 2 / temp[0];
1181  b = bm - (u * (em - e)) / v / mu;
1182  }
1183  f = LikeMod3(mu, b, e, x, bm, em, u, v);
1184  }
1185  }
1186 
1187  return f;
1188 }
1189 
1190 ////////////////////////////////////////////////////////////////////////////////
1191 /// Profile Likelihood function for MODEL 3:
1192 /// Gauss background/Gauss Efficiency
1193 
1195 {
1196  double s = e*mu+b;
1197  double lls = - s;
1198  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1199  double llb = 0;
1200  if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;
1201  double lle = 0;
1202  if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1203 
1204  return 2*( lls + llb + lle);
1205 
1206 }
1207 
1208 ////////////////////////////////////////////////////////////////////////////////
1209 /// Calculates the Profile Likelihood for MODEL 4:
1210 /// Poiss background/Efficiency known
1211 /// - what = 1: Maximum likelihood estimate is returned
1212 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1213 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1214 /// otherwise parameters as described in the beginning of the class)
1215 
1217 {
1218  Double_t f = 0.0;
1219 
1220  if (what == 1) f = x - y / tau;
1221  if (what == 2) {
1222  mu = x - y / tau;
1223  Double_t b = y / tau;
1224  f = LikeMod4(mu, b, x, y, tau);
1225  }
1226  if (what == 3) {
1227  if (mu == 0.0) {
1228  Double_t b = Double_t(x + y) / (1 + tau);
1229  f = LikeMod4(mu, b, x, y, tau);
1230  } else {
1231  Double_t b = (x + y - (1 + tau) * mu + sqrt((x + y - (1 + tau) * mu) * (x + y - (1 + tau) * mu) + 4 * (1 + tau) * y * mu)) / 2 / (1 + tau);
1232  f = LikeMod4(mu, b, x, y, tau);
1233  }
1234  }
1235  return f;
1236 }
1237 
1238 ////////////////////////////////////////////////////////////////////////////////
1239 /// Profile Likelihood function for MODEL 4:
1240 /// Poiss background/Efficiency known
1241 
1243 {
1244  double s = mu+b;
1245  double lls = - s;
1246  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1247  double bg = tau*b;
1248  double llb = -bg;
1249  if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y);
1250 
1251  return 2*( lls + llb);
1252 }
1253 
1254 ////////////////////////////////////////////////////////////////////////////////
1255 /// Calculates the Profile Likelihood for MODEL 5:
1256 /// Gauss background/Efficiency known
1257 /// - what = 1: Maximum likelihood estimate is returned
1258 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1259 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1260 /// otherwise parameters as described in the beginning of the class)
1261 
1263 {
1264  Double_t u = sdb * sdb;
1265  Double_t f = 0;
1266 
1267  if (what == 1) {
1268  f = x - bm;
1269  }
1270  if (what == 2) {
1271  mu = x - bm;
1272  Double_t b = bm;
1273  f = LikeMod5(mu, b, x, bm, u);
1274  }
1275 
1276  if (what == 3) {
1277  Double_t b = ((bm - u - mu) + TMath::Sqrt((bm - u - mu) * (bm - u - mu) - 4 * (mu * u - mu * bm - u * x))) / 2;
1278  f = LikeMod5(mu, b, x, bm, u);
1279  }
1280  return f;
1281 }
1282 
1283 ////////////////////////////////////////////////////////////////////////////////
1284 /// Profile Likelihood function for MODEL 5:
1285 /// Gauss background/Efficiency known
1286 
1288 {
1289  double s = mu+b;
1290  double lls = - s;
1291  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1292  double llb = 0;
1293  if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;
1294 
1295  return 2*( lls + llb);
1296 }
1297 
1298 ////////////////////////////////////////////////////////////////////////////////
1299 /// Calculates the Profile Likelihood for MODEL 6:
1300 /// Background known/Efficiency binomial
1301 /// - what = 1: Maximum likelihood estimate is returned
1302 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1303 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1304 /// otherwise parameters as described in the beginning of the class)
1305 
1307 {
1308  Double_t coef[4], roots[3];
1309  Double_t f = 0.;
1310  Double_t zm = Double_t(z) / m;
1311 
1312  if (what == 1) {
1313  f = (x - b) / zm;
1314  }
1315 
1316  if (what == 2) {
1317  mu = (x - b) / zm;
1318  Double_t e = zm;
1319  f = LikeMod6(mu, b, e, x, z, m);
1320  }
1321  if (what == 3) {
1322  Double_t e;
1323  if (mu == 0) {
1324  e = zm;
1325  } else {
1326  coef[3] = mu * mu;
1327  coef[2] = mu * b - mu * x - mu * mu - mu * m;
1328  coef[1] = mu * x - mu * b + mu * z - m * b;
1329  coef[0] = b * z;
1330  TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);
1331  e = roots[1];
1332  }
1333  f = LikeMod6(mu, b, e, x, z, m);
1334  }
1335  return f;
1336 }
1337 
1338 ////////////////////////////////////////////////////////////////////////////////
1339 /// Profile Likelihood function for MODEL 6:
1340 /// background known/ Efficiency binomial
1341 
1343 {
1344  double s = e*mu+b;
1345  double lls = - s;
1346  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1347 
1348  double lle = 0;
1349  if (z == 0) lle = m * TMath::Log(1-e);
1350  else if ( z == m) lle = m * TMath::Log(e);
1351  else lle = z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z);
1352 
1353  return 2* (lls + lle);
1354 }
1355 
1356 
1357 ////////////////////////////////////////////////////////////////////////////////
1358 /// Calculates the Profile Likelihood for MODEL 7:
1359 /// background known/Efficiency Gauss
1360 /// - what = 1: Maximum likelihood estimate is returned
1361 /// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
1362 /// - what = 3: Profile Likelihood of Test hypothesis is returned
1363 /// otherwise parameters as described in the beginning of the class)
1364 
1366 {
1367  Double_t v = sde * sde;
1368  Double_t f = 0.;
1369 
1370  if (what == 1) {
1371  f = (x - b) / em;
1372  }
1373 
1374  if (what == 2) {
1375  mu = (x - b) / em;
1376  Double_t e = em;
1377  f = LikeMod7(mu, b, e, x, em, v);
1378  }
1379 
1380  if (what == 3) {
1381  Double_t e;
1382  if (mu == 0) {
1383  e = em;
1384  } else {
1385  e = (-(mu * em - b - mu * mu * v) - TMath::Sqrt((mu * em - b - mu * mu * v) * (mu * em - b - mu * mu * v) + 4 * mu * (x * mu * v - mu * b * v + b * em))) / (- mu) / 2;
1386  }
1387  f = LikeMod7(mu, b, e, x, em, v);
1388  }
1389 
1390  return f;
1391 }
1392 
1393 ////////////////////////////////////////////////////////////////////////////////
1394 /// Profile Likelihood function for MODEL 6:
1395 /// background known/ Efficiency gaussian
1396 
1398 {
1399  double s = e*mu+b;
1400  double lls = - s;
1401  if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
1402 
1403  double lle = 0;
1404  if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;
1405 
1406  return 2*( lls + lle);
1407 }
1408 
1409 ////////////////////////////////////////////////////////////////////////////////
1410 /// Evaluate polynomial
1411 
1413 {
1414  const Int_t *p;
1415  p = coef;
1416  Double_t ans = *p++;
1417  Int_t i = N;
1418 
1419  do
1420  ans = ans * x + *p++;
1421  while (--i);
1422 
1423  return ans;
1424 }
1425 
1426 ////////////////////////////////////////////////////////////////////////////////
1427 /// Evaluate mononomial
1428 
1430 {
1431  Double_t ans;
1432  const Int_t *p;
1433 
1434  p = coef;
1435  ans = x + *p++;
1436  Int_t i = N - 1;
1437 
1438  do
1439  ans = ans * x + *p++;
1440  while (--i);
1441 
1442  return ans;
1443 }
1444 
1445 ////////////////////////////////////////////////////////////////////////////////
1446 /// LogFactorial function (use the logGamma function via the relation Gamma(n+1) = n!
1447 
1449 {
1450  return TMath::LnGamma(n+1);
1451 }
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)
compute the Poisson distribution function for (x,par) This is a non-smooth function.
Definition: TMath.cxx:602
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:648
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
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)
Evaluate the quantiles of the chi-squared probability distribution function.
Definition: TMath.cxx:2148
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)
Gradient model likelihood.
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:489
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)
Calculates roots of polynomial of 3rd order a*x^3 + b*x^2 + c*x + d, where a == coef[3], b == coef[2], c == coef[1], d == coef[0] coef[3] must be different from 0 If the boolean returned by the method is false: ==> there are 3 real roots a,b,c If the boolean returned by the method is true: ==> there is one real root a and 2 complex conjugates roots (b+i*c,b-i*c)
Definition: TMath.cxx:1082
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
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