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, 33 insertions, 1 deletions

diff --git a/src/proton.c b/src/proton.c
index 922c552..f5a6b61 100644
--- a/src/proton.c
+++ b/src/proton.c
@@ -18,9 +18,12 @@
 
 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;
 
@@ -71,6 +74,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;
@@ -103,6 +107,9 @@ static int init()
             case 0:
                 x[n] = value;
                 break;
+            case 2:
+                dEdx_rad[n] = value;
+                break;
             case 3:
                 dEdx[n] = value;
                 break;
@@ -119,6 +126,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);
@@ -158,6 +169,27 @@ double proton_get_range(double T, double rho)
     return gsl_spline_eval(spline_range, T, acc_range)/rho;
 }
 
+double proton_get_dEdx_rad(double T, double rho)
+{
+    /* Returns the approximate radiative 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_rad->x[0]) return spline_dEdx_rad->y[0];
+
+    return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
+}
+
 double proton_get_dEdx(double T, double rho)
 {
     /* Returns the approximate dE/dx for a proton in water with kinetic energy