speed up likelihood function and switch to using fixed dx

This commit speeds up the likelihood function by about ~20% by using the
precomputed track positions, directions, times, etc. instead of interpolating
them on the fly.

It also switches to computing the number of points to integrate along the track
by dividing the track length by a specified distance, currently set to 1 cm.
This should hopefully speed things up for lower energies and result in more
stable fits at high energies.

1 files changed, 25 insertions, 26 deletions

diff --git a/src/path.c b/src/path.c
index 85233d6..fd5d2a3 100644
--- a/src/path.c
+++ b/src/path.c
@@ -10,8 +10,6 @@
 #include <gsl/gsl_sf_psi.h> /* for gsl_sf_psi_n() */
 #include "scattering.h"
 
-#define N 1000
-
 static double foo(double s, double range, double theta0, int k)
 {
     return sqrt(2*pow(theta0,2)*range)*sin(M_PI*s*(k-0.5)/range)/(M_PI*(k-0.5));
@@ -32,7 +30,7 @@ double path_get_coefficient(unsigned int k, double *s, double *x, double theta0,
     return sum*pow(k-0.5,2)*pow(M_PI,2)/pow(theta0*range,2);
 }
 
-path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m)
+path *path_init(double *pos, double *dir, double T0, double range, size_t len, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m)
 {
     size_t i, j;
     double T, E, mom, theta, phi;
@@ -49,22 +47,23 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
     p->dir[2] = dir[2];
 
     p->n = n;
-
-    double *s = malloc(sizeof(double)*N);
-    double *theta1 = calloc(N,sizeof(double));
-    double *theta2 = calloc(N,sizeof(double));
-    double *x = calloc(N,sizeof(double));
-    double *y = calloc(N,sizeof(double));
-    double *z = calloc(N,sizeof(double));
-    double *beta = calloc(N,sizeof(double));
-    double *t = calloc(N,sizeof(double));
-    double *dx = calloc(N,sizeof(double));
-    double *dy = calloc(N,sizeof(double));
-    double *dz = calloc(N,sizeof(double));
+    p->len = len;
+
+    double *s = calloc(len,sizeof(double));
+    double *theta1 = calloc(len,sizeof(double));
+    double *theta2 = calloc(len,sizeof(double));
+    double *x = calloc(len,sizeof(double));
+    double *y = calloc(len,sizeof(double));
+    double *z = calloc(len,sizeof(double));
+    double *beta = calloc(len,sizeof(double));
+    double *t = calloc(len,sizeof(double));
+    double *dx = calloc(len,sizeof(double));
+    double *dy = calloc(len,sizeof(double));
+    double *dz = calloc(len,sizeof(double));
 
     dz[0] = 1.0;
-    for (i = 0; i < N; i++) {
-        s[i] = range*i/(N-1);
+    for (i = 0; i < len; i++) {
+        s[i] = range*i/(len-1);
         for (j = 0; j < n; j++) {
             theta1[i] += z1[j]*foo(s[i],range,theta0,j+1);
             theta2[i] += z2[j]*foo(s[i],range,theta0,j+1);
@@ -125,7 +124,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
 
     /* Now we rotate and translate all the positions and rotate all the
      * directions. */
-    for (i = 0; i < N; i++) {
+    for (i = 0; i < len; i++) {
         tmp[0] = x[i];
         tmp[1] = y[i];
         tmp[2] = z[i];
@@ -186,21 +185,21 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
 
 void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0)
 {
-    pos[0] = interp1d(s,p->s,p->x,N);
-    pos[1] = interp1d(s,p->s,p->y,N);
-    pos[2] = interp1d(s,p->s,p->z,N);
+    pos[0] = interp1d(s,p->s,p->x,p->len);
+    pos[1] = interp1d(s,p->s,p->y,p->len);
+    pos[2] = interp1d(s,p->s,p->z,p->len);
 
-    *beta = interp1d(s,p->s,p->beta,N);
-    *t = interp1d(s,p->s,p->t,N);
+    *beta = interp1d(s,p->s,p->beta,p->len);
+    *t = interp1d(s,p->s,p->t,p->len);
 
     /* Technically we should be interpolating between the direction vectors
      * using an algorithm like Slerp (see https://en.wikipedia.org/wiki/Slerp),
      * but since we expect the direction vectors to be pretty close to each
      * other, we just linearly interpolate and then normalize it. */
 
-    dir[0] = interp1d(s,p->s,p->dx,N);
-    dir[1] = interp1d(s,p->s,p->dy,N);
-    dir[2] = interp1d(s,p->s,p->dz,N);
+    dir[0] = interp1d(s,p->s,p->dx,p->len);
+    dir[1] = interp1d(s,p->s,p->dy,p->len);
+    dir[2] = interp1d(s,p->s,p->dz,p->len);
 
     normalize(dir);