180 fNumWarningsDeprecated1(0),
181 fNumWarningsDeprecated2(0)
376 std::cerr <<
"TRolke - Error: Model id "<<
f_mid<<std::endl;
378 std::cerr <<
"TRolke - Please specify a model with e.g. 'SetGaussBkgGaussEff' (read the docs in Rolke.cxx )"<<std::endl;
383 ComputeInterval(
f_x,
f_y,
f_z,
f_bm,
f_em,
f_e,
f_mid,
f_sde,
f_sdb,
f_tau,
f_b,
f_m);
389 std::cerr <<
"TRolke - Warning: no limits found" <<std::endl;
435 std::cerr <<
"TRolke::GetBackground(): Model NR: " <<
436 f_mid <<
" unknown"<<std::endl;
456 ComputeInterval(loop_x,
f_y,
f_z,
f_bm,
f_em,
f_e,
f_mid,
f_sde,
f_sdb,
f_tau,
f_b,
f_m);
461 if (loop_x > (background + 1)) {
462 if ((weightSum > (1 - pPrecision)) || (weight < 1
e-12))
break;
491 if (weightSum >= integral) {
499 ComputeInterval(loop_x,
f_y,
f_z,
f_bm,
f_em,
f_e,
f_mid,
f_sde,
f_sdb,
f_tau,
f_b,
f_m);
519 while (loop_x <= loop_max) {
526 if (loop_x >= loop_max) {
527 std::cout <<
"internal error finding maximum of distribution" << std::endl;
533 ComputeInterval(loop_x,
f_y,
f_z,
f_bm,
f_em,
f_e,
f_mid,
f_sde,
f_sdb,
f_tau,
f_b,
f_m);
551 int rolke_ncrit = -1;
554 maxj = 1000 + (
Int_t)background;
556 for (j = 0;j < maxj;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);
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;
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;
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;
637Double_t TRolke::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){
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;
658 return ComputeInterval(
f_x,
f_y,
f_z,
f_bm,
f_em,
f_e,
f_mid,
f_sde,
f_sdb,
f_tau,
f_b,
f_m);
675void TRolke::SetModelParameters(
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)
694 SetModelParameters(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
713Double_t TRolke::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)
719 limit[1] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
731 limit[1] =
Interval(trial_x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
732 if (limit[1] > 0) done = 1;
754Double_t TRolke::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)
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;
760 Double_t maxiter = 1000, acc = 0.00001;
764 if ((mid != 3) && (mid != 5)) bm =
y;
765 if ((mid == 3) || (mid == 5)) {
766 if (bm == 0) bm = 0.00001;
769 if ((mid == 6) || (mid == 7)) {
770 if (bm == 0) bm = 0.00001;
773 if ((mid <= 2) || (mid == 4)) bp = 1;
776 if (bp == 1 &&
x == 0 && bm > 0) {
777 for (i = 0; i < 2; i++) {
779 tempxy[i] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
782 slope = tempxy[1] - tempxy[0];
783 limits[1] = tempxy[0] - slope;
785 if (limits[1] < 0) limits[1] = 0.0;
789 if (bp != 1 &&
x == 0) {
791 for (i = 0; i < 2; i++) {
793 tempxy[i] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
795 slope = tempxy[1] - tempxy[0];
796 limits[1] = tempxy[0] - slope;
798 if (limits[1] < 0) limits[1] = 0.0;
802 if (bp != 1 && bm == 0) {
803 for (i = 0; i < 2; i++) {
805 limits[1] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
806 tempxy[i] = limits[1];
808 slope = tempxy[1] - tempxy[0];
809 limits[1] = tempxy[0] - slope;
810 if (limits[1] < 0) limits[1] = 0;
814 if (
x == 0 && bm == 0) {
817 limits[1] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
818 tempxy[0] = limits[1];
821 limits[1] =
Interval(
x,
y, z, bm, em,
e, mid, sde, sdb, tau,
b,
m);
822 tempxy[1] = limits[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;
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);
834 f0 =
Likelihood(
test,
x,
y, z, bm, em, mid, sde, sdb, tau,
b,
m, 3);
836 if (
mu0 < 0) maximum = f0;
839 target = maximum - dchi2;
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;
861 fmid =
Likelihood(med,
x,
y, z, bm, em, mid, sde, sdb, tau,
b,
m, 3);
869 if ((high - low) < acc*high)
break;
883 ftest =
Likelihood(
test,
x,
y, z, bm, em, mid, sde, sdb, tau,
b,
m, 3);
884 if (ftest < target) {
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);
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);
908 if (high - low < acc*high)
break;
916 if ((mid == 4) || (mid == 5)) {
933Double_t TRolke::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)
951 std::cerr <<
"TRolke::Likelihood(...): Model NR: " <<
952 f_mid <<
" unknown"<<std::endl;
973 f = (
x -
y / tau) / zm;
977 mu = (
x -
y / tau) / zm;
1017 double f = 2*( lls + llb + lle);
1023struct LikeFunction1 {
1032 Int_t maxiter = 1000;
1034 Double_t emin = ((
m + mu * tau) -
TMath::Sqrt((
m + mu * tau) * (
m + mu * tau) - 4 * mu * tau * z)) / 2 / mu / tau;
1039 for (
Int_t i = 0; i < maxiter; i++) {
1040 med = (low + high) / 2.;
1044 if (high < 0.5) acc = 0.00001 * high;
1045 else acc = 0.00001 * (1 - high);
1047 if ((high - low) < acc*high)
break;
1049 if (fmid > 0) low = med;
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));
1068 bprime = (
b *
b * etaprime) / mu /
y;
1069 f = (mu + bprime) * (
x / (
e * mu +
b) - 1) + (
y /
b - tau) * bprime + eta;
1088 f = (
x -
y / tau) / em;
1092 mu = (
x -
y / tau) / em;
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;
1113 if (
v > 0)
b =
y / (tau + (em -
e) / mu /
v);
1135 if (
v > 0) lle = - 0.9189385 -
TMath::Log(
v) / 2 - (em -
e)*(em -
e) /
v / 2;
1137 return 2*( lls + llb + lle);
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;
1200 if ( u > 0) llb = - 0.9189385 -
TMath::Log(u) / 2 - (bm -
b)*(bm -
b) / u / 2;
1202 if (
v > 0) lle = - 0.9189385 -
TMath::Log(
v) / 2 - (em -
e)*(em -
e) /
v / 2;
1204 return 2*( lls + llb + lle);
1220 if (
what == 1)
f =
x -
y / tau;
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);
1251 return 2*( lls + llb);
1277 Double_t b = ((bm - u - mu) +
TMath::Sqrt((bm - u - mu) * (bm - u - mu) - 4 * (mu * u - mu * bm - u *
x))) / 2;
1293 if ( u > 0) llb = - 0.9189385 -
TMath::Log(u) / 2 - (bm -
b)*(bm -
b) / u / 2;
1295 return 2*( lls + llb);
1327 coef[2] = mu *
b - mu *
x - mu * mu - mu *
m;
1328 coef[1] = mu *
x - mu *
b + mu * z -
m *
b;
1353 return 2* (lls + lle);
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;
1404 if (
v > 0) lle = - 0.9189385 -
TMath::Log(
v) / 2 - (em -
e)*(em -
e) /
v / 2;
1406 return 2*( lls + lle);
1420 ans = ans *
x + *p++;
1439 ans = ans *
x + *p++;
This class computes confidence intervals for the rate of a Poisson process in the presence of uncerta...
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.
void SetGaussBkgKnownEff(Int_t x, Double_t bm, Double_t sdb, Double_t e)
Model 5: Background - Gaussian, Efficiency - known (x,bm,sdb,e.
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 'Interval'.
bool GetLimitsML(Double_t &low, Double_t &high, Int_t &out_x)
get the upper and lower limits for the most likely outcome.
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.
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.
bool GetCriticalNumber(Int_t &ncrit, Int_t maxtry=-1)
get the value of x corresponding to rejection of the null hypothesis.
static Double_t EvalPolynomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate polynomial.
void SetKnownBkgBinomEff(Int_t x, Int_t z, Int_t m, Double_t b)
Model 6: Background - known, Efficiency - Binomial (x,z,m,b)
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...
Int_t fNumWarningsDeprecated1
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.
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.
void SetPoissonBkgKnownEff(Int_t x, Int_t y, Double_t tau, Double_t e)
Model 4: Background - Poisson, Efficiency - known (x,y,tau,e)
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.
Double_t GetUpperLimit()
Calculate and get upper limit for the pre-specified model.
void SetBounding(const bool bnd)
void SetKnownBkgGaussEff(Int_t x, Double_t em, Double_t sde, Double_t b)
Model 7: Background - known, Efficiency - Gaussian (x,em,sde,b)
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)
void SetPoissonBkgGaussEff(Int_t x, Int_t y, Double_t em, Double_t tau, Double_t sde)
Model 2: Background - Poisson, Efficiency - Gaussian.
static Double_t EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
Evaluate mononomial.
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.
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.
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.
TRolke(Double_t CL=0.9, Option_t *option="")
Constructor with optional Confidence Level argument.
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.
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.
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.
void SetSwitch(bool bnd)
Deprecated name for SetBounding.
void SetModelParameters()
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.
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)
Deprecated and error prone model selection interface.
bool GetLimits(Double_t &low, Double_t &high)
Calculate and get the upper and lower limits for the pre-specified model.
void SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
Model 1: Background - Poisson, Efficiency - Binomial.
Double_t LogFactorial(Int_t n)
LogFactorial function (use the logGamma function via the relation Gamma(n+1) = n!
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.
virtual ~TRolke()
Destructor.
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.
Double_t GetLowerLimit()
Calculate and get lower limit for the pre-specified model.
Int_t fNumWarningsDeprecated2
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.
Double_t GetBackground()
Return a simple background value (estimate/truth) given the pre-specified model.
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.
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.
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.
void Print(Option_t *) const
Dump internals. Print members.
static constexpr double s
static constexpr double mu0
Short_t Max(Short_t a, Short_t b)
Double_t PoissonI(Double_t x, Double_t par)
Compute the Discrete Poisson distribution function for (x,par).
Double_t Sqrt(Double_t x)
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
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.
Double_t ChisquareQuantile(Double_t p, Double_t ndf)
Evaluate the quantiles of the chi-squared probability distribution function.