path: root/src/muon.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 <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 "muon.h"
#include "sno.h"
#include "scattering.h"
#include "pmt_response.h"
#include "misc.h"
#include "util.h"

static int initialized = 0;

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;

/* Muon critical energy in H2O and D2O. These values are used in computing the
 * kinetic energy of the muon as a function of distance in get_T().
 *
 * These values come from
 * http://pdgprod.lbl.gov/~deg/AtomicNuclearProperties/HTML/deuterium_oxide_liquid.html for D2O and
 * http://pdg.lbl.gov/2018/AtomicNuclearProperties/HTML/water_liquid.html for H2O.
 */
static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6;
static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3;

/* Returns the average number of Cerenkov photons in the range 200-800 nm
 * produced by secondary particles in a muon shower.
 *
 * This comes from fitting the ratio # shower photons/rad loss to the function:
 *
 *     c0*(1-exp(-T/c1))
 *
 * I don't really have any good theoretical motivation for this. My initial
 * thought was that the number of photons should be roughly proportional to the
 * energy lost due to radiation which is why I chose to fit the ratio.
 *
 * This functional form just happened to fit the ratio as a function of energy
 * pretty well from 300 MeV to 10 GeV. At 10 GeV, it looks like the ratio is
 * starting to decrease so a different form for energies past that is probably
 * needed.
 *
 * `T0` is the initial kinetic energy of the electron in MeV and `rad` is the
 * energy lost due to radiation in MeV. */
double muon_get_shower_photons(double T0, double rad)
{
    return rad*(9.288929e+03*(1 - exp(-T0/8.403863e+02)));
}

void muon_get_position_distribution_parameters(double T0, double *a, double *b)
{
    /* Computes the gamma distribution parameters describing the longitudinal
     * profile of the photons generated from an electromagnetic shower.
     *
     * From the PDG "Passage of Particles" section 32.5:
     *
     * "The mean longitudinal profile of the energy deposition in an
     * electromagnetic cascade is reasonably well described by a gamma
     * distribution."
     *
     * Here we use a slightly different form of the gamma distribution:
     *
     *     f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a)
     *
     * I determined a and b by simulating high energy muons using
     * RAT-PAC and fitting the histogram of the position of all photons as a
     * function of the distance along the track length.
     *
     * Note: Unlike the case of a shower produced by an electron, the
     * distribution of photons from high energy muons does not seem to follow a
     * gamma distribution very well. In addition, the numbers I use here are
     * really approximate. The b parameter was obtained by a single degree
     * polynomial fit because it looked pretty good, but for the a parameter, I
     * couldn't find any functional form that would describe it well as a
     * function of energy and so I decided to just approximate it by a
     * constant.
     *
     * FIXME: Should update this in the future.
     *
     * `T` should be in units of MeV.
     *
     * Example:
     *
     *     double a, b;
     *     muon_get_position_distribution_parameters(1000.0, &a, &b);
     *     double pdf = gamma_pdf(x, a, b);
     *
     * See http://pdg.lbl.gov/2014/reviews/rpp2014-rev-passage-particles-matter.pdf.
     *
     * FIXME: Double check that this is correct for muons. */
    *b = fmax(20.0,-7.8 + 0.118928*T0);
    *a = 1.5;
}

double muon_get_angular_distribution_alpha(double T0)
{
    /* To describe the angular distribution of photons in an electromagnetic
     * shower I came up with the heuristic form:
     *
     *     f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * I simulated high energy muons using RAT-PAC in heavy water to fit
     * for the alpha and beta parameters as a function of energy and determined
     * that a reasonably good fit to the data was:
     *
     *     alpha = c0 + c1/log(c2*T0 + c3)
     *
     * where T0 is the initial energy of the muon in MeV and c0, c1, c2,
     * and c3 are constants which I fit out. */
    return fmin(100.0,8.238633e-01 + 3.896665e-03/log(1.581060e-05*T0 + 9.991576e-01));
}

double muon_get_angular_distribution_beta(double T0)
{
    /* To describe the angular distribution of photons in an electromagnetic
     * shower I came up with the heuristic form:
     *
     * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * I simulated high energy muons using RAT-PAC in heavy water to fit
     * for the alpha and beta parameters as a function of energy and determined
     * that a reasonably good fit to the data was:
     *
     * beta = c0 + c1/log(c2*T0 + c3)
     *
     * where T0 is the initial energy of the muon in MeV and c0, c1, c2,
     * and c3 are constants which I fit out. */
    return 2.236058e-01 + 5.376336e-03/log(1.200215e-05*T0 + 1.002540e+00);
}

double muon_get_delta_ray_photons(double T0)
{
    /* Returns the number of photons in the range 200-800 nm produced by delta
     * rays.
     *
     * This comes from a simple polynomial fit to the number of photons
     * generated by simulating muons using RAT-PAC in heavy water. */
    return fmax(0.0,-7532.39 + 39.4548*T0);
}

void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b)
{
    /* To describe the angular distribution of photons from delta rays I use
     * the same heuristic form as used for electromagnetic showers:
     *
     * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
     *
     * I simulated high energy muons using RAT-PAC in heavy water to fit
     * for the alpha and beta parameters as a function of energy and determined
     * that a reasonably good fit to the data was:
     *
     * beta = c0 + c1/log(c2*T0 + c3)
     *
     * where T0 is the initial energy of the muon in MeV and c0, c1, c2,
     * and c3 are constants which I fit out. */
    *a = fmin(100.0,3.463093e-01 + 1.110835e-02/log(5.662948e-06*T0 + 1.009215e+00));
    *b = 2.297358e-01 + 4.085721e-03/log(8.218774e-06*T0 + 1.007135e+00);

    return;
}

static int init()
{
    int i, j;
    char line[256];
    char *str;
    double value;
    int n;

    FILE *f = open_file("muE_deuterium_oxide_liquid.txt", "r");

    if (!f) {
        fprintf(stderr, "failed to open muE_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 10 lines since it's just a header. */
        if (i <= 10) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;
        else if (strstr(line, "Minimum ionization")) continue;
        else if (strstr(line, "Muon critical energy")) 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 10 lines since it's just a header. */
        if (i <= 10) continue;

        if (!len) continue;
        else if (line[0] == '#') continue;
        else if (strstr(line, "Minimum ionization")) continue;
        else if (strstr(line, "Muon critical energy")) continue;

        str = strtok(line," \n");

        j = 0;
        while (str) {
            value = strtod(str, NULL);
            switch (j) {
            case 0:
                x[n] = value;
                break;
            case 6:
                dEdx_rad[n] = value;
                break;
            case 7:
                dEdx[n] = value;
                break;
            case 8:
                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);

    initialized = 1;

    return 0;

err:
    fclose(f);

    return -1;
}

/* Returns the maximum kinetic energy for a muon in the range tables.
 *
 * If you call muon_get_range() or muon_get_dEdx() with a kinetic energy higher
 * you will get a GSL interpolation error. */
double muon_get_max_energy(void)
{
    if (!initialized) {
        if (init()) {
            exit(1);
        }
    }

    return x[size-1];
}

double muon_get_range(double T, double rho)
{
    /* Returns the approximate range a muon 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 muon_get_dEdx_rad(double T, double rho)
{
    /* Returns the approximate radiative dE/dx for a muon 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_rad->x[0]) return spline_dEdx_rad->y[0];

    return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
}

double muon_get_dEdx(double T, double rho)
{
    /* Returns the approximate dE/dx for a muon 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;
}