update likelihood function to fit electrons!

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(c2*T0 + c3)
    beta  = c0 + c1/log(c2*T0 + 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.

1 files changed, 194 insertions, 1 deletions

diff --git a/src/electron.c b/src/electron.c
index e8e24d3..261ba39 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -15,12 +15,17 @@
 #include "scattering.h"
 #include "pmt_response.h"
 #include "misc.h"
+#include <gsl/gsl_integration.h>
+#include <math.h> /* for fmax() */
 
 static int initialized = 0;
 
-static double *x, *dEdx, *csda_range;
+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;
 
@@ -36,6 +41,165 @@ static gsl_spline *spline_range;
 static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33;
 static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33;
 
+void electron_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 the b parameter by simulating high energy electrons using
+     * RAT-PAC and determined that it's roughly equal to the radiation length.
+     * To calculate the a parameter we use the formula from the PDG, i.e.
+     *
+     *     tmax = (a-1)/b = ln(E/E_C) - 0.5
+     *
+     * Therefore, we calculate a as:
+     *
+     *     a = tmax*b+1.
+     *
+     * `T` should be in units of MeV.
+     *
+     * Example:
+     *
+     *     double a, b;
+     *     electron_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. */
+    double tmax;
+
+    *b = RADIATION_LENGTH;
+    tmax = log(T0/ELECTRON_CRITICAL_ENERGY_D2O) - 0.5;
+    *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1);
+}
+
+double electron_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 electrons 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 electron in MeV and c0, c1, c2,
+     * and c3 are constants which I fit out. */
+    return 3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00);
+}
+
+double electron_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 electrons 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 electron in MeV and c0, c1, c2,
+     * and c3 are constants which I fit out. */
+    return 1.35372e-01 + 2.22344e-01/log(1.96249e-02*T0 + 1.23912e+00);
+}
+
+double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double beta, double mu)
+{
+    return exp(-pow(fabs(cos_theta-mu),alpha)/beta);
+}
+
+static double gsl_electron_get_angular_pdf(double cos_theta, void *params)
+{
+    double alpha = ((double *) params)[0];
+    double beta = ((double *) params)[1];
+    double mu = ((double *) params)[2];
+    return electron_get_angular_pdf_no_norm(cos_theta,alpha,beta,mu);
+}
+
+static double electron_get_angular_pdf_norm(double alpha, double beta, double mu)
+{
+    double params[3];
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
+    double result, error;
+    int status;
+    size_t nevals;
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &gsl_electron_get_angular_pdf;
+    F.params = params;
+
+    params[0] = alpha;
+    params[1] = beta;
+    params[2] = mu;
+
+    status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);
+
+    if (status) {
+        fprintf(stderr, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    return result;
+}
+
+double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu)
+{
+    /* Returns the probability density that a photon in the wavelength range
+     * 200 nm - 800 nm is emitted at an angle cos_theta.
+     *
+     * The angular distribution is modelled by the pdf:
+     *
+     *     f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
+     *
+     * where alpha and beta are constants which depend on the initial energy of
+     * the particle.
+     *
+     * There is no nice closed form solution for the integral of this function,
+     * so we just compute it on the fly. To make things faster, we keep track
+     * of the last alpha, beta, and mu parameters that were passed in so we
+     * avoid recomputing the normalization factor. */
+    size_t i;
+    static double last_alpha, last_beta, last_mu, norm;
+    static int first = 1;
+    static double x[10000], y[10000];
+
+    if (first || alpha != last_alpha || beta != last_beta || mu != last_mu) {
+        norm = electron_get_angular_pdf_norm(alpha, beta, mu);
+
+        last_alpha = alpha;
+        last_beta = beta;
+        last_mu = mu;
+
+        for (i = 0; i < LEN(x); i++) {
+            x[i] = -1.0 + 2.0*i/(LEN(x)-1);
+            y[i] = electron_get_angular_pdf_no_norm(x[i], alpha, beta, mu)/norm;
+        }
+
+        first = 0;
+    }
+
+    return interp1d(cos_theta,x,y,LEN(x));
+}
+
 static int init()
 {
     int i, j;
@@ -80,6 +244,7 @@ static int init()
     }
 
     x = malloc(sizeof(double)*n);
+    dEdx_rad = malloc(sizeof(double)*n);
     dEdx = malloc(sizeof(double)*n);
     csda_range = malloc(sizeof(double)*n);
     size = n;
@@ -112,6 +277,9 @@ static int init()
             case 0:
                 x[n] = value;
                 break;
+            case 2:
+                dEdx_rad[n] = value;
+                break;
             case 3:
                 dEdx[n] = value;
                 break;
@@ -128,6 +296,10 @@ static int init()
 
     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);
@@ -167,6 +339,27 @@ double electron_get_range(double T, double rho)
     return gsl_spline_eval(spline_range, T, acc_range)/rho;
 }
 
+double electron_get_dEdx_rad(double T, double rho)
+{
+    /* Returns the approximate radiative 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 (!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 electron_get_dEdx(double T, double rho)
 {
     /* Returns the approximate dE/dx for a electron in water with kinetic energy