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, 51 insertions, 20 deletions

diff --git a/src/muon.c b/src/muon.c
index f2a7733..5eb1f54 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -14,6 +14,7 @@
 #include "sno.h"
 #include "scattering.h"
 #include "pmt_response.h"
+#include "misc.h"
 
 static int initialized = 0;
 
@@ -171,41 +172,69 @@ double get_range(double T, double rho)
     return gsl_spline_eval(spline_range, T, acc_range)/rho;
 }
 
-double get_T(double T0, double x, double rho)
+muon_energy *muon_init_energy(double T0, double rho, size_t n)
 {
-    /* Returns the approximate kinetic energy of a muon in water after
-     * travelling `x` cm with an initial kinetic energy `T`.
+    /* Returns a struct with arrays of the muon position and kinetic energy.
+     * This struct can then be passed to muon_get_energy() to interpolate the
+     * muon's kinetic energy at any point along the track. For example:
      *
-     * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3.
+     *     muon_energy *m = muon_init_energy(1000.0, HEAVY_WATER_DENSITY, 100);
+     *     double T = muon_get_energy(x, m);
      *
-     * Return value is in MeV.
+     * To compute the kinetic energy as a function of distance we need to solve
+     * the following integral equation:
      *
-     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
-    double a, b, range, T;
+     *     T(x) = T0 - \int_0^x dT(T(x))/dx
+     *
+     * which we solve by dividing the track up into `n` segments and then
+     * numerically computing the energy at each step. */
+    size_t i;
+    double range, dEdx;
 
-    if (!initialized) {
-        if (init()) {
-            exit(1);
-        }
+    /* Get the range of the muon. */
+    range = get_range(T0, rho);
+
+    muon_energy *m = malloc(sizeof(muon_energy));
+    m->x = malloc(sizeof(double)*n);
+    m->T = malloc(sizeof(double)*n);
+    m->n = n;
+
+    m->x[0] = 0;
+    m->T[0] = T0;
+
+    /* Now we loop over the points along the track and assume dE/dx is constant
+     * between points. */
+    for (i = 1; i < n; i++) {
+        m->x[i] = range*i/(n-1);
+        dEdx = get_dEdx(m->T[i-1], rho);
+        m->T[i] = m->T[i-1] - dEdx*(m->x[i]-m->x[i-1]);
     }
 
-    range = get_range(T0, rho);
+    return m;
+}
 
-    /* This comes from Equation 33.42 in the PDG Passage of Particles Through
-     * Matter article. */
-    b = log(1 + T0/MUON_CRITICAL_ENERGY_D2O)/range;
-    /* Now we compute the ionization energy loss from the known range and b. */
-    a = b*T0/(exp(b*range)-1.0);
+double muon_get_energy(double x, muon_energy *m)
+{
+    /* Returns the approximate kinetic energy of a muon in water after
+     * travelling `x` cm with an initial kinetic energy `T`.
+     *
+     * Return value is in MeV. */
+    double T;
 
-    /* Compute the kinetic energy after travelling a distance `x` in the
-     * continuous slowing down approximation. */
-    T = -a/b + (T0+a/b)*exp(-b*x);
+    T = interp1d(x,m->x,m->T,m->n);
 
     if (T < 0) return 0;
 
     return T;
 }
 
+void muon_free_energy(muon_energy *m)
+{
+    free(m->x);
+    free(m->T);
+    free(m);
+}
+
 double get_dEdx(double T, double rho)
 {
     /* Returns the approximate dE/dx for a muon in water with kinetic energy
@@ -222,6 +251,8 @@ double get_dEdx(double T, double rho)
         }
     }
 
+    if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];
+
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
 }