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

[MartinDuyTat]

This Pull request:

Fix bug in calculation of Crystal Ball integrals used in RooFit

Changes

• Analytical integral of RooCBShape
• Analytical integral of RooCrystalBall

Why?

In the integral of the Crystal Ball function, the analytical expression results in 0/0 when n = 1. To circumvent this, when abs(n - 1) < 1e-5 an alternative calculation is used which can be derived by a Taylor expansion.

In particular, the expression

r / (1.0 - n) * (std::pow(r / (b - tmin), n - 1.0) - std::pow(r / (b - tmax), n - 1.0));

is approximated to

std::pow(r, n) * (log(b - tmin) - log(b - tmax));

This is done by expanding the quantity inside the parentheses to first order in n - 1.

However, since there is a factor 1.0 - n in front of the parentheses, this expression effectively becomes a zeroth order expansion in n - 1. This is problematic because the pre-factor std::pow(r, n) is in fact a function of n, with a non-zero first order term in n - 1.

Therefore, the quantity inside the parentheses really should be expanded to second order in n - 1, otherwise the whole expression will have an incorrect value at first order in n - 1.

Here are plots that show the issue, and also the C++ scripts that produced the plots. I have scanned the relevant functions across n = 1 and one can clearly see that the current implementation is wrong, and this PR ensures that the gradient is correct near n = 1.

RooCrystalBall_fix.pdf
RooCBShape_fix.pdf
RooCBShape_fix.txt
RooCrystalBall_fix.txt

Checklist:

• tested changes locally

[ferdymercury]

Thanks!

• Where is the check abs(n-1)<1e-5 ?
• Can you also update math/mathcore/src/ProbFuncMathCore.cxx

[MartinDuyTat]

* Where is the check `abs(n-1)<1e-5` ?

The check is done on lines 270 and 274 in roofit/roofit/src/RooCrystalBall.cxx and line 691 in roofit/roofitcore/inc/RooFit/Detail/MathFuncs.h.

* Can you also update math/mathcore/src/ProbFuncMathCore.cxx

I don't think this one is necessary? On line 113 I only see bool useLog = (n == 1.0), and at this single point n = 1 I believe the logarithmic approximation is only needed at zeroth order.

[ferdymercury]

On line 113 I only see bool useLog = (n == 1.0), and at this single point n = 1 I believe the logarithmic approximation is only needed at zeroth order.

Thanks a lot for clarifying!
If you generate with your script the same canvas you attached, but using the function from ProbFundMathCore, does the function behave correctly in the range n = 1 - 1e-6 or does it show some artefacts ?
If it shows artefacts, we should try to make everything consistent by also changing bool useLog = abs(n-1) < 1e-5 and adding 2nd order, as in the other two files.

[MartinDuyTat]

On line 113 I only see bool useLog = (n == 1.0), and at this single point n = 1 I believe the logarithmic approximation is only needed at zeroth order.

Thanks a lot for clarifying! If you generate with your script the same canvas you attached, but using the function from ProbFundMathCore, does the function behave correctly in the range n = 1 - 1e-6 or does it show some artefacts ? If it shows artefacts, we should try to make everything consistent by also changing bool useLog = abs(n-1) < 1e-5 and adding 2nd order, as in the other two files.

Hi, for the function ROOT::Math::crystalball_integral I just scanned n - 1 across the range [-1e-6, 1e-6], in 1e-8 increments, using these parameters:

x = -2.0
alpha = 0.5
sigma = 1.0
mean = 0.0

In the first attachment you can see that there's no weird behaviour, other than what looks like float-point errors. I assume these small float-point errors is the reason why the RooFit implementation uses the logarithmic approximation in some finite range around n = 1. It looks like very small noise, as you can see in the second attachment where I've scanned around a larger range [-1e-5, 1e-5]. But I haven't tested the function for different shape parameters and I don't know if there are any cases where these float-point errors become large.

MathCore_crystalball_fix.pdf
MathCore_crystalball_fix_wide_range.pdf

[ferdymercury]

Thanks a lot!
Leaving @guitargeek comment on whether the Math:: function should also be made modified (switchToLogThreshold --> 1e-5 range) so that useLog is consistent across ROOT.

[guitargeek]

Thank you @MartinDuyTat for this great PR!

For consistency, I would indeed suggest to use the same tolerance for the TMath function before merging this PR:

• https://github.com/root-project/root/blob/master/math/mathcore/src/ProbFuncMathCore.cxx#L110

[github-actions]

Test Results

    15 files      15 suites   2d 19h 34m 15s ⏱️
 3 363 tests  3 344 ✅ 0 💤 19 ❌
49 110 runs  49 020 ✅ 0 💤 90 ❌

For more details on these failures, see this check.

Results for commit b464b36.

♻️ This comment has been updated with latest results.

[MartinDuyTat]

Thank you @MartinDuyTat for this great PR!

For consistency, I would indeed suggest to use the same tolerance for the TMath function before merging this PR:

* https://github.com/root-project/root/blob/master/math/mathcore/src/ProbFuncMathCore.cxx#L110

Hi @guitargeek I've just updated the TMath function to use the same tolerance and I've put in the same fix there as well. I'm also attaching a comparison between the original and new calculation.

MathCore_crystalball_fix.pdf