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.

1 files changed, 9 insertions, 8 deletions

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