1 files changed, 18 insertions, 9 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index f426271..1bf9549 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -27,14 +27,14 @@ typedef struct intParams {
     int i;
 } intParams;
 
-particle *particle_init(int id, double T0, double rho, size_t n)
+particle *particle_init(int id, double T0, size_t n)
 {
     /* Returns a struct with arrays of the particle position and kinetic
      * energy. This struct can then be passed to particle_get_energy() to
      * interpolate the particle's kinetic energy at any point along the track.
      * For example:
      *
-     *     particle *p = particle_init(IDP_MU_MINUS, 1000.0, HEAVY_WATER_DENSITY, 100);
+     *     particle *p = particle_init(IDP_MU_MINUS, 1000.0, 100);
      *     double T = particle_get_energy(x, p);
      *
      * To compute the kinetic energy as a function of distance we need to solve
@@ -59,13 +59,15 @@ particle *particle_init(int id, double T0, double rho, size_t n)
     case IDP_E_MINUS:
     case IDP_E_PLUS:
         p->mass = ELECTRON_MASS;
-        p->range = electron_get_range(T0, rho);
+        /* We use the density of light water here since we don't have the
+         * tables for heavy water for electrons. */
+        p->range = electron_get_range(T0, WATER_DENSITY);
 
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
-            dEdx = electron_get_dEdx(p->T[i-1], rho);
+            dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
         }
 
@@ -73,26 +75,33 @@ particle *particle_init(int id, double T0, double rho, size_t n)
     case IDP_MU_MINUS:
     case IDP_MU_PLUS:
         p->mass = MUON_MASS;
-        p->range = muon_get_range(T0, rho);
+        /* We use the density of heavy water here since we only load the heavy
+         * water table for muons.
+         *
+         * FIXME: It would be nice in the future to load both tables and then
+         * load the correct table based on the media. */
+        p->range = muon_get_range(T0, HEAVY_WATER_DENSITY);
 
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
-            dEdx = muon_get_dEdx(p->T[i-1], rho);
+            dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
         }
 
         break;
     case IDP_PROTON:
         p->mass = PROTON_MASS;
-        p->range = proton_get_range(T0, rho);
+        /* We use the density of light water here since we don't have the
+         * tables for heavy water for protons. */
+        p->range = proton_get_range(T0, WATER_DENSITY);
 
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
-            dEdx = proton_get_dEdx(p->T[i-1], rho);
+            dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
         }
 
@@ -521,7 +530,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     gsl_function F;
     F.params = &params;
 
-    p = particle_init(id, T0, HEAVY_WATER_DENSITY, 10000);
+    p = particle_init(id, T0, 10000);
     range = p->range;
 
     /* Calculate total energy */