update likelihood to make sure we integrate over at least 100 points

2 files changed, 14 insertions, 0 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 724cc94..ab39649 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -761,6 +761,8 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     /* Number of points to sample along the particle track. */
     npoints = (size_t) (range/dx + 0.5);
 
+    if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS;
+
     /* Calculate total energy */
     E0 = T0 + p->mass;
     p0 = sqrt(E0*E0 - p->mass*p->mass);
@@ -784,12 +786,21 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
              * the electromagnetic shower. */
             npoints_shower = (size_t) ((xhi-xlo)/dx + 0.5);
 
+            if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS;
+
             x = malloc(sizeof(double)*npoints_shower);
             pdf = malloc(sizeof(double)*npoints_shower);
             for (i = 0; i < npoints_shower; i++) {
                 x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1);
                 pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
             }
+
+            /* Make sure the PDF is normalized. */
+            double norm = trapz(pdf,x[1]-x[0],npoints_shower);
+            for (i = 0; i < npoints_shower; i++) {
+                pdf[i] /= norm;
+            }
+
         }
     }
 
diff --git a/src/likelihood.h b/src/likelihood.h
index 129c313..f087e90 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -4,6 +4,9 @@
 #include "event.h"
 #include <stddef.h> /* for size_t */
 
+/* Minimum number of points in the track integral. */
+#define MIN_NPOINTS 100
+
 /* To avoid having to allocate new arrays every time we evaluate the likelihood
  * function, we just allocate static arrays with `MAX_POINTS` elements. */
 #define MAX_NPOINTS 1000