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