hardcode the density when computing dE/dx

Since we only have the range and dE/dx tables for light water for electrons and
protons it's not correct to use the heavy water density. Also, even though we
have both tables for muons, currently we only load the heavy water table, so we
hardcode the density to that of heavy water.

In the future, it would be nice to load both tables and use the correct one
depending on if we are fitting in the heavy or light water.

3 files changed, 22 insertions, 13 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 */
diff --git a/src/likelihood.h b/src/likelihood.h
index 945fc46..71f77f0 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -41,7 +41,7 @@ typedef struct particle {
     size_t n;
 } particle;
 
-particle *particle_init(int id, double T0, double rho, size_t n);
+particle *particle_init(int id, double T0, size_t n);
 double particle_get_energy(double x, particle *p);
 void particle_free(particle *p);
 double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected);
diff --git a/src/test.c b/src/test.c
index 20a66ec..c1da272 100644
--- a/src/test.c
+++ b/src/test.c
@@ -124,7 +124,7 @@ int test_proton_get_energy(char *err)
 
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
-    p = particle_init(IDP_PROTON, T0,1.0,10000);
+    p = particle_init(IDP_PROTON, T0, 10000);
     T = particle_get_energy(1e-9,p);
 
     /* At the beginning of the track we should have roughly the same energy. */
@@ -160,7 +160,7 @@ int test_electron_get_energy(char *err)
 
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
-    p = particle_init(IDP_E_MINUS, T0,1.0,10000);
+    p = particle_init(IDP_E_MINUS, T0, 10000);
     T = particle_get_energy(1e-9,p);
 
     /* At the beginning of the track we should have roughly the same energy. */
@@ -196,7 +196,7 @@ int test_muon_get_energy(char *err)
 
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
-    p = particle_init(IDP_MU_MINUS, T0,1.0,10000);
+    p = particle_init(IDP_MU_MINUS, T0, 10000);
     T = particle_get_energy(1e-9,p);
 
     /* At the beginning of the track we should have roughly the same energy. */