add path to the likelihood fit

This commit updates the likelihood fit to use the KL path expansion. Currently,
I'm just using one coefficient for the path in both x and y.

8 files changed, 82 insertions, 74 deletions

diff --git a/src/Makefile b/src/Makefile
index 1ea1685..75a75cb 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -24,7 +24,7 @@ test-charge: test-charge.o sno_charge.o misc.o
 
 test-zebra: test-zebra.o zebra.o pack2b.o
 
-fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o
+fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o
 
 clean:
 	rm -f *.o calculate_limits test Makefile.dep
diff --git a/src/fit.c b/src/fit.c
index 2a86456..7deec9d 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -20,6 +20,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
     double pos[3], dir[3];
     static size_t iter;
     double fval;
+    double z1[1], z2[1];
 
     T = x[0];
 
@@ -35,9 +36,12 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
 
     t0 = x[6];
 
-    fval = nll_muon((event *) params, T, pos, dir, t0);
+    z1[0] = x[7];
+    z2[0] = x[8];
 
-    printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f f() = %7.3e\n",
+    fval = nll_muon((event *) params, T, pos, dir, t0, z1, z2, 1);
+
+    printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %4.2f %4.2f f() = %7.3e\n",
            iter++,
            x[0],
            x[1],
@@ -46,6 +50,8 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
            x[4],
            x[5],
            x[6],
+           x[7],
+           x[8],
            fval);
 
     return fval;
@@ -53,20 +59,22 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
 
 int fit_event(event *ev, double *fmin)
 {
-    double x[7], ss[7], lb[7], ub[7];
+    double x[9], ss[9], lb[9], ub[9];
     int rv;
 
-    nlopt_opt opt = nlopt_create(NLOPT_LN_SBPLX, 7);
+    nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
     nlopt_set_min_objective(opt,nll,ev);
     nlopt_set_ftol_abs(opt, 1e-5);
 
-    x[0] = 1000.0;
-    x[1] = 100.0;
-    x[2] = 100.0;
-    x[3] = 100.0;
+    x[0] = 1200.0;
+    x[1] = 0.0;
+    x[2] = 0.0;
+    x[3] = 0.0;
     x[4] = 1.57;
     x[5] = 0.0;
     x[6] = 130.0;
+    x[7] = 0.0;
+    x[8] = 0.0;
 
     ss[0] = 10.0;
     ss[1] = 10.0;
@@ -75,6 +83,8 @@ int fit_event(event *ev, double *fmin)
     ss[4] = 0.1;
     ss[5] = 0.1;
     ss[6] = 1.0;
+    ss[7] = 0.1;
+    ss[8] = 0.1;
 
     lb[0] = 1.0;
     lb[1] = -1000.0;
@@ -83,6 +93,8 @@ int fit_event(event *ev, double *fmin)
     lb[4] = 0.0;
     lb[5] = 0.0;
     lb[6] = 0.0;
+    lb[7] = -10.0;
+    lb[8] = -10.0;
 
     ub[0] = 100000.0;
     ub[1] = 1000.0;
@@ -91,6 +103,8 @@ int fit_event(event *ev, double *fmin)
     ub[4] = M_PI;
     ub[5] = 2*M_PI;
     ub[6] = 400.0;
+    ub[7] = 10.0;
+    ub[8] = 10.0;
 
     nlopt_set_lower_bounds(opt, lb);
     nlopt_set_upper_bounds(opt, ub);
diff --git a/src/likelihood.c b/src/likelihood.c
index a4864a8..f06fade 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -11,6 +11,13 @@
 #include "optics.h"
 #include "sno_charge.h"
 #include "pdg.h"
+#include "path.h"
+#include <stddef.h> /* for size_t */
+
+typedef struct intParams {
+    path *p;
+    int i;
+} intParams;
 
 double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
 {
@@ -76,28 +83,16 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
 
 static double gsl_muon_time(double x, void *params)
 {
-    double *params2 = (double *) params;
-    double T0 = params2[0];
-    double pos0[3], dir[3], pos[3], pmt_dir[3], R, n, wavelength0;
-    int i;
-    double t;
-    i = (int) params2[1];
-    pos0[0] = params2[2];
-    pos0[1] = params2[3];
-    pos0[2] = params2[4];
-    dir[0] = params2[5];
-    dir[1] = params2[6];
-    dir[2] = params2[7];
-
-    pos[0] = pos0[0] + dir[0]*x;
-    pos[1] = pos0[1] + dir[1]*x;
-    pos[2] = pos0[2] + dir[2]*x;
+    intParams *pars = (intParams *) params;
+    double dir[3], pos[3], pmt_dir[3], R, n, wavelength0, T, t, theta0, distance;
+
+    path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
 
     R = NORM(pos);
 
-    SUB(pmt_dir,pmts[i].pos,pos);
+    SUB(pmt_dir,pmts[pars->i].pos,pos);
 
-    double distance = NORM(pmt_dir);
+    distance = NORM(pmt_dir);
 
     /* FIXME: I just calculate delta assuming 400 nm light. */
     wavelength0 = 400.0;
@@ -107,38 +102,32 @@ static double gsl_muon_time(double x, void *params)
         n = get_index_snoman_h2o(wavelength0);
     }
 
-    t = x/SPEED_OF_LIGHT + distance*n/SPEED_OF_LIGHT;
+    t += distance*n/SPEED_OF_LIGHT;
 
-    return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+    return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
 static double gsl_muon_charge(double x, void *params)
 {
-    double *params2 = (double *) params;
-    double T0 = params2[0];
-    double pos0[3], dir[3], pos[3];
-    int i;
-    i = (int) params2[1];
-    pos0[0] = params2[2];
-    pos0[1] = params2[3];
-    pos0[2] = params2[4];
-    dir[0] = params2[5];
-    dir[1] = params2[6];
-    dir[2] = params2[7];
-
-    pos[0] = pos0[0] + dir[0]*x;
-    pos[1] = pos0[1] + dir[1]*x;
-    pos[2] = pos0[2] + dir[2]*x;
-
-    return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+    intParams *pars = (intParams *) params;
+    double dir[3], pos[3], T, t, theta0;
+
+    path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+
+    return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
-double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
+double getKineticEnergy(double x, double T0)
+{
+    return get_T(T0, x, HEAVY_WATER_DENSITY);
+}
+
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n)
 {
     size_t i, j;
-    double params[8];
+    intParams params;
     double total_charge;
-    double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov;
+    double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0;
     double tmean = 0.0;
     int npmt = 0;
 
@@ -155,21 +144,24 @@ double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
     gsl_function F;
     F.params = &params;
 
-    range = get_range(T, HEAVY_WATER_DENSITY);
+    range = get_range(T0, HEAVY_WATER_DENSITY);
+
+    /* Calculate total energy */
+    E0 = T0 + MUON_MASS;
+    p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+    beta0 = p0/E0;
+
+    /* FIXME: is this formula valid for muons? */
+    theta0 = get_scattering_rms(range/2,p0,beta0,1.0,HEAVY_WATER_DENSITY)/sqrt(range/2);
+
+    params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, z1, z2, n, MUON_MASS);
 
     total_charge = 0.0;
     npmt = 0;
     for (i = 0; i < MAX_PMTS; i++) {
         if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue;
 
-        params[0] = T;
-        params[1] = i;
-        params[2] = pos[0];
-        params[3] = pos[1];
-        params[4] = pos[2];
-        params[5] = dir[0];
-        params[6] = dir[1];
-        params[7] = dir[2];
+        params.i = i;
 
         /* First, we try to compute the distance along the track where the
          * PMT is at the Cerenkov angle. The reason for this is because for
@@ -245,6 +237,8 @@ double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
         }
     }
 
+    path_free(params.p);
+
     tmean /= npmt;
 
     gsl_integration_cquad_workspace_free(w);
diff --git a/src/likelihood.h b/src/likelihood.h
index 43b0823..9267ef8 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -2,6 +2,7 @@
 #define LIKELIHOOD_H
 
 #include "event.h"
+#include <stddef.h> /* for size_t */
 
 #define MAX_PE 100
 #define CHARGE_FRACTION 60000.0
@@ -11,6 +12,6 @@
 /* Event window (ns) */
 #define GTVALID 400.0
 
-double nll_muon(event *ev, double T, double *pos, double *dir, double t0);
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n);
 
 #endif
diff --git a/src/muon.c b/src/muon.c
index 7012214..8db6237 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -216,9 +216,9 @@ double get_dEdx(double T, double rho)
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
 }
 
-double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
+double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
 {
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R, E0, p0, beta0;
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R;
 
     z = 1.0;
 
@@ -244,23 +244,13 @@ double get_expected_charge(double x, double T, double T0, double *pos, double *d
     wavelength0 = 400.0;
 
     if (R <= AV_RADIUS) {
-        rho = HEAVY_WATER_DENSITY;
         n = get_index_snoman_d2o(wavelength0);
     } else {
-        rho = WATER_DENSITY;
         n = get_index_snoman_h2o(wavelength0);
     }
 
     if (beta < 1/n) return 0;
 
-    /* Calculate total energy */
-    E0 = T0 + MUON_MASS;
-    p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
-    beta0 = p0/E0;
-
-    /* FIXME: is this formula valid for muons? */
-    theta0 = get_scattering_rms(x,p0,beta0,z,rho);
-
     /* FIXME: add angular response and scattering/absorption. */
     return omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
 }
diff --git a/src/path.c b/src/path.c
index f09ec45..dff1dc7 100644
--- a/src/path.c
+++ b/src/path.c
@@ -36,6 +36,8 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
 
     path *p = malloc(sizeof(path));
 
+    p->theta0 = theta0;
+
     p->pos[0] = pos[0];
     p->pos[1] = pos[1];
     p->pos[2] = pos[2];
@@ -114,10 +116,13 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
     return p;
 }
 
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t)
+void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0)
 {
     double normal[3], k[3], tmp[3], phi;
 
+    if (s > p->spline[0]->x[p->spline[0]->size-1])
+        s = p->spline[0]->x[p->spline[0]->size-1];
+
     tmp[0] = gsl_spline_eval(p->spline[0],s,p->acc[0]);
     tmp[1] = gsl_spline_eval(p->spline[1],s,p->acc[1]);
     tmp[2] = gsl_spline_eval(p->spline[2],s,p->acc[2]);
@@ -146,6 +151,9 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t
     normalize(tmp);
 
     rotate(dir,tmp,normal,phi);
+
+    /* FIXME: This should be the *residual* scattering RMS. */
+    *theta0 = p->theta0*sqrt(s);
 }
 
 void path_free(path *p)
diff --git a/src/path.h b/src/path.h
index 8aaa024..6e73d10 100644
--- a/src/path.h
+++ b/src/path.h
@@ -15,7 +15,7 @@ typedef struct path {
 
 double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, size_t n);
 path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, double *z1, double *z2, size_t n, double m);
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t);
+void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0);
 void path_free(path *p);
 
 #endif
diff --git a/src/test.c b/src/test.c
index 3ef6531..baf0914 100644
--- a/src/test.c
+++ b/src/test.c
@@ -342,6 +342,7 @@ int test_path(char *err)
     double pos[3], dir[3];
     double E, mom, beta;
     double s;
+    double theta0;
 
     T0 = 1.0;
     range = 1.0;
@@ -371,7 +372,7 @@ int test_path(char *err)
 
             t_expected = s/(beta*SPEED_OF_LIGHT);
 
-            path_eval(p,s,pos,dir,&T,&t);
+            path_eval(p,s,pos,dir,&T,&t,&theta0);
 
             for (k = 0; k < 3; k++) {
                 if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {