update muon kinetic energy calculation

This commit updates the calculation of the muon kinetic energy as a function of
distance along the track. Previously I was using an approximation from the PDG,
but it doesn't seem to be very accurate and won't generalize to the case of
electrons. The kinetic energy is now calculated using the tabulated values of
dE/dx as a function of energy.

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;