make sure that the kinetic energy is zero at the last step

Occasionally when fitting electrons the kinetic energy at the last step would
be high enough that the electron never crossed the BETA_MIN threshold which
would cause the gsl routine to throw an error.

This commit updates particle_init() to set the kinetic energy at the last
step to zero to make sure that we can bisect the point along the track where
the speed drops to BETA_MIN.

1 files changed, 24 insertions, 0 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 1bf9549..b06ebc4 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -70,6 +70,14 @@ particle *particle_init(int id, double T0, size_t n)
             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]);
         }
+        /* Make sure that the energy is zero at the last step. This is so that
+         * when we try to bisect the point along the track where the speed of
+         * the particle is equal to BETA_MIN, we guarantee that there is a
+         * point along the track where the speed drops below BETA_MIN.
+         *
+         * A possible future improvement would be to dynamically compute the
+         * range here using the dE/dx table instead of reading in the range. */
+        p->T[n-1] = 0;
 
         break;
     case IDP_MU_MINUS:
@@ -89,6 +97,14 @@ particle *particle_init(int id, double T0, size_t n)
             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]);
         }
+        /* Make sure that the energy is zero at the last step. This is so that
+         * when we try to bisect the point along the track where the speed of
+         * the particle is equal to BETA_MIN, we guarantee that there is a
+         * point along the track where the speed drops below BETA_MIN.
+         *
+         * A possible future improvement would be to dynamically compute the
+         * range here using the dE/dx table instead of reading in the range. */
+        p->T[n-1] = 0;
 
         break;
     case IDP_PROTON:
@@ -104,6 +120,14 @@ particle *particle_init(int id, double T0, size_t n)
             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]);
         }
+        /* Make sure that the energy is zero at the last step. This is so that
+         * when we try to bisect the point along the track where the speed of
+         * the particle is equal to BETA_MIN, we guarantee that there is a
+         * point along the track where the speed drops below BETA_MIN.
+         *
+         * A possible future improvement would be to dynamically compute the
+         * range here using the dE/dx table instead of reading in the range. */
+        p->T[n-1] = 0;
 
         break;
     default: