From b52f09bd049aa122d0c15498de454555538b34f0 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Sun, 25 Nov 2018 11:38:19 -0600 Subject: update likelihood to make sure we integrate over at least 100 points --- src/likelihood.c | 11 +++++++++++ src/likelihood.h | 3 +++ 2 files changed, 14 insertions(+) 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 /* 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 -- cgit