add a fast likelihood function

This commit adds a fast function to calculate the expected number of PE at a
PMT without numerically integrating over the track. This calculation is *much*
faster than integrating over the track (~30 ms compared to several seconds) and
so we use it during the "quick" minimization phase of the fit to quickly find
the best position.

1 files changed, 64 insertions, 2 deletions

diff --git a/src/scattering.c b/src/scattering.c
index 66e8398..433387d 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -10,6 +10,9 @@
 #include <stdlib.h> /* for exit() */
 #include <stdio.h> /* for fprintf() */
 #include <gsl/gsl_errno.h> /* for gsl_strerror() */
+#include "pdg.h"
+#include <gsl/gsl_interp.h>
+#include <gsl/gsl_spline.h>
 
 static double xlo = -1.0;
 static double xhi = 1.0;
@@ -26,6 +29,16 @@ static gsl_spline2d *spline;
 static gsl_interp_accel *xacc;
 static gsl_interp_accel *yacc;
 
+static double x2lo = 0.0;
+static double x2hi = 1.0;
+static size_t nx2 = 1000;
+
+static double *x2;
+static double *y2;
+
+static gsl_spline *spline2;
+static gsl_interp_accel *xacc2;
+
 static double prob_scatter(double wavelength, void *params)
 {
     /* Calculate the number of photons emitted per unit solid angle per cm at
@@ -45,6 +58,21 @@ static double prob_scatter(double wavelength, void *params)
     return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI);
 }
 
+static double prob_scatter2(double wavelength, void *params)
+{
+    /* Calculate the number of photons emitted per per cm for a particle
+     * travelling with velocity `beta`. */
+    double qe, index;
+
+    double beta = ((double *) params)[0];
+
+    qe = get_quantum_efficiency(wavelength);
+
+    index = get_index_snoman_d2o(wavelength);
+
+    return 2*M_PI*FINE_STRUCTURE_CONSTANT*(1-(1/(beta*beta*index*index)))*qe/pow(wavelength,2)*1e7;
+}
+
 void init_interpolation(void)
 {
     size_t i, j;
@@ -52,6 +80,8 @@ void init_interpolation(void)
     double result, error;
     size_t nevals;
     int status;
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
 
     x = malloc(nx*sizeof(double));
     y = malloc(ny*sizeof(double));
@@ -69,9 +99,8 @@ void init_interpolation(void)
         y[i] = ylo + (yhi-ylo)*i/(ny-1);
     }
 
-    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100);
+    w = gsl_integration_cquad_workspace_alloc(100);
 
-    gsl_function F;
     F.function = &prob_scatter;
     F.params = params;
 
@@ -92,6 +121,34 @@ void init_interpolation(void)
 
     gsl_spline2d_init(spline, x, y, z, nx, ny);
 
+    x2 = malloc(nx2*sizeof(double));
+    y2 = malloc(nx2*sizeof(double));
+
+    spline2 = gsl_spline_alloc(gsl_interp_linear, nx2);
+    xacc2 = gsl_interp_accel_alloc();
+
+    for (i = 0; i < nx2; i++) {
+        x2[i] = x2lo + (x2hi-x2lo)*i/(nx2-1);
+    }
+
+    F.function = &prob_scatter2;
+    F.params = params;
+
+    for (i = 0; i < nx2; i++) {
+        params[0] = x2[i];
+
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+        if (status) {
+            fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
+            exit(1);
+        }
+
+        y2[i] = result;
+    }
+
+    gsl_spline_init(spline2, x2, y2, nx2);
+
     gsl_integration_cquad_workspace_free(w);
 }
 
@@ -106,6 +163,11 @@ double get_probability(double beta, double cos_theta, double theta0)
     return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta);
 }
 
+double get_probability2(double beta)
+{
+    return gsl_spline_eval(spline2, beta, xacc2);
+}
+
 void free_interpolation(void)
 {
     free(x);