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, 219 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
index f8cedb9..7179f6b 100644
--- a/src/test.c
+++ b/src/test.c
@@ -19,6 +19,8 @@
 #include "likelihood.h"
 #include "id_particles.h"
 #include <gsl/gsl_spline2d.h>
+#include <gsl/gsl_errno.h> /* for gsl_strerror() */
+#include <gsl/gsl_randist.h>
 
 typedef int testFunction(char *err);
 
@@ -970,6 +972,218 @@ err:
     return 1;
 }
 
+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(cos_theta,alpha,beta,mu);
+}
+
+int test_electron_get_angular_pdf(char *err)
+{
+    /* Tests that the electron_get_angular_pdf() function integrates to 1. */
+    size_t i;
+    double params[3];
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
+    double result, error;
+    int status;
+    size_t nevals;
+    double T0;
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &gsl_electron_get_angular_pdf;
+    F.params = params;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        T0 = genrand_real2()*1000;
+
+        params[0] = electron_get_angular_distribution_alpha(T0);
+        params[1] = electron_get_angular_distribution_beta(T0);
+        params[2] = genrand_real2();
+
+        status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);
+
+        if (status) {
+            sprintf(err, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
+            goto err;
+        }
+
+        if (!isclose(result, 1.0, 0, 1e-3)) {
+            sprintf(err, "integral of electron_get_angular_pdf() returned %.5g, but expected %.5g", result, 1.0);
+            goto err;
+        }
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    return 0;
+
+err:
+    gsl_integration_cquad_workspace_free(w);
+
+    return 1;
+}
+
+int test_gamma_pdf(char *err)
+{
+    /* Tests the gamma_pdf() function. */
+    size_t i;
+    double x, k, theta, expected;
+    double result;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        x = genrand_real2();
+        k = genrand_real2();
+        theta = genrand_real2();
+
+        result = gamma_pdf(x,k,theta);
+
+        expected = gsl_ran_gamma_pdf(x,k,theta);
+
+        if (!isclose(result, expected, 0, 1e-9)) {
+            sprintf(err, "gamma_pdf() returned %.5g, but expected %.5g", result, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
+int test_solid_angle2(char *err)
+{
+    /* Tests the get_solid_angle()2 function. */
+    size_t i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+    double dir[3];
+    double L, r0, R;
+    double expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        rand_sphere(pmt);
+        rand_sphere(n);
+        r = genrand_real2();
+
+        dir[0] = pos[0] - pmt[0];
+        dir[1] = pos[1] - pmt[1];
+        dir[2] = pos[2] - pmt[2];
+
+        L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+        R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+        r0 = sqrt(R*R - L*L);
+
+        solid_angle = get_solid_angle2(L/r,r0/r);
+
+        expected = get_solid_angle(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, expected, 1e-4, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+int test_solid_angle_lookup(char *err)
+{
+    /* Tests the get_solid_angle_lookup() function. */
+    size_t i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+    double dir[3];
+    double L, r0, R;
+    double expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        rand_sphere(pmt);
+        MUL(pmt,800);
+        rand_sphere(n);
+        r = genrand_real2();
+
+        dir[0] = pos[0] - pmt[0];
+        dir[1] = pos[1] - pmt[1];
+        dir[2] = pos[2] - pmt[2];
+
+        L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+        R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+        r0 = sqrt(R*R - L*L);
+
+        solid_angle = get_solid_angle_lookup(L/r,r0/r);
+
+        expected = get_solid_angle(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, expected, 1e-2, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+int test_solid_angle_fast(char *err)
+{
+    /* Tests the get_solid_angle_lookup() function. */
+    size_t i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+    double dir[3];
+    double L, r0, R;
+    double expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        rand_sphere(pmt);
+        MUL(pmt,800);
+        rand_sphere(n);
+        r = genrand_real2();
+
+        dir[0] = pos[0] - pmt[0];
+        dir[1] = pos[1] - pmt[1];
+        dir[2] = pos[2] - pmt[2];
+
+        L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+        R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+        r0 = sqrt(R*R - L*L);
+
+        solid_angle = get_solid_angle_fast(pos,pmt,n,r);
+
+        expected = get_solid_angle(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, expected, 1e-2, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -999,6 +1213,11 @@ struct tests {
     {test_log_norm, "test_log_norm"},
     {test_trapz, "test_trapz"},
     {test_interp2d, "test_interp2d"},
+    {test_electron_get_angular_pdf, "test_electron_get_angular_pdf"},
+    {test_gamma_pdf, "test_gamma_pdf"},
+    {test_solid_angle2, "test_solid_angle2"},
+    {test_solid_angle_lookup, "test_solid_angle_lookup"},
+    {test_solid_angle_fast, "test_solid_angle_fast"},
 };
 
 int main(int argc, char **argv)