From 27fe3ff5fd6fea66d06ade5b21bea879745c9226 Mon Sep 17 00:00:00 2001 From: "Ryan M. L. McFadden" Date: Mon, 10 Jan 2022 17:13:11 +0100 Subject: [PATCH] fix the discontinuity encountered at large t when sigma+/- is large - 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. --- src/classes/PTheory.cpp | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 42ba4228..ca64bdcc 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -2213,18 +2213,20 @@ Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector ¶mValues, // 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::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::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); } //--------------------------------------------------------------------------