Changeset 8105
- Timestamp:
- 10/17/06 17:32:52 (18 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtools/MRolke.cc
r8087 r8105 1 // @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1. 1 2006-10-17 12:51:31meyer Exp $1 // @(#)root/physics:$Name: not supported by cvs2svn $:$Id: MRolke.cc,v 1.2 2006-10-17 16:32:52 meyer Exp $ 2 2 // Author: Jan Conrad 9/2/2004 3 3 … … 115 115 MRolke::MRolke(Double_t CL, Option_t * /*option*/) 116 116 { 117 118 119 120 121 122 // 1: bounded likelihood117 //constructor 118 fUpperLimit = 0.0; 119 fLowerLimit = 0.0; 120 fCL = CL; 121 fSwitch = 0; // 0: unbounded likelihood 122 // 1: bounded likelihood 123 123 } 124 124 … … 132 132 Double_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 133 { 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 }153 154 155 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 156 } 157 157 … … 162 162 Double_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 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) { 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) { 324 283 high = med; 325 284 fhigh = fmid; 326 } else {285 } else { 327 286 low = med; 328 287 flow = fmid; 329 } 330 if (high-low < acc*high) break; 331 } 332 limits[1] = med; 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 333 334 334 done: 335 335 336 337 338 339 340 341 342 343 344 345 346 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 347 } 348 348 … … 351 351 Double_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 352 { 353 354 355 356 357 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 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 368 } 369 369 … … 371 371 Double_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 372 { 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 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 407 } 408 408 … … 410 410 Double_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 411 { 412 413 414 415 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 416 } 417 417 … … 419 419 void 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 420 { 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 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 449 } 450 450 … … 453 453 { 454 454 //gradient model 455 456 457 458 459 460 461 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 462 } 463 463 … … 465 465 Double_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 466 { 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 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 493 } else { 494 495 496 497 498 499 500 501 502 503 504 } 505 506 507 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 508 } 509 509 … … 512 512 { 513 513 // Profile Likelihood function for MODEL 2: 514 515 516 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 517 } 518 518 … … 521 521 Double_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 522 { 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 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 551 f = LikeMod3(mu,b,e,x,bm,em,u,v); 552 553 554 555 556 557 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 558 b = bm-(u*(em-e))/v/mu; 559 559 f = LikeMod3(mu,b,e,x,bm,em,u,v); 560 561 562 563 560 } 561 } 562 563 return f; 564 564 } 565 565 … … 567 567 Double_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 568 { 569 570 571 572 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 573 } 574 574 … … 576 576 Double_t MRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Double_t b, Int_t what) 577 577 { 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 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 603 } 604 604 … … 606 606 Double_t MRolke::LikeMod4(Double_t mu,Double_t b,Int_t x,Int_t y,Double_t tau) 607 607 { 608 609 610 611 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 612 } 613 613 … … 615 615 Double_t MRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Double_t b, Int_t what) 616 616 { 617 618 619 620 621 622 623 624 625 626 627 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 628 f = x - bm; 629 630 631 632 633 634 635 636 637 638 639 640 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 641 } 642 642 … … 644 644 Double_t MRolke::LikeMod5(Double_t mu,Double_t b,Int_t x,Double_t bm,Double_t u) 645 645 { 646 647 648 649 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 650 } 651 651 … … 653 653 Double_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 654 { 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 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 689 } 690 690 … … 692 692 Double_t MRolke::LikeMod6(Double_t mu,Double_t b,Double_t e,Int_t x,Int_t z,Int_t m) 693 693 { 694 695 696 697 698 699 700 701 702 703 704 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 705 } 706 706 … … 709 709 Double_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 710 { 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 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 741 } 742 742 … … 744 744 Double_t MRolke::LikeMod7(Double_t mu,Double_t b,Double_t e,Int_t x,Double_t em,Double_t v) 745 745 { 746 747 748 749 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 750 } 751 751 … … 753 753 Double_t MRolke::EvalPolynomial(Double_t x, const Int_t coef[], Int_t N) 754 754 { 755 756 757 758 759 760 761 762 763 764 765 766 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 767 } 768 768 … … 770 770 Double_t MRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N) 771 771 { 772 773 774 775 776 777 778 779 780 781 782 783 784 785 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 786 } 787 787 //-------------------------------------------------------------------------- … … 793 793 Double_t MRolke::LogFactorial(Int_t n) 794 794 { 795 796 797 798 799 } 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 TracChangeset
for help on using the changeset viewer.