LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

fix the discontinuity encountered at large t when sigma+/- is large

Browse Source 
- Though the discontinuity encountered in the previous version is small (because of the large Gaussain damping terms), it introduced considerable numeric instability when fitting. Consequently, simply zeroing the already heavily damped is sensible and yeilds smooth behaviour of the function.
- Also ensure that the return value for the odd frequency component is always finite.
This commit is contained in:
Ryan M. L. McFadden 2022-01-10 17:13:11 +01:00 committed by Andreas Suter
parent aacad8e7ad
commit 0c4b33b056
1 changed files with 11 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/classes/PTheory.cpp
View File

@ -2213,18 +2213,20 @@ Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector &paramValues,
// Note: the check against z_max is needed to prevent floating-point overflow
// in the return value of ROOT::Math::conf_hyperg(1/2, 3/2, z * z)
// (i.e., confluent hypergeometric function of the first kind, 1F1).
// In the case that z > z_max, use the maximum value of "Double_t".
// In the case that z > z_max, return zero (otherwise there is some
// numeric discontinuity at later times).
Double_t skg_sin =
TMath::Sin(phase + freq * tt) *
((w_m * g_m * 2.0 * z_m / SQRT_PI) *
(z_m > z_max ? std::numeric_limits<Double_t>::max()
: ROOT::Math::conf_hyperg(0.5, 1.5, z_m * z_m)) -
(w_p * g_p * 2.0 * z_p / SQRT_PI) *
(z_p > z_max ? std::numeric_limits<Double_t>::max()
: ROOT::Math::conf_hyperg(0.5, 1.5, z_p * z_p)));
((z_m > z_max) or (z_p > z_max)
? 0.0
: (w_m * g_m * 2.0 * z_m / SQRT_PI) *
ROOT::Math::conf_hyperg(0.5, 1.5, z_m * z_m) -
(w_p * g_p * 2.0 * z_p / SQRT_PI) *
ROOT::Math::conf_hyperg(0.5, 1.5, z_p * z_p));
// return the skewed Gaussian: skg = skg_cos + skg_sin.
return skg_cos + skg_sin;
// Return the skewed Gaussian: skg = skg_cos + skg_sin.
// Also check that skg_sin is finite!
return skg_cos + (std::isfinite(skg_sin) ? skg_sin : 0.0);
}
//--------------------------------------------------------------------------