From c274fc18f49b70880f33cf762d0d1f78f037be26 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Thu, 23 May 2019 10:28:53 -0400 Subject: add 1e-10 to the angular pdfs to ensure that the CDF is strictly increasing --- src/likelihood.c | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) (limited to 'src/likelihood.c') diff --git a/src/likelihood.c b/src/likelihood.c index 61034a3..9a98562 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -211,8 +211,11 @@ particle *particle_init(int id, double T0, size_t n) * seems to be an OK assumption for high energy showers, but is not * exactly correct for showers from lower energy electrons or for * delta rays from lower energy muons. It seems good enough, but in - * the future it would be nice to parameterize this. */ - pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm; + * the future it would be nice to parameterize this. + * + * Note: We add EPSILON to the pdf since the cdf values are + * required to be strictly increasing in order to use gsl splines. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm + EPSILON; if (i > 0) p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; else @@ -242,8 +245,11 @@ particle *particle_init(int id, double T0, size_t n) * seems to be an OK assumption for high energy showers, but is not * exactly correct for showers from lower energy electrons or for * delta rays from lower energy muons. It seems good enough, but in - * the future it would be nice to parameterize this. */ - pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm; + * the future it would be nice to parameterize this. + * + * Note: We add EPSILON to the pdf since the cdf values are + * required to be strictly increasing in order to use gsl splines. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm + EPSILON; if (i > 0) p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; else -- cgit