From fdc9ac8e7e257e080bd30e3d1dcd7e888de0f329 Mon Sep 17 00:00:00 2001 From: "Ryan M. L. McFadden" Date: Mon, 10 Jan 2022 14:20:08 +0100 Subject: [PATCH 1/2] improvements to PTheory::SkewedGauss - Divide the function evaluation into even/odd frequency components. Some additional "helper" terms have been added to aid in this. - Use a better value to check for floating-point overflow in the evaluation of ROOT::conf_hyperg and return std::numeric_limits::max() when necessary. This extends the function's range of validity to arbitrary time (whereas the previous implementation would fail loudly for large sigma+/-). - Format the src with clang-format. --- src/classes/PTheory.cpp | 82 +++++++++++++++++++++++++---------------- 1 file changed, 51 insertions(+), 31 deletions(-) diff --git a/src/classes/PTheory.cpp b/src/classes/PTheory.cpp index 8279e55f..42ba4228 100644 --- a/src/classes/PTheory.cpp +++ b/src/classes/PTheory.cpp @@ -2160,51 +2160,71 @@ Double_t PTheory::InternalBessel(Double_t t, const PDoubleVector& paramValues, c * \param paramValues parameter values * \param funcValues vector with the functions (i.e. functions of the parameters) */ -Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const -{ - // expected parameters: phase frequency sigma- sigma+ [tshift] - +Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector ¶mValues, + const PDoubleVector &funcValues) const { + // Expected parameters: phase, frequency, sigma-, sigma+, [tshift]. + // To be stored in the array "val" as: + // val[0] = phase + // val[1] = frequency + // val[2] = sigma- + // val[3] = sigma+ + // val[4] = tshift [optional] Double_t val[5]; - Double_t skg; + // Check that we have the correct number of fit parameters. assert(fParamNo.size() <= 5); - // check if FUNCTIONS are used - for (UInt_t i=0; i= 25.0) || (zm >= 25.0)) // needed to prevent crash of 1F1 - skg = 2.0e6; - else if (fabs(val[2]) == fabs(val[3])) // sigma+ == sigma- -> Gaussian - skg = TMath::Cos(phase+freq*tt) * gp; - else - skg = TMath::Cos(phase+freq*tt) * (wm*gm + wp*gp) + - TMath::Sin(phase+freq*tt) * (wm*gm*2.0*zm/SQRT_PI*ROOT::Math::conf_hyperg(0.5,1.5,zm*zm) - - wp*gp*2.0*zp/SQRT_PI*ROOT::Math::conf_hyperg(0.5,1.5,zp*zp)); + // Evalute the EVEN frequency component of the skewed Gaussian. + Double_t skg_cos = TMath::Cos(phase + freq * tt) * (w_m * g_m + w_p * g_p); - return skg; + // Evalute the ODD frequency component of the skewed Gaussian. + constexpr Double_t z_max = 26.7776; + // 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". + 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))); + + // return the skewed Gaussian: skg = skg_cos + skg_sin. + return skg_cos + skg_sin; } //-------------------------------------------------------------------------- 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 2/2] 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); } //--------------------------------------------------------------------------