1 files changed, 27 insertions, 19 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 76364fd..6881e8f 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -122,7 +122,7 @@ static double gsl_muon_charge(double x, void *params)
     return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
 }
 
-double get_total_charge_approx(double T0, double *pos, double *dir, int i, double smax, double theta0, double *t)
+double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t)
 {
     /* Returns the approximate expected number of photons seen by PMT `i` using
      * an analytic formula.
@@ -218,7 +218,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
     b = (s-z)/a;
 
     /* Compute the kinetic energy at the point `s` along the track. */
-    T = get_T(T0,s,HEAVY_WATER_DENSITY);
+    T = muon_get_energy(s,m);
 
     /* Calculate the particle velocity at the point `s`. */
     E = T + MUON_MASS;
@@ -258,36 +258,36 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl
     return f*exp(-a/absorption_length)*n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI);
 }
 
-double getKineticEnergy(double x, double T0)
-{
-    return get_T(T0, x, HEAVY_WATER_DENSITY);
-}
+typedef struct betaRootParams {
+    muon_energy *m;
+    double beta_min;
+} betaRootParams;
 
 static double beta_root(double x, void *params)
 {
     /* Function used to find at what point along a track a particle crosses
      * some threshold in beta. */
-    double T, E, p, beta, T0, beta_min;
+    double T, E, p, beta;
+    betaRootParams *pars;
 
-    T0 = ((double *) params)[0];
-    beta_min = ((double *) params)[1];
+    pars = (betaRootParams *) params;
 
-    T = get_T(T0, x, HEAVY_WATER_DENSITY);
+    T = muon_get_energy(x, pars->m);
 
     /* Calculate total energy */
     E = T + MUON_MASS;
     p = sqrt(E*E - MUON_MASS*MUON_MASS);
     beta = p/E;
 
-    return beta - beta_min;
+    return beta - pars->beta_min;
 }
 
-static int get_smax(double T0, double beta_min, double range, double *smax)
+static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
 {
     /* Find the point along the track at which the particle's velocity drops to
      * `beta_min`. */
     int status;
-    double params[2];
+    betaRootParams pars;
     gsl_root_fsolver *s;
     gsl_function F;
     int iter = 0, max_iter = 100;
@@ -295,11 +295,11 @@ static int get_smax(double T0, double beta_min, double range, double *smax)
 
     s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
 
-    params[0] = T0;
-    params[1] = beta_min;
+    pars.m = m;
+    pars.beta_min = beta_min;
 
     F.function = &beta_root;
-    F.params = params;
+    F.params = &pars;
 
     gsl_root_fsolver_set(s, &F, 0.0, range);
 
@@ -323,6 +323,11 @@ static int get_smax(double T0, double beta_min, double range, double *smax)
     return status;
 }
 
+double getKineticEnergy(double x, void *params)
+{
+    return muon_get_energy(x,(muon_energy *) params);
+}
+
 double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
 {
     size_t i, j, nhit;
@@ -330,6 +335,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     double total_charge;
     double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
     double tmean = 0.0;
+    muon_energy *m;
 
     double mu_direct[MAX_PMTS];
     double ts[MAX_PMTS];
@@ -354,10 +360,12 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
     /* 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);
+    m = muon_init_energy(T0,HEAVY_WATER_DENSITY,10000);
+
+    params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, m, z1, z2, n, MUON_MASS);
 
     if (beta0 > BETA_MIN)
-        get_smax(T0, BETA_MIN, range, &smax);
+        get_smax(m, BETA_MIN, range, &smax);
     else
         smax = 0.0;
 
@@ -368,7 +376,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
         params.i = i;
 
         if (fast) {
-            mu_direct[i] = get_total_charge_approx(T0, pos, dir, i, smax, theta0, &ts[i]);
+            mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i]);
             ts[i] += t0;
         } else {
             F.function = &gsl_muon_charge;