diff --git a/math/mathcore/src/ProbFuncMathCore.cxx b/math/mathcore/src/ProbFuncMathCore.cxx index 7ac7d63f0e5e4..201fe10728b0d 100644 --- a/math/mathcore/src/ProbFuncMathCore.cxx +++ b/math/mathcore/src/ProbFuncMathCore.cxx @@ -107,7 +107,7 @@ namespace Math { MATH_ERROR_MSG("crystalball_integral", "CrystalBall function not defined at alpha=0"); return 0.; } - bool useLog = (n == 1.0); + bool useLog = (std::abs(n - 1.) < 1.0e-5); if (n <= 0) MATH_WARN_MSG("crystalball_integral", "No physical meaning when n<=0"); @@ -132,11 +132,13 @@ namespace Math { double B = r - abs_alpha; if (!useLog) { - intpow = A * (1 - std::pow(r / (B - z), n - 1)) / (n - 1); + intpow = A * (1. - std::pow(r / (B - z), n - 1.)) / (n - 1.); } else { // for n=1 the primitive of 1/x is log(x) - intpow = A * std::pow(r, n - 1) * (std::log(B - z) - std::log(r)); + double log_B_z = std::log(B - z); + double log_r = std::log(r); + intpow = A * std::pow(r, n - 1.) * (log_B_z - log_r + 0.5 * (1. - n) * (log_B_z * log_B_z - log_r * log_r)); } intgaus = sqrtpiover2 * (1. + ROOT::Math::erf(abs_alpha * oneoversqrt2)); } diff --git a/roofit/roofit/src/RooCrystalBall.cxx b/roofit/roofit/src/RooCrystalBall.cxx index bb0dc10722d8c..1a8803003aa68 100644 --- a/roofit/roofit/src/RooCrystalBall.cxx +++ b/roofit/roofit/src/RooCrystalBall.cxx @@ -185,7 +185,10 @@ inline double integrateTailLogVersion(double sigma, double alpha, double n, doub double a = std::pow(r, n) * exp(-0.5 * alpha * alpha); double b = r - alpha; - return a * sigma * (log(b - tmin) - log(b - tmax)); + double log_b_tmin = log(b - tmin); + double log_b_tmax = log(b - tmax); + + return a * sigma * (log_b_tmin - log_b_tmax + 0.5 * (1.0 - n) * (log_b_tmin * log_b_tmin - log_b_tmax * log_b_tmax)); } inline double integrateTailRegular(double sigma, double alpha, double n, double tmin, double tmax) diff --git a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h index a573e5c38ff53..a06933d074c3c 100644 --- a/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h +++ b/roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h @@ -712,7 +712,10 @@ inline double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double b = r - absAlpha; if (useLog) { - result += a * std::pow(r, n - 1) * sig * (std::log(b - tmin) - std::log(b - tmax)); + double log_b_tmin = std::log(b - tmin); + double log_b_tmax = std::log(b - tmax); + result += a * std::pow(r, n - 1) * sig * + (log_b_tmin - log_b_tmax + 0.5 * (1.0 - n) * (log_b_tmin * log_b_tmin - log_b_tmax * log_b_tmax)); } else { result += a * sig / (1.0 - n) * (std::pow(r / (b - tmin), n - 1.0) - std::pow(r / (b - tmax), n - 1.0)); } @@ -723,7 +726,10 @@ inline double cbShapeIntegral(double mMin, double mMax, double m0, double sigma, double term1 = 0.0; if (useLog) { - term1 = a * std::pow(r, n - 1) * sig * (std::log(b - tmin) - std::log(r)); + double log_b_tmin = std::log(b - tmin); + double log_r = std::log(r); + term1 = a * std::pow(r, n - 1) * sig * + (log_b_tmin - log_r + 0.5 * (1.0 - n) * (log_b_tmin * log_b_tmin - log_r * log_r)); } else { term1 = a * sig / (1.0 - n) * (std::pow(r / (b - tmin), n - 1.0) - 1.0); }