update fit to fit for electrons and protons

1 files changed, 180 insertions, 0 deletions

diff --git a/src/proton.c b/src/proton.c
new file mode 100644
index 0000000..922c552
--- /dev/null
+++ b/src/proton.c
@@ -0,0 +1,180 @@
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+#include <math.h>
+#include "optics.h"
+#include "quantum_efficiency.h"
+#include "solid_angle.h"
+#include "pdg.h"
+#include "vector.h"
+#include "proton.h"
+#include "sno.h"
+#include "scattering.h"
+#include "pmt_response.h"
+#include "misc.h"
+
+static int initialized = 0;
+
+static double *x, *dEdx, *csda_range;
+static size_t size;
+
+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("proton_water_liquid.txt", "r");
+
+    if (!f) {
+        fprintf(stderr, "failed to open proton_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 = 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 3:
+                dEdx[n] = value;
+                break;
+            case 4:
+                csda_range[n] = value;
+                break;
+            }
+            j += 1;
+            str = strtok(NULL," \n");
+        }
+
+        n += 1;
+    }
+
+    fclose(f);
+
+    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);
+
+    initialized = 1;
+
+    return 0;
+
+err:
+    fclose(f);
+
+    return -1;
+}
+
+double proton_get_range(double T, double rho)
+{
+    /* Returns the approximate range a proton with kinetic energy `T` will travel
+     * in water before losing all of its energy. This range is interpolated
+     * based on data from the PDG which uses the continuous slowing down
+     * approximation.
+     *
+     * `T` should be in MeV, and `rho` should be in g/cm^3.
+     *
+     * Return value is in cm.
+     *
+     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+    if (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    return gsl_spline_eval(spline_range, T, acc_range)/rho;
+}
+
+double proton_get_dEdx(double T, double rho)
+{
+    /* Returns the approximate dE/dx for a proton 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 (!initialized) {
+        if (init()) {
+            exit(1);
+        }
+    }
+
+    if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];
+
+    return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
+}