Skip to content

Fix Crystal Ball integral near n = 1 for RooCBShape and RooCrystalBall #19602

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Fix Crystal Ball integral near n = 1 for RooCBShape and RooCrystalBall #19602

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
7f96d74
[RF] Fix RooCrystalBall integral near n = 1
MartinDuyTat Aug 11, 2025
a34b672
fixup! [RF] Fix RooCrystalBall integral near n = 1
MartinDuyTat Aug 13, 2025
382066e
[RF] Fix RooCBShape integral near n = 1
MartinDuyTat Aug 11, 2025
ba957cc
fixup! [RF] Fix RooCBShape integral near n = 1
MartinDuyTat Aug 13, 2025
b464b36
[RF] Fix ROOT::Math::crystalball_integral near n = 1
MartinDuyTat Aug 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions math/mathcore/src/ProbFuncMathCore.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");

Expand All @@ -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));
}
Expand Down
5 changes: 4 additions & 1 deletion roofit/roofit/src/RooCrystalBall.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
10 changes: 8 additions & 2 deletions roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
Expand All @@ -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);
}
Expand Down
Loading