From d319be4e9710a4c87a6217ba0a7b06a73b9e65cc Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 26 Aug 2019 10:25:41 -0500 Subject: fix a bug in get_{shower,delta_ray}_weights() causing crash This commit fixes a bug in get_shower_weights() and get_delta_ray_weights() which was causing an inf value to propagate and cause the fitter to crash. The problem came because due to floating point roundoff the cdf value at the end of the loop was slightly greater than the last cdf value we wanted which was causing it to get mapped to cos(theta) = -1 (I think?) and then subsequently get interpolated to an infinite value for xcdf. The fix is just to make sure that the x coordinate is always between x1 and x2. --- src/likelihood.c | 14 ++++++++++++++ 1 file changed, 14 insertions(+) (limited to 'src') diff --git a/src/likelihood.c b/src/likelihood.c index 4a0e755..5ca9c8d 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -358,6 +358,13 @@ int get_delta_ray_weights(particle *p, double range, double distance_to_psup, do cdf = cdf1 + i*delta; cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + + /* Occasionally due to floating point rounding error, it's possible + * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, + * we just make sure that xcdf is bounded by x1 and x2. */ + if (xcdf < x1) xcdf = x1; + if (xcdf > x2) xcdf = x2; + if (i > 0) { /* For each interval we approximate the integral as: * @@ -442,6 +449,13 @@ int get_shower_weights(particle *p, double distance_to_psup, double *x, double * cdf = cdf1 + i*delta; cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + + /* Occasionally due to floating point rounding error, it's possible + * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, + * we just make sure that xcdf is bounded by x1 and x2. */ + if (xcdf < x1) xcdf = x1; + if (xcdf > x2) xcdf = x2; + if (i > 0) { /* For each interval we approximate the integral as: * -- cgit