extend electron range tables up to 1 TeV

4 files changed, 259 insertions, 1 deletions

diff --git a/src/.gitignore b/src/.gitignore
index 47de265..209802f 100644
--- a/src/.gitignore
+++ b/src/.gitignore
@@ -10,3 +10,4 @@ release.h
 calculate-csda-range
 test-time-pdf
 zdab-cat
+extend-electron-range-table
diff --git a/src/Makefile b/src/Makefile
index cb701cb..834bdb4 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -7,7 +7,7 @@ PREFIX?=$(HOME)/local
 INSTALL_BIN=$(PREFIX)/bin
 INSTALL=install
 
-all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks zdab-cat
+all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf test-zebra test-find-peaks zdab-cat extend-electron-range-table
 
 Makefile.dep:
 	-$(CC) -MM *.c > Makefile.dep
@@ -40,6 +40,8 @@ calculate-csda-range: calculate-csda-range.o
 
 zdab-cat: zdab-cat.o zebra.o pmt.o vector.o misc.o zdab_utils.o pack2b.o db.o dqxx.o dict.o siphash.o release.o dc.o sort.o util.o hdf5_utils.o sno_charge.o
 
+extend-electron-range-table: extend-electron-range-table.o misc.o
+
 install:
 	@mkdir -p $(INSTALL_BIN)
 	$(INSTALL) fit $(INSTALL_BIN)
diff --git a/src/e_water_liquid.txt b/src/e_water_liquid.txt
index a4e3a8f..633440f 100644
--- a/src/e_water_liquid.txt
+++ b/src/e_water_liquid.txt
@@ -100,3 +100,23 @@ MeV       MeV cm2/g MeV cm2/g MeV cm2/g g/cm2
 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
+# Note: The values below were extrapolated by me from the actual ESTAR data above. They are most likely not 100% correct!
+2.000E+04 -         5.492E+02 5.518E+02 2.071E+02 -         -
+3.000E+04 -         8.242E+02 8.268E+02 2.218E+02 -         -
+4.000E+04 -         1.099E+03 1.102E+03 2.323E+02 -         -
+5.000E+04 -         1.374E+03 1.377E+03 2.404E+02 -         -
+6.000E+04 -         1.649E+03 1.652E+03 2.470E+02 -         -
+7.000E+04 -         1.924E+03 1.927E+03 2.526E+02 -         -
+8.000E+04 -         2.199E+03 2.202E+03 2.574E+02 -         -
+9.000E+04 -         2.474E+03 2.477E+03 2.617E+02 -         -
+1.000E+05 -         2.749E+03 2.752E+03 2.656E+02 -         -
+2.000E+05 -         5.499E+03 5.502E+03 2.907E+02 -         -
+3.000E+05 -         8.249E+03 8.252E+03 3.055E+02 -         -
+4.000E+05 -         1.100E+04 1.100E+04 3.159E+02 -         -
+5.000E+05 -         1.375E+04 1.375E+04 3.241E+02 -         -
+6.000E+05 -         1.650E+04 1.650E+04 3.307E+02 -         -
+7.000E+05 -         1.925E+04 1.925E+04 3.363E+02 -         -
+8.000E+05 -         2.200E+04 2.200E+04 3.412E+02 -         -
+9.000E+05 -         2.475E+04 2.475E+04 3.454E+02 -         -
+1.000E+06 -         2.750E+04 2.750E+04 3.493E+02 -         -
+
diff --git a/src/extend-electron-range-table.c b/src/extend-electron-range-table.c
new file mode 100644
index 0000000..89901df
--- /dev/null
+++ b/src/extend-electron-range-table.c
@@ -0,0 +1,235 @@
+/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
+ *
+ * This program is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option)
+ * any later version.
+
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
+ * more details.
+
+ * You should have received a copy of the GNU General Public License along with
+ * this program. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#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"
+#include "misc.h"
+
+static double *x, *dEdx_rad, *dEdx;
+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 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);
+    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;
+            }
+            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);
+
+    return 0;
+
+err:
+    fclose(f);
+
+    return -1;
+}
+
+/* Returns the approximate dE/dx for a electron in water with kinetic energy
+ * `T`. If `T` is outside the bounds of the range table, extrapolate linearly
+ * using the last two points.
+ *
+ * `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. */
+double electron_get_dEdx(double T, double rho)
+{
+    if (T < spline_dEdx->x[0]) {
+        return spline_dEdx->y[0];
+    } else if (T < spline_dEdx->x[size-1]) {
+        return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
+    }
+
+    /* We extrapolate using the last two points. */
+    return spline_dEdx->y[size-1] + (spline_dEdx->y[size-1]-spline_dEdx->y[size-2])*(T-spline_dEdx->x[size-1])/(spline_dEdx->x[size-1]-spline_dEdx->x[size-2]);
+}
+
+/* Returns the approximate radiative dE/dx for a electron in water with kinetic
+ * energy `T`. If `T` is outside the bounds of the range table, extrapolate
+ * linearly using the last two points.
+ *
+ * `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. */
+double electron_get_dEdx_rad(double T, double rho)
+{
+    if (T < spline_dEdx_rad->x[0]) {
+        return spline_dEdx_rad->y[0];
+    } else if (T < spline_dEdx_rad->x[size-1]) {
+        return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
+    }
+
+    /* We extrapolate using the last two points. */
+    return spline_dEdx_rad->y[size-1] + (spline_dEdx_rad->y[size-1]-spline_dEdx_rad->y[size-2])*(T-spline_dEdx_rad->x[size-1])/(spline_dEdx_rad->x[size-1]-spline_dEdx_rad->x[size-2]);
+}
+
+double electron_get_range(double T0)
+{
+    double T, dx, range;
+
+    T = T0;
+
+    dx = 1e-5;
+
+    range = 0.0;
+    while (T > 0) {
+        double dEdx2 = electron_get_dEdx(T, WATER_DENSITY);
+        T -= dEdx2*dx;
+        range += dx;
+    }
+
+    return range;
+}
+
+int main(int argc, char **argv)
+{
+    int i;
+    double T0;
+    double Ts[18] = {
+        2.000E+04,
+        3.000E+04,
+        4.000E+04,
+        5.000E+04,
+        6.000E+04,
+        7.000E+04,
+        8.000E+04,
+        9.000E+04,
+        1.000E+05,
+        2.000E+05,
+        3.000E+05,
+        4.000E+05,
+        5.000E+05,
+        6.000E+05,
+        7.000E+05,
+        8.000E+05,
+        9.000E+05,
+        1.000E+06};
+
+    init();
+
+    for (i = 0; i < LEN(Ts); i++) {
+        T0 = Ts[i];
+
+        printf("%.3E - %.3E %.3E %.3E - -\n", T0, electron_get_dEdx_rad(T0,WATER_DENSITY), electron_get_dEdx(T0,WATER_DENSITY), electron_get_range(T0));
+    }
+
+    return 0;
+}