diff --git a/math/physics/inc/TRolke.h b/math/physics/inc/TRolke.h index 801d68e6f34bf..78ce61b24cfb4 100644 --- a/math/physics/inc/TRolke.h +++ b/math/physics/inc/TRolke.h @@ -65,11 +65,7 @@ class TRolke : public TObject // LIKELIHOOD ROUTINE 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); - //MODEL 1 - Double_t EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what); - 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); - 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); - Double_t LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); + //MODEL 1 is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm //MODEL 2 Double_t EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t sde, Double_t tau, Int_t what); @@ -131,8 +127,7 @@ class TRolke : public TObject } // The Set methods for the different models are described in Rolke.cxx - // model 1 - void SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m); + // model 1 is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm // model 2 void SetPoissonBkgGaussEff(Int_t x, Int_t y, Double_t em, Double_t tau, Double_t sde); diff --git a/math/physics/src/TRolke.cxx b/math/physics/src/TRolke.cxx index 39535830ed29e..b673e39d4b1e6 100644 --- a/math/physics/src/TRolke.cxx +++ b/math/physics/src/TRolke.cxx @@ -25,6 +25,7 @@ If unsure, first consider model 3, 4 or 5. 1: SetPoissonBkgBinomEff(x,y,z,tau,m) +\warning This model is removed due to an implementation error that may lead to wrong limits estimation and can't be fixed without massive changes in the algorithm ~~~ Background: Poisson Efficiency: Binomial @@ -190,30 +191,6 @@ TRolke::~TRolke() { } -//////////////////////////////////////////////////////////////////////////////// -/// Model 1: Background - Poisson, Efficiency - Binomial -/// - x : number of observed events in the experiment -/// - y : number of observed events in background region -/// - z : number of MC events observed -/// - tau : ratio parameter (read TRolke.cxx for details) -/// - m : number of MC events generated - -void TRolke::SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - SetModelParameters( - x , // Int_t x, - y , // Int_t y, - z , // Int_t z, - 0 , // Double_t bm, - 0 , // Double_t em, - 0 , // Double_t e, - 1 , // Int_t mid, - 0 , // Double_t sde, - 0 , // Double_t sdb, - tau, // Double_t tau, - 0 , // Double_t b, - m); // Int_t m -} //////////////////////////////////////////////////////////////////////////////// /// Model 2: Background - Poisson, Efficiency - Gaussian @@ -934,7 +911,10 @@ Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, { switch (mid) { case 1: - return EvalLikeMod1(mu, x, y, z, tau, m, what); + std::cerr << "TRolke is removed due to an implementation error " << + "that may lead to wrong limits estimation and can't be " << + "fixed without massive changes in the algorithm" << std::endl; + return 0; case 2: return EvalLikeMod2(mu, x, y, em, sde, tau, what); case 3: @@ -956,120 +936,10 @@ Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, return 0; } -//////////////////////////////////////////////////////////////////////////////// -/// Calculates the Profile Likelihood for MODEL 1: -/// Poisson background/ Binomial Efficiency -/// - what = 1: Maximum likelihood estimate is returned -/// - what = 2: Profile Likelihood of Maximum Likelihood estimate is returned. -/// - what = 3: Profile Likelihood of Test hypothesis is returned -/// otherwise parameters as described in the beginning of the class) - -Double_t TRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what) -{ - Double_t f = 0; - Double_t zm = Double_t(z) / m; - - if (what == 1) { - f = (x - y / tau) / zm; - } - - if (what == 2) { - mu = (x - y / tau) / zm; - Double_t b = y / tau; - Double_t e = zm; - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } - - if (what == 3) { - if (mu == 0) { - Double_t b = (x + y) / (1.0 + tau); - Double_t e = zm; - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } else { - Double_t e = 0; - Double_t b = 0; - ProfLikeMod1(mu, b, e, x, y, z, tau, m); - f = LikeMod1(mu, b, e, x, y, z, tau, m); - } - } - - return f; -} - -//////////////////////////////////////////////////////////////////////////////// -/// Profile Likelihood function for MODEL 1: -/// Poisson background/ Binomial Efficiency - -Double_t TRolke::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) -{ - double s = e*mu+b; - double lls = - s; - if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x); - double bg = tau*b; - double llb = -bg; - if ( y > 0) llb = y*TMath::Log( bg) - bg - LogFactorial(y); - - double lle = 0; // binomial log-like - if (z == 0) lle = m * TMath::Log(1-e); - else if ( z == m) lle = m * TMath::Log(e); - else lle = z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z); - - double f = 2*( lls + llb + lle); - return f; -} - - // this code is non-sense - // need to solve using Minuit struct LikeFunction1 { }; -//////////////////////////////////////////////////////////////////////////////// -/// Helper for calculation of estimates of efficiency and background for model 1 - -void TRolke::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) -{ - Double_t med = 0.0, fmid; - Int_t maxiter = 1000; - Double_t acc = 0.00001; - Double_t emin = ((m + mu * tau) - TMath::Sqrt((m + mu * tau) * (m + mu * tau) - 4 * mu * tau * z)) / 2 / mu / tau; - - Double_t low = TMath::Max(1e-10, emin + 1e-10); - Double_t high = 1 - 1e-10; - - for (Int_t i = 0; i < maxiter; i++) { - med = (low + high) / 2.; - - fmid = LikeGradMod1(med, mu, x, y, z, tau, m); - - if (high < 0.5) acc = 0.00001 * high; - else acc = 0.00001 * (1 - high); - - if ((high - low) < acc*high) break; - - if (fmid > 0) low = med; - else high = med; - } - - e = med; - Double_t eta = Double_t(z) / e - Double_t(m - z) / (1 - e); - - b = Double_t(y) / (tau - eta / mu); -} - -//////////////////////////////////////////////////////////////////////////////// -/// Gradient model likelihood - -Double_t TRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m) -{ - Double_t eta, etaprime, bprime, f; - eta = static_cast(z) / e - static_cast(m - z) / (1.0 - e); - etaprime = (-1) * (static_cast(m - z) / ((1.0 - e) * (1.0 - e)) + static_cast(z) / (e * e)); - Double_t b = y / (tau - eta / mu); - bprime = (b * b * etaprime) / mu / y; - f = (mu + bprime) * (x / (e * mu + b) - 1) + (y / b - tau) * bprime + eta; - return f; -} - //////////////////////////////////////////////////////////////////////////////// /// Calculates the Profile Likelihood for MODEL 2: /// Poisson background/ Gauss Efficiency diff --git a/tutorials/legacy/math/Rolke.C b/tutorials/legacy/math/Rolke.C index eb78e0a322d06..2a42ac11bcf79 100644 --- a/tutorials/legacy/math/Rolke.C +++ b/tutorials/legacy/math/Rolke.C @@ -41,31 +41,6 @@ void Rolke() Double_t ul ; // upper limit Double_t ll ; // lower limit -//----------------------------------------------- -// Model 1 assumes: -// -// Poisson uncertainty in the background estimate -// Binomial uncertainty in the efficiency estimate -// - cout << endl<<" ======================================================== " <