add photons from delta rays to likelihood calculation

This commit updates the likelihood function to take into account Cerenkov light
produced from delta rays produced by muons. The angular distribution of this
light is currently assumed to be constant along the track and parameterized in
the same way as the Cerenkov light from an electromagnetic shower. Currently I
assume the light is produced uniformly along the track which isn't exactly
correct, but should be good enough.

1 files changed, 39 insertions, 0 deletions

diff --git a/src/electron.c b/src/electron.c
index 261ba39..c6b5ed8 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -161,6 +161,45 @@ static double electron_get_angular_pdf_norm(double alpha, double beta, double mu
     return result;
 }
 
+double electron_get_angular_pdf_delta_ray(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));
+}
+
 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