set the maximum kinetic energy in the fit dynamically based on particle ID

The range and energy loss tables have different maximum values for electrons,
muons, and protons so we have to dynamically set the maximum energy of the fit
in order to avoid a GSL interpolation error.

This commit adds {electron,muon,proton}_get_max_energy() functions to return
the maximum energy in the tables and that is then used to set the maximum value
in the fit.

7 files changed, 57 insertions, 8 deletions

diff --git a/src/electron.c b/src/electron.c
index 5b4133a..3f5d4d8 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -365,6 +365,21 @@ err:
     return -1;
 }
 
+/* Returns the maximum kinetic energy for an electron in the range tables.
+ *
+ * If you call electron_get_range() or electron_get_dEdx() with a kinetic
+ * energy higher you will get a GSL interpolation error. */
+double electron_get_max_energy(void)
+{
+    if (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    return x[size-1];
+}
+
 double electron_get_range(double T, double rho)
 {
     /* Returns the approximate range a electron with kinetic energy `T` will travel
diff --git a/src/electron.h b/src/electron.h
index cc550a7..01a8915 100644
--- a/src/electron.h
+++ b/src/electron.h
@@ -19,6 +19,7 @@
 
 #include <stddef.h> /* for size_t */
 
+double electron_get_max_energy(void);
 double electron_get_shower_photons(double T0, double rad);
 void electron_get_position_distribution_parameters(double T0, double *a, double *b);
 double electron_get_angular_distribution_alpha(double T0);
diff --git a/src/fit.c b/src/fit.c
index 84ebbdc..fd62c45 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -44,16 +44,14 @@
 #include "find_peaks.h"
 #include "util.h"
 #include "dc.h"
+#include "electron.h"
+#include "muon.h"
+#include "proton.h"
 
 /* Maximum number of fit parameters. Should be at least 4 + 3*MAX_VERTICES. */
 #define MAX_PARS 100
 /* Maximum number of peaks to search for in Hough transform. */
 #define MAX_NPEAKS 10
-/* Maximum kinetic energy for any particle. Currently set to 10^4 GeV which is
- * approximately the energy of the most energetic muons in SNO (see the paper
- * "Measurement of the Cosmic Ray and Neutrino-Induced Muon Flux at the Sudbury
- * Neutrino Observatory".*/
-#define MAX_ENERGY 1e7
 
 char *GitSHA1(void);
 char *GitDirty(void);
@@ -5263,7 +5261,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
      * press ctrl-c. */
     size_t i, j;
     fitParams fpars;
-    double x[100], ss[100], lb[100], ub[100], fval, n_d2o, x0[100], T0, Tmin, mass, pos[3], t0, dir[3];
+    double x[100], ss[100], lb[100], ub[100], fval, n_d2o, x0[100], T0, Tmin, mass, pos[3], t0, dir[3], Tmax;
     struct timeval tv_start, tv_stop;
     double time_elapsed;
     double peak_theta[MAX_NPEAKS];
@@ -5371,13 +5369,16 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
             case IDP_E_MINUS:
             case IDP_E_PLUS:
                 mass = ELECTRON_MASS;
+                Tmax = electron_get_max_energy();
                 break;
             case IDP_MU_MINUS:
             case IDP_MU_PLUS:
                 mass = MUON_MASS;
+                Tmax = muon_get_max_energy();
                 break;
             case IDP_PROTON:
                 mass = PROTON_MASS;
+                Tmax = proton_get_max_energy();
                 break;
             default:
                 fprintf(stderr, "unknown particle id %i\n", id[j]);
@@ -5391,7 +5392,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
              * threshold. */
             double Tmin2 = sqrt(mass*mass/(1-pow(0.9,2)))-mass;
             if (T0 < Tmin2) T0 = Tmin2;
-            if (T0 > MAX_ENERGY) T0 = MAX_ENERGY;
+            if (T0 > Tmax) T0 = Tmax;
 
             x[4+3*j] = T0;
             ss[4+3*j] = x[4+3*j]*0.1;
@@ -5400,7 +5401,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double
             lb[4+3*j] = Tmin;
             lb[5+3*j] = -INFINITY;
             lb[6+3*j] = -INFINITY;
-            ub[4+3*j] = MAX_ENERGY;
+            ub[4+3*j] = Tmax;
             ub[5+3*j] = +INFINITY;
             ub[6+3*j] = +INFINITY;
         }
diff --git a/src/muon.c b/src/muon.c
index aea2309..9d6c91a 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -312,6 +312,21 @@ err:
     return -1;
 }
 
+/* Returns the maximum kinetic energy for a muon in the range tables.
+ *
+ * If you call muon_get_range() or muon_get_dEdx() with a kinetic energy higher
+ * you will get a GSL interpolation error. */
+double muon_get_max_energy(void)
+{
+    if (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    return x[size-1];
+}
+
 double muon_get_range(double T, double rho)
 {
     /* Returns the approximate range a muon with kinetic energy `T` will travel
diff --git a/src/muon.h b/src/muon.h
index a8d2ea1..8e1d39e 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -30,6 +30,7 @@
  * FIXME: Actually determine what this is. */
 #define MUON_PHOTONS_PER_MEV 7368.0
 
+double muon_get_max_energy(void);
 double muon_get_shower_photons(double T0, double rad);
 void muon_get_position_distribution_parameters(double T0, double *a, double *b);
 double muon_get_angular_distribution_alpha(double T0);
diff --git a/src/proton.c b/src/proton.c
index 664d930..4714d76 100644
--- a/src/proton.c
+++ b/src/proton.c
@@ -165,6 +165,21 @@ err:
     return -1;
 }
 
+/* Returns the maximum kinetic energy for a proton in the range tables.
+ *
+ * If you call proton_get_range() or proton_get_dEdx() with a kinetic
+ * energy higher you will get a GSL interpolation error. */
+double proton_get_max_energy(void)
+{
+    if (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    return x[size-1];
+}
+
 double proton_get_range(double T, double rho)
 {
     /* Returns the approximate range a proton with kinetic energy `T` will travel
diff --git a/src/proton.h b/src/proton.h
index f658c09..c402229 100644
--- a/src/proton.h
+++ b/src/proton.h
@@ -28,6 +28,7 @@
  * FIXME: Actually determine what this is. */
 #define PROTON_PHOTONS_PER_MEV 400.0
 
+double proton_get_max_energy(void);
 double proton_get_range(double T, double rho);
 double proton_get_dEdx_rad(double T, double rho);
 double proton_get_dEdx(double T, double rho);