1 files changed, 43 insertions, 49 deletions

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);