diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 11 | ||||
| -rw-r--r-- | src/likelihood.h | 3 |
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 |