source: branches/Mars_McMismatchStudy/mtools/MRolke.cc@ 18679

Last change on this file since 18679 was 8105, checked in by meyer, 18 years ago
*** empty log message ***
File size: 24.0 KB
Line 
1// @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1.2 2006-10-17 16:32:52 meyer Exp $
2// Author: Jan Conrad 9/2/2004
3
4/*************************************************************************
5 * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12//////////////////////////////////////////////////////////////////////////////
13//
14// MRolke
15//
16// This class computes confidence intervals for the rate of a Poisson
17// in the presence of background and efficiency with a fully frequentist
18// treatment of the uncertainties in the efficiency and background estimate
19// using the profile likelihood method.
20//
21// The signal is always assumed to be Poisson.
22//
23// The method is very similar to the one used in MINUIT (MINOS).
24//
25// Two options are offered to deal with cases where the maximum likelihood
26// estimate (MLE) is not in the physical region. Version "bounded likelihood"
27// is the one used by MINOS if bounds for the physical region are chosen. Versi// on "unbounded likelihood (the default) allows the MLE to be in the
28// unphysical region. It has however better coverage.
29// For more details consult the reference (see below).
30//
31//
32// It allows the following Models:
33//
34// 1: Background - Poisson, Efficiency - Binomial (cl,x,y,z,tau,m)
35// 2: Background - Poisson, Efficiency - Gaussian (cl,xd,y,em,tau,sde)
36// 3: Background - Gaussian, Efficiency - Gaussian (cl,x,bm,em,sd)
37// 4: Background - Poisson, Efficiency - known (cl,x,y,tau,e)
38// 5: Background - Gaussian, Efficiency - known (cl,x,y,z,sdb,e)
39// 6: Background - known, Efficiency - Binomial (cl,x,z,m,b)
40// 7: Background - known, Efficiency - Gaussian (cl,x,em,sde,b)
41//
42// Parameter definition:
43//
44// cl = Confidence level
45//
46// x = number of observed events
47//
48// y = number of background events
49//
50// z = number of simulated signal events
51//
52// em = measurement of the efficiency.
53//
54// bm = background estimate
55//
56// tau = ratio between signal and background region (in case background is
57// observed) ratio between observed and simulated livetime in case
58// background is determined from MC.
59//
60// sd(x) = sigma of the Gaussian
61//
62// e = true efficiency (in case known)
63//
64// b = expected background (in case known)
65//
66// m = number of MC runs
67//
68// mid = ID number of the model ...
69//
70// For a description of the method and its properties:
71//
72// W.Rolke, A. Lopez, J. Conrad and Fred James
73// "Limits and Confidence Intervals in presence of nuisance parameters"
74// http://lanl.arxiv.org/abs/physics/0403059
75// Nucl.Instrum.Meth.A551:493-503,2005
76//
77// Should I use MRolke, TFeldmanCousins, TLimit?
78// ============================================
79// 1. I guess MRolke makes TFeldmanCousins obsolete?
80//
81// Certainly not. TFeldmanCousins is the fully frequentist construction and
82// should be used in case of no (or negligible uncertainties). It is however
83// not capable of treating uncertainties in nuisance parameters.
84// MRolke is desined for this case and it is shown in the reference above
85// that it has good coverage properties for most cases, ie it might be
86// used where FeldmannCousins can't.
87//
88// 2. What are the advantages of MRolke over TLimit?
89//
90// MRolke is fully frequentist. TLimit treats nuisance parameters Bayesian.
91// For a coverage study of a Bayesian method refer to
92// physics/0408039 (Tegenfeldt & J.C). However, this note studies
93// the coverage of Feldman&Cousins with Bayesian treatment of nuisance
94// parameters. To make a long story short: using the Bayesian method you
95// might introduce a small amount of over-coverage (though I haven't shown it
96// for TLimit). On the other hand, coverage of course is a not so interesting
97// when you consider yourself a Bayesian.
98//
99// Author: Jan Conrad (CERN)
100//
101// see example in tutorial Rolke.C
102//
103// Copyright CERN 2004 Jan.Conrad@cern.ch
104//
105///////////////////////////////////////////////////////////////////////////
106
107
108#include "MRolke.h"
109#include "TMath.h"
110#include "Riostream.h"
111
112ClassImp(MRolke)
113
114//__________________________________________________________________________
115MRolke::MRolke(Double_t CL, Option_t * /*option*/)
116{
117 //constructor
118 fUpperLimit = 0.0;
119 fLowerLimit = 0.0;
120 fCL = CL;
121 fSwitch = 0; // 0: unbounded likelihood
122 // 1: bounded likelihood
123}
124
125//___________________________________________________________________________
126MRolke::~MRolke()
127{
128}
129
130
131//___________________________________________________________________________
132Double_t MRolke::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)
133{
134 //calculate interval
135 Int_t done = 0;
136 Double_t limit[2];
137
138 limit[1] = Interval(x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
139
140 if (limit[1] > 0) {
141 done = 1;
142 }
143
144 if (fSwitch == 0) {
145
146 Int_t trial_x = x;
147
148 while (done == 0) {
149 trial_x++;
150 limit[1] = Interval(trial_x,y,z,bm,em,e,mid, sde,sdb,tau,b,m);
151 if (limit[1] > 0) done = 1;
152 }
153 }
154
155 return limit[1];
156}
157
158
159
160
161//_____________________________________________________________________
162Double_t MRolke::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)
163{
164 // Calculates the Confidence Interval
165
166 //Double_t dchi2 = Chi2Percentile(1,1-fCL);
167 Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1);
168
169 Double_t tempxy[2],limits[2] = {0,0};
170 Double_t slope,fmid,low,flow,high,fhigh,test,ftest,mu0,maximum,target,l,f0;
171 Double_t med = 0;
172 Double_t maxiter=1000, acc = 0.00001;
173 Int_t i;
174 Int_t bp = 0;
175
176 if ((mid != 3) && (mid != 5)) bm = (Double_t)y;
177
178 if ((mid == 3) || (mid == 5)) {
179 if (bm == 0) bm = 0.00001;
180 }
181
182 if ((mid <= 2) || (mid == 4)) bp = 1;
183
184
185 if (bp == 1 && x == 0 && bm > 0 ){
186
187 for(Int_t i = 0; i < 2; i++) {
188 x++;
189 tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
190 }
191 slope = tempxy[1] - tempxy[0];
192 limits[1] = tempxy[0] - slope;
193 limits[0] = 0.0;
194 if (limits[1] < 0) limits[1] = 0.0;
195 goto done;
196 }
197
198 if (bp != 1 && x == 0){
199
200 for(Int_t i = 0; i < 2; i++) {
201 x++;
202 tempxy[i] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
203 }
204 slope = tempxy[1] - tempxy[0];
205 limits[1] = tempxy[0] - slope;
206 limits[0] = 0.0;
207 if (limits[1] < 0) limits[1] = 0.0;
208 goto done;
209 }
210
211 if (bp != 1 && bm == 0){
212 for(Int_t i = 0; i < 2; i++) {
213 bm++;
214 limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
215 tempxy[i] = limits[1];
216 }
217 slope = tempxy[1] - tempxy[0];
218 limits[1] = tempxy[0] - slope;
219 if (limits[1] < 0) limits[1] = 0;
220 goto done;
221 }
222
223
224 if (x == 0 && bm == 0){
225 x++;
226 bm++;
227
228 limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
229 tempxy[0] = limits[1];
230 x = 1;
231 bm = 2;
232 limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
233 tempxy[1] = limits[1];
234 x = 2;
235 bm = 1;
236 limits[1] = Interval(x,y,z,bm,em,e,mid,sde,sdb,tau,b,m);
237 limits[1] = 3*tempxy[0] -tempxy[1] - limits[1];
238 if (limits[1] < 0) limits[1] = 0;
239 goto done;
240 }
241
242 mu0 = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,1);
243
244 maximum = Likelihood(0,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,2);
245
246 test = 0;
247
248 f0 = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
249
250 if ( fSwitch == 1 ) { // do this only for the unbounded likelihood case
251 if ( mu0 < 0 ) maximum = f0;
252 }
253
254 target = maximum - dchi2;
255
256 if (f0 > target) {
257 limits[0] = 0;
258 } else {
259 if (mu0 < 0){
260 limits[0] = 0;
261 limits[1] = 0;
262 }
263
264 low = 0;
265 flow = f0;
266 high = mu0;
267 fhigh = maximum;
268
269 for(Int_t i = 0; i < maxiter; i++) {
270 l = (target-fhigh)/(flow-fhigh);
271 if (l < 0.2) l = 0.2;
272 if (l > 0.8) l = 0.8;
273
274 med = l*low + (1-l)*high;
275 if(med < 0.01){
276 limits[1]=0.0;
277 goto done;
278 }
279
280 fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
281
282 if (fmid > target) {
283 high = med;
284 fhigh = fmid;
285 } else {
286 low = med;
287 flow = fmid;
288 }
289 if ((high-low) < acc*high) break;
290 }
291 limits[0] = med;
292 }
293
294
295 if(mu0 > 0) {
296 low = mu0;
297 flow = maximum;
298 } else {
299 low = 0;
300 flow = f0;
301 }
302
303 test = low +1 ;
304
305 ftest = Likelihood(test,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
306
307 if (ftest < target) {
308 high = test;
309 fhigh = ftest;
310 } else {
311 slope = (ftest - flow)/(test - low);
312 high = test + (target -ftest)/slope;
313 fhigh = Likelihood(high,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
314 }
315
316 for(i = 0; i < maxiter; i++) {
317 l = (target-fhigh)/(flow-fhigh);
318 if (l < 0.2) l = 0.2;
319 if (l > 0.8) l = 0.8;
320 med = l * low + (1.-l)*high;
321 fmid = Likelihood(med,x,y,z,bm,em,e,mid,sde,sdb,tau,b,m,3);
322
323 if (fmid < target) {
324 high = med;
325 fhigh = fmid;
326 } else {
327 low = med;
328 flow = fmid;
329 }
330 if (high-low < acc*high) break;
331 }
332 limits[1] = med;
333
334done:
335
336 if ( (mid == 4) || (mid==5) ) {
337 limits[0] /= e;
338 limits[1] /= e;
339 }
340
341
342 fUpperLimit = limits[1];
343 fLowerLimit = TMath::Max(limits[0],0.0);
344
345
346 return limits[1];
347}
348
349
350//___________________________________________________________________________
351Double_t MRolke::Likelihood(Double_t mu, 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, Int_t what)
352{
353 // Chooses between the different profile likelihood functions to use for the
354 // different models.
355 // Returns evaluation of the profile likelihood functions.
356
357 switch (mid) {
358 case 1: return EvalLikeMod1(mu,x,y,z,e,tau,b,m,what);
359 case 2: return EvalLikeMod2(mu,x,y,em,e,sde,tau,b,what);
360 case 3: return EvalLikeMod3(mu,x,bm,em,e,sde,sdb,b,what);
361 case 4: return EvalLikeMod4(mu,x,y,tau,b,what);
362 case 5: return EvalLikeMod5(mu,x,bm,sdb,b,what);
363 case 6: return EvalLikeMod6(mu,x,z,e,b,m,what);
364 case 7: return EvalLikeMod7(mu,x,em,e,sde,b,what);
365 }
366
367 return 0;
368}
369
370//_________________________________________________________________________
371Double_t MRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t e, Double_t tau, Double_t b, Int_t m, Int_t what)
372{
373 // Calculates the Profile Likelihood for MODEL 1:
374 // Poisson background/ Binomial Efficiency
375 // what = 1: Maximum likelihood estimate is returned
376 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
377 // what = 3: Profile Likelihood of Test hypothesis is returned
378 // otherwise parameters as described in the beginning of the class)
379
380 Double_t f = 0;
381 Double_t zm = Double_t(z)/m;
382
383 if (what == 1) {
384 f = (x-y/tau)/zm;
385 }
386
387 if (what == 2) {
388 mu = (x-y/tau)/zm;
389 b = y/tau;
390 Double_t e = zm;
391 f = LikeMod1(mu,b,e,x,y,z,tau,m);
392 }
393
394 if (what == 3) {
395 if (mu == 0){
396 b = (x+y)/(1.0+tau);
397 e = zm;
398 f = LikeMod1(mu,b,e,x,y,z,tau,m);
399 } else {
400 MRolke g;
401 g.ProfLikeMod1(mu,b,e,x,y,z,tau,m);
402 f = LikeMod1(mu,b,e,x,y,z,tau,m);
403 }
404 }
405
406 return f;
407}
408
409//________________________________________________________________________
410Double_t MRolke::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)
411{
412 // Profile Likelihood function for MODEL 1:
413 // Poisson background/ Binomial Efficiency
414
415 return 2*(x*TMath::Log(e*mu+b)-(e*mu +b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) + TMath::Log(TMath::Factorial(m)) - TMath::Log(TMath::Factorial(m-z)) - TMath::Log(TMath::Factorial(z))+ z * TMath::Log(e) + (m-z)*TMath::Log(1-e));
416}
417
418//________________________________________________________________________
419void MRolke::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)
420{
421 // Void needed to calculate estimates of efficiency and background for model 1
422
423 Double_t med = 0.0,fmid;
424 Int_t maxiter =1000;
425 Double_t acc = 0.00001;
426 Double_t emin = ((m+mu*tau)-TMath::Sqrt((m+mu*tau)*(m+mu*tau)-4 * mu* tau * z))/2/mu/tau;
427
428 Double_t low = TMath::Max(1e-10,emin+1e-10);
429 Double_t high = 1 - 1e-10;
430
431 for(Int_t i = 0; i < maxiter; i++) {
432 med = (low+high)/2.;
433
434 fmid = LikeGradMod1(med,mu,x,y,z,tau,m);
435
436 if(high < 0.5) acc = 0.00001*high;
437 else acc = 0.00001*(1-high);
438
439 if ((high - low) < acc*high) break;
440
441 if(fmid > 0) low = med;
442 else high = med;
443 }
444
445 e = med;
446 Double_t eta = Double_t(z)/e -Double_t(m-z)/(1-e);
447
448 b = Double_t(y)/(tau -eta/mu);
449}
450
451//___________________________________________________________________________
452Double_t MRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x,Int_t y,Int_t z,Double_t tau,Int_t m)
453{
454 //gradient model
455 Double_t eta, etaprime, bprime,f;
456 eta = static_cast<double>(z)/e - static_cast<double>(m-z)/(1.0 - e);
457 etaprime = (-1) * (static_cast<double>(m-z)/((1.0 - e)*(1.0 - e)) + static_cast<double>(z)/(e*e));
458 Double_t b = y/(tau - eta/mu);
459 bprime = (b*b * etaprime)/mu/y;
460 f = (mu + bprime) * (x/(e * mu + b) - 1)+(y/b - tau) * bprime + eta;
461 return f;
462}
463
464//___________________________________________________________________________
465Double_t MRolke::EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t e,Double_t sde, Double_t tau, Double_t b, Int_t what)
466{
467 // Calculates the Profile Likelihood for MODEL 2:
468 // Poisson background/ Gauss Efficiency
469 // what = 1: Maximum likelihood estimate is returned
470 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
471 // what = 3: Profile Likelihood of Test hypothesis is returned
472 // otherwise parameters as described in the beginning of the class)
473
474 Double_t v = sde*sde;
475 Double_t coef[4],roots[3];
476 Double_t f = 0;
477
478 if (what == 1) {
479 f = (x-y/tau)/em;
480 }
481
482 if (what == 2) {
483 mu = (x-y/tau)/em;
484 b = y/tau;
485 e = em;
486 f = LikeMod2(mu,b,e,x,y,em,tau,v);
487 }
488
489 if (what == 3) {
490 if (mu == 0 ) {
491 b = (x+y)/(1+tau);
492 f = LikeMod2(mu,b,e,x,y,em,tau,v);
493 } else {
494 coef[3] = mu;
495 coef[2] = mu*mu*v-2*em*mu-mu*mu*v*tau;
496 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;
497 coef[0] = x*mu*mu*v*v*tau + x*em*mu*v - y*mu*mu*v*v + y*em*mu*v;
498
499 TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
500
501 e = roots[1];
502 b = y/(tau + (em - e)/mu/v);
503 f = LikeMod2(mu,b,e,x,y,em,tau,v);
504 }
505 }
506
507 return f;
508}
509
510//_________________________________________________________________________
511Double_t MRolke::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)
512{
513 // Profile Likelihood function for MODEL 2:
514 // Poisson background/Gauss Efficiency
515
516 return 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
517}
518
519//_____________________________________________________________________
520
521Double_t MRolke::EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t e, Double_t sde, Double_t sdb, Double_t b, Int_t what)
522{
523 // Calculates the Profile Likelihood for MODEL 3:
524 // Gauss background/ Gauss Efficiency
525 // what = 1: Maximum likelihood estimate is returned
526 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
527 // what = 3: Profile Likelihood of Test hypothesis is returned
528 // otherwise parameters as described in the beginning of the class)
529
530 Double_t f = 0.;
531 Double_t v = sde*sde;
532 Double_t u = sdb*sdb;
533
534 if (what == 1) {
535 f = (x-bm)/em;
536 }
537
538
539 if (what == 2) {
540 mu = (x-bm)/em;
541 b = bm;
542 e = em;
543 f = LikeMod3(mu,b,e,x,bm,em,u,v);
544 }
545
546
547 if(what == 3) {
548 if(mu == 0.0){
549 b = ((bm-u)+TMath::Sqrt((bm-u)*(bm-u)+4*x*u))/2.;
550 e = em;
551 f = LikeMod3(mu,b,e,x,bm,em,u,v);
552 } else {
553 Double_t temp[3];
554 temp[0] = mu*mu*v+u;
555 temp[1] = mu*mu*mu*v*v+mu*v*u-mu*mu*v*em+mu*v*bm-2*u*em;
556 temp[2] = mu*mu*v*v*bm-mu*v*u*em-mu*v*bm*em+u*em*em-mu*mu*v*v*x;
557 e = (-temp[1]+TMath::Sqrt(temp[1]*temp[1]-4*temp[0]*temp[2]))/2/temp[0];
558 b = bm-(u*(em-e))/v/mu;
559 f = LikeMod3(mu,b,e,x,bm,em,u,v);
560 }
561 }
562
563 return f;
564}
565
566//____________________________________________________________________
567Double_t MRolke::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)
568{
569 // Profile Likelihood function for MODEL 3:
570 // Gauss background/Gauss Efficiency
571
572 return 2*(x * TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)-1.837877-TMath::Log(u)/2-(bm-b)*(bm-b)/u/2-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
573}
574
575//____________________________________________________________________
576Double_t MRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what)
577{
578 // Calculates the Profile Likelihood for MODEL 4:
579 // Poiss background/Efficiency known
580 // what = 1: Maximum likelihood estimate is returned
581 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
582 // what = 3: Profile Likelihood of Test hypothesis is returned
583 // otherwise parameters as described in the beginning of the class)
584
585 Double_t f = 0.0;
586
587 if (what == 1) f = x-y/tau;
588 if (what == 2) {
589 mu = x-y/tau;
590 b = Double_t(y)/tau;
591 f = LikeMod4(mu,b,x,y,tau);
592 }
593 if (what == 3) {
594 if (mu == 0.0) {
595 b = Double_t(x+y)/(1+tau);
596 f = LikeMod4(mu,b,x,y,tau);
597 } else {
598 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);
599 f = LikeMod4(mu,b,x,y,tau);
600 }
601 }
602 return f;
603}
604
605//___________________________________________________________________
606Double_t MRolke::LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau)
607{
608 // Profile Likelihood function for MODEL 4:
609 // Poiss background/Efficiency known
610
611 return 2*(x*TMath::Log(mu+b)-(mu+b)-LogFactorial(x)+y*TMath::Log(tau*b)-tau*b-LogFactorial(y) );
612}
613
614//___________________________________________________________________
615Double_t MRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what)
616{
617 // Calculates the Profile Likelihood for MODEL 5:
618 // Gauss background/Efficiency known
619 // what = 1: Maximum likelihood estimate is returned
620 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
621 // what = 3: Profile Likelihood of Test hypothesis is returned
622 // otherwise parameters as described in the beginning of the class)
623
624 Double_t u=sdb*sdb;
625 Double_t f = 0;
626
627 if(what == 1) {
628 f = x - bm;
629 }
630 if(what == 2) {
631 mu = x-bm;
632 b = bm;
633 f = LikeMod5(mu,b,x,bm,u);
634 }
635
636 if (what == 3) {
637 b = ((bm-u-mu)+TMath::Sqrt((bm-u-mu)*(bm-u-mu)-4*(mu*u-mu*bm-u*x)))/2;
638 f = LikeMod5(mu,b,x,bm,u);
639 }
640 return f;
641}
642
643//_______________________________________________________________________
644Double_t MRolke::LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u)
645{
646 // Profile Likelihood function for MODEL 5:
647 // Gauss background/Efficiency known
648
649 return 2*(x*TMath::Log(mu+b)-(mu + b)-LogFactorial(x)-0.9189385-TMath::Log(u)/2-((bm-b)*(bm-b))/u/2);
650}
651
652//_______________________________________________________________________
653Double_t MRolke::EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t e, Double_t b, Int_t m, Int_t what)
654{
655 // Calculates the Profile Likelihood for MODEL 6:
656 // Gauss known/Efficiency binomial
657 // what = 1: Maximum likelihood estimate is returned
658 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
659 // what = 3: Profile Likelihood of Test hypothesis is returned
660 // otherwise parameters as described in the beginning of the class)
661
662 Double_t coef[4],roots[3];
663 Double_t f = 0.;
664 Double_t zm = Double_t(z)/m;
665
666 if(what==1){
667 f = (x-b)/zm;
668 }
669
670 if(what==2){
671 mu = (x-b)/zm;
672 e = zm;
673 f = LikeMod6(mu,b,e,x,z,m);
674 }
675 if(what == 3){
676 if(mu==0){
677 e = zm;
678 } else {
679 coef[3] = mu*mu;
680 coef[2] = mu * b - mu * x - mu*mu - mu * m;
681 coef[1] = mu * x - mu * b + mu * z - m * b;
682 coef[0] = b * z;
683 TMath::RootsCubic(coef,roots[0],roots[1],roots[2]);
684 e = roots[1];
685 }
686 f =LikeMod6(mu,b,e,x,z,m);
687 }
688 return f;
689}
690
691//_______________________________________________________________________
692Double_t MRolke::LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m)
693{
694 // Profile Likelihood function for MODEL 6:
695 // background known/ Efficiency binomial
696
697 Double_t f = 0.0;
698
699 if (z > 100 || m > 100) {
700 f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+(m*TMath::Log(m) - m)-(z*TMath::Log(z) - z) - ((m-z)*TMath::Log(m-z) - m + z)+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
701 } else {
702 f = 2*(x*TMath::Log(e*mu+b)-(e*mu+b)-LogFactorial(x)+TMath::Log(TMath::Factorial(m))-TMath::Log(TMath::Factorial(z))-TMath::Log(TMath::Factorial(m-z))+z*TMath::Log(e)+(m-z)*TMath::Log(1-e));
703 }
704 return f;
705}
706
707
708//___________________________________________________________________________
709Double_t MRolke::EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t e, Double_t sde, Double_t b, Int_t what)
710{
711 // Calculates the Profile Likelihood for MODEL 7:
712 // background known/Efficiency Gauss
713 // what = 1: Maximum likelihood estimate is returned
714 // what = 2: Profile Likelihood of Maxmimum Likelihood estimate is returned.
715 // what = 3: Profile Likelihood of Test hypothesis is returned
716 // otherwise parameters as described in the beginning of the class)
717
718 Double_t v=sde*sde;
719 Double_t f = 0.;
720
721 if(what == 1) {
722 f = (x-b)/em;
723 }
724
725 if(what == 2) {
726 mu = (x-b)/em;
727 e = em;
728 f = LikeMod7(mu, b, e, x, em, v);
729 }
730
731 if(what == 3) {
732 if(mu==0) {
733 e = em;
734 } else {
735 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;
736 }
737 f = LikeMod7(mu, b, e, x, em, v);
738 }
739
740 return f;
741}
742
743//___________________________________________________________________________
744Double_t MRolke::LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v)
745{
746 // Profile Likelihood function for MODEL 6:
747 // background known/ Efficiency binomial
748
749 return 2*(x*TMath::Log(e*mu+b)-(e*mu + b)-LogFactorial(x)-0.9189385-TMath::Log(v)/2-(em-e)*(em-e)/v/2);
750}
751
752//______________________________________________________________________
753Double_t MRolke::EvalPolynomial(Double_t x, const Int_t coef[], Int_t N)
754{
755 // evaluate polynomial
756
757 const Int_t *p;
758 p = coef;
759 Double_t ans = *p++;
760 Int_t i = N;
761
762 do
763 ans = ans * x + *p++;
764 while( --i );
765
766 return ans;
767}
768
769//______________________________________________________________________
770Double_t MRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
771{
772 // evaluate mononomial
773
774 Double_t ans;
775 const Int_t *p;
776
777 p = coef;
778 ans = x + *p++;
779 Int_t i = N-1;
780
781 do
782 ans = ans * x + *p++;
783 while( --i );
784
785 return ans;
786}
787//--------------------------------------------------------------------------
788//
789// Uses Stirling-Wells formula: ln(N!) ~ N*ln(N) - N + 0.5*ln(2piN)
790// if N exceeds 70, otherwise use the TMath functions
791//
792//
793Double_t MRolke::LogFactorial(Int_t n)
794{
795 if (n>69)
796 return n*TMath::Log(n)-n + 0.5*TMath::Log(TMath::TwoPi()*n);
797 else
798 return TMath::Log(TMath::Factorial(n));
799}
Note: See TracBrowser for help on using the repository browser.