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

diff --git a/src/misc.c b/src/misc.c
index 605f3ea..d5667e8 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -498,3 +498,11 @@ double std(const double *x, size_t n)
 
     return sqrt(sum/n);
 }
+
+double gamma_pdf(double x, double k, double theta)
+{
+    /* Returns the PDF for the gamma distribution.
+     *
+     * See https://en.wikipedia.org/wiki/Gamma_distribution. */
+    return pow(x,k-1)*exp(-x/theta)/(gsl_sf_gamma(k)*pow(theta,k));
+}