update electron range tables

This commit also adds a script to calculate the CSDA range for electrons from a
table of the stopping power as a function of energy. We need this script since
the NIST ESTAR website will only output the stopping power above a certain
energy threshold (1 GeV for electrons).

See https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html.

4 files changed, 193 insertions, 1 deletions

diff --git a/src/.gitignore b/src/.gitignore
new file mode 100644
index 0000000..2750fd2
--- /dev/null
+++ b/src/.gitignore
@@ -0,0 +1 @@
+calculate-csda-range
diff --git a/src/Makefile b/src/Makefile
index 62a95af..80eef49 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -3,7 +3,7 @@ release_hdr := $(shell sh -c './mkreleasehdr.sh')
 CFLAGS=-O2 -Wall -g -DSWAP_BYTES
 LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++
 
-all: test test-vector test-likelihood fit test-charge test-path
+all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range
 
 Makefile.dep:
 	-$(CC) -MM *.c > Makefile.dep
@@ -28,6 +28,8 @@ test-zebra: test-zebra.o zebra.o pack2b.o
 
 fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o
 
+calculate-csda-range: calculate-csda-range.o
+
 clean:
 	rm -f *.o calculate_limits test Makefile.dep
 
diff --git a/src/calculate-csda-range.c b/src/calculate-csda-range.c
new file mode 100644
index 0000000..7d71096
--- /dev/null
+++ b/src/calculate-csda-range.c
@@ -0,0 +1,176 @@
+#include <stdio.h>
+#include <gsl/gsl_spline.h>
+#include <stdlib.h> /* for size_t, strtod() */
+#include <errno.h> /* for errno */
+#include <string.h> /* for strerror(), strtok() */
+#include "sno.h"
+
+static double *x, *dEdx_rad, *dEdx, *csda_range;
+static size_t size;
+
+static gsl_interp_accel *acc_dEdx_rad;
+static gsl_spline *spline_dEdx_rad;
+
+static gsl_interp_accel *acc_dEdx;
+static gsl_spline *spline_dEdx;
+
+static gsl_interp_accel *acc_range;
+static gsl_spline *spline_range;
+
+static int init()
+{
+    int i, j;
+    char line[256];
+    char *str;
+    double value;
+    int n;
+
+    FILE *f = fopen("e_water_liquid.txt", "r");
+
+    if (!f) {
+        fprintf(stderr, "failed to open e_water_liquid.txt: %s\n", strerror(errno));
+        return -1;
+    }
+
+    i = 0;
+    n = 0;
+    /* For the first pass, we just count how many values there are. */
+    while (fgets(line, sizeof(line), f)) {
+        size_t len = strlen(line);
+        if (len && (line[len-1] != '\n')) {
+            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
+            goto err;
+        }
+
+        i += 1;
+
+        /* Skip the first 8 lines since it's just a header. */
+        if (i <= 8) continue;
+
+        if (!len) continue;
+        else if (line[0] == '#') continue;
+
+        str = strtok(line," \n");
+
+        while (str) {
+            value = strtod(str, NULL);
+            str = strtok(NULL," \n");
+        }
+
+        n += 1;
+    }
+
+    x = malloc(sizeof(double)*n);
+    dEdx_rad = malloc(sizeof(double)*n);
+    dEdx = malloc(sizeof(double)*n);
+    csda_range = malloc(sizeof(double)*n);
+    size = n;
+
+    i = 0;
+    n = 0;
+    /* Now, we actually store the values. */
+    rewind(f);
+    while (fgets(line, sizeof(line), f)) {
+        size_t len = strlen(line);
+        if (len && (line[len-1] != '\n')) {
+            fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line);
+            goto err;
+        }
+
+        i += 1;
+
+        /* Skip the first 8 lines since it's just a header. */
+        if (i <= 8) continue;
+
+        if (!len) continue;
+        else if (line[0] == '#') continue;
+
+        str = strtok(line," \n");
+
+        j = 0;
+        while (str) {
+            value = strtod(str, NULL);
+            switch (j) {
+            case 0:
+                x[n] = value;
+                break;
+            case 2:
+                dEdx_rad[n] = value;
+                break;
+            case 3:
+                dEdx[n] = value;
+                break;
+            case 4:
+                csda_range[n] = value;
+                break;
+            }
+            j += 1;
+            str = strtok(NULL," \n");
+        }
+
+        n += 1;
+    }
+
+    fclose(f);
+
+    acc_dEdx_rad = gsl_interp_accel_alloc();
+    spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
+    gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);
+
+    acc_dEdx = gsl_interp_accel_alloc();
+    spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
+    gsl_spline_init(spline_dEdx, x, dEdx, size);
+
+    acc_range = gsl_interp_accel_alloc();
+    spline_range = gsl_spline_alloc(gsl_interp_linear, size);
+    gsl_spline_init(spline_range, x, csda_range, size);
+
+    return 0;
+
+err:
+    fclose(f);
+
+    return -1;
+}
+
+double electron_get_dEdx(double T, double rho)
+{
+    /* Returns the approximate dE/dx for a electron in water with kinetic energy
+     * `T`.
+     *
+     * `T` should be in MeV and `rho` in g/cm^3.
+     *
+     * Return value is in MeV/cm.
+     *
+     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+    if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];
+
+    return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
+}
+
+int main(int argc, char **argv)
+{
+    int i;
+    double T0, T, dx, range;
+
+    init();
+
+    dx = 1e-5;
+
+    for (i = 0; i < size; i++) {
+        T0 = x[i];
+
+        T = T0;
+
+        range = 0.0;
+        while (T > 0) {
+            double dEdx2 = electron_get_dEdx(T, WATER_DENSITY);
+            T -= dEdx2*dx;
+            range += dx;
+        }
+
+        printf("%.3E %.3E\n", T0, range);
+    }
+
+    return 0;
+}
diff --git a/src/e_water_liquid.txt b/src/e_water_liquid.txt
index 4e29104..a4e3a8f 100644
--- a/src/e_water_liquid.txt
+++ b/src/e_water_liquid.txt
@@ -87,3 +87,16 @@ MeV       MeV cm2/g MeV cm2/g MeV cm2/g g/cm2
 8.000E+02 2.381E+00 2.133E+01 2.371E+01 9.258E+01 7.425E-01 1.121E+01 
 9.000E+02 2.391E+00 2.406E+01 2.645E+01 9.657E+01 7.605E-01 1.145E+01 
 1.000E+03 2.400E+00 2.679E+01 2.919E+01 1.002E+02 7.759E-01 1.166E+01 
+1.500E+03 2.435E+00 4.047E+01 4.291E+01 1.142E+02 -         1.247E+01
+2.000E+03 2.459E+00 5.418E+01 5.664E+01 1.244E+02 -         1.304E+01
+2.500E+03 2.479E+00 6.790E+01 7.038E+01 1.323E+02 -         1.349E+01
+3.000E+03 2.494E+00 8.163E+01 8.413E+01 1.388E+02 -         1.385E+01
+3.500E+03 2.507E+00 9.537E+01 9.788E+01 1.443E+02 -         1.416E+01
+4.000E+03 2.519E+00 1.091E+02 1.116E+02 1.490E+02 -         1.443E+01
+4.500E+03 2.529E+00 1.229E+02 1.254E+02 1.533E+02 -         1.466E+01
+5.000E+03 2.538E+00 1.366E+02 1.391E+02 1.570E+02 -         1.488E+01
+6.000E+03 2.553E+00 1.641E+02 1.667E+02 1.636E+02 -         1.524E+01
+7.000E+03 2.566E+00 1.916E+02 1.942E+02 1.692E+02 -         1.555E+01
+8.000E+03 2.578E+00 2.192E+02 2.217E+02 1.740E+02 -         1.582E+01
+9.000E+03 2.588E+00 2.467E+02 2.493E+02 1.782E+02 -         1.605E+01
+1.000E+04 2.597E+00 2.742E+02 2.768E+02 1.820E+02 -         1.626E+01