src/electron.c

/* 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"

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;
}
Age	Commit message (Collapse)	Author
2019-06-14	set the maximum kinetic energy in the fit dynamically based on particle ID	tlatorre
	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.
2019-03-26	add energy dependent number of shower photons	tlatorre
	This commit updates the code to calculate the number of Cerenkov photons from secondary particles produced in an electromagnetic shower from electrons to use an energy dependent formula I fit to data simulated with RAT-PAC.
2019-03-16	add GPLv3 license	tlatorre

2019-03-07	update code to allow you to run the fit outside of the src directory	tlatorre
	To enable the fitter to run outside of the src directory, I created a new function open_file() which works exactly like fopen() except that it searches for the file in both the current working directory and the path specified by an environment variable.
2019-01-31	small updates to make sure we don't calculate nans	tlatorre

2019-01-29	add a function to compute the angular distribution normalization	tlatorre
	This seems to speed things up a little bit.
2019-01-27	add photons from delta rays to likelihood calculation	tlatorre
	This commit updates the likelihood function to take into account Cerenkov light produced from delta rays produced by muons. The angular distribution of this light is currently assumed to be constant along the track and parameterized in the same way as the Cerenkov light from an electromagnetic shower. Currently I assume the light is produced uniformly along the track which isn't exactly correct, but should be good enough.
2018-11-11	update likelihood function to fit electrons!	tlatorre
	To characterize the angular distribution of photons from an electromagnetic shower I came up with the following functional form: f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) and fit this to data simulated using RAT-PAC at several different energies. I then fit the alpha and beta coefficients as a function of energy to the functional form: alpha = c0 + c1/log(c2T0 + c3) beta = c0 + c1/log(c2T0 + c3). where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3 are parameters which I fit. The longitudinal distribution of the photons generated from an electromagnetic shower is described by a gamma distribution: f(x) = x*(a-1)exp(-x/b)/(Gamma(a)b*a). This parameterization comes from the PDG "Passage of particles through matter" section 32.5. I also fit the data from my RAT-PAC simulation, but currently I am not using it, and instead using a simpler form to calculate the coefficients from the PDG (although I estimated the b parameter from the RAT-PAC data). I also sped up the calculation of the solid angle by making a lookup table since it was taking a significant fraction of the time to compute the likelihood function.
2018-10-18	update fit to fit for electrons and protons	tlatorre