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.

7 files changed, 281 insertions, 60 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
diff --git a/src/electron.h b/src/electron.h
index cb3f24e..20e5ea2 100644
--- a/src/electron.h
+++ b/src/electron.h
@@ -3,9 +3,17 @@
 
 #include <stddef.h> /* for size_t */
 
+/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy
+ * lost to radiation for electrons.
+ *
+ * FIXME: This is just a rough estimate, should use an energy dependent
+ * quantity from simulation. */
+#define ELECTRON_PHOTONS_PER_MEV 400.0
+
 void electron_get_position_distribution_parameters(double T0, double *a, double *b);
 double electron_get_angular_distribution_alpha(double T0);
 double electron_get_angular_distribution_beta(double T0);
+double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu);
 double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu);
 double electron_get_range(double T, double rho);
 double electron_get_dEdx_rad(double T, double rho);
diff --git a/src/likelihood.c b/src/likelihood.c
index 786b221..8e12f57 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -42,16 +42,19 @@ particle *particle_init(int id, double T0, size_t n)
      * which we solve by dividing the track up into `n` segments and then
      * numerically computing the energy at each step. */
     size_t i;
-    double dEdx;
+    double dEdx, rad;
 
     particle *p = malloc(sizeof(particle));
     p->id = id;
     p->x = malloc(sizeof(double)*n);
     p->T = malloc(sizeof(double)*n);
     p->n = n;
-    p->rad = 0.0;
     p->a = 0.0;
     p->b = 0.0;
+    p->shower_photons = 0.0;
+    p->delta_ray_a = 0.0;
+    p->delta_ray_b = 0.0;
+    p->delta_ray_photons = 0.0;
 
     p->x[0] = 0;
     p->T[0] = T0;
@@ -67,14 +70,17 @@ particle *particle_init(int id, double T0, size_t n)
         p->a = electron_get_angular_distribution_alpha(T0);
         p->b = electron_get_angular_distribution_beta(T0);
 
+        p->delta_ray_photons = 0.0;
+
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
+        rad = 0.0;
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
             dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
             dEdx = electron_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            rad += dEdx*(p->x[i]-p->x[i-1]);
         }
         /* Make sure that the energy is zero at the last step. This is so that
          * when we try to bisect the point along the track where the speed of
@@ -85,6 +91,8 @@ particle *particle_init(int id, double T0, size_t n)
          * range here using the dE/dx table instead of reading in the range. */
         p->T[n-1] = 0;
 
+        p->shower_photons = rad*ELECTRON_PHOTONS_PER_MEV;
+
         break;
     case IDP_MU_MINUS:
     case IDP_MU_PLUS:
@@ -96,14 +104,21 @@ particle *particle_init(int id, double T0, size_t n)
          * load the correct table based on the media. */
         p->range = muon_get_range(T0, HEAVY_WATER_DENSITY);
 
+        p->a = muon_get_angular_distribution_alpha(T0);
+        p->b = muon_get_angular_distribution_beta(T0);
+
+        muon_get_delta_ray_distribution_parameters(T0, &p->delta_ray_a, &p->delta_ray_b);
+        p->delta_ray_photons = muon_get_delta_ray_photons(T0);
+
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
+        rad = 0.0;
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
             dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
             dEdx = muon_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            rad += dEdx*(p->x[i]-p->x[i-1]);
         }
         /* Make sure that the energy is zero at the last step. This is so that
          * when we try to bisect the point along the track where the speed of
@@ -114,6 +129,8 @@ particle *particle_init(int id, double T0, size_t n)
          * range here using the dE/dx table instead of reading in the range. */
         p->T[n-1] = 0;
 
+        p->shower_photons = rad*MUON_PHOTONS_PER_MEV;
+
         break;
     case IDP_PROTON:
         p->mass = PROTON_MASS;
@@ -121,14 +138,17 @@ particle *particle_init(int id, double T0, size_t n)
          * tables for heavy water for protons. */
         p->range = proton_get_range(T0, WATER_DENSITY);
 
+        p->delta_ray_photons = 0.0;
+
         /* Now we loop over the points along the track and assume dE/dx is
          * constant between points. */
+        rad = 0.0;
         for (i = 1; i < n; i++) {
             p->x[i] = p->range*i/(n-1);
             dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY);
             p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
             dEdx = proton_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
-            p->rad += dEdx*(p->x[i]-p->x[i-1]);
+            rad += dEdx*(p->x[i]-p->x[i-1]);
         }
         /* Make sure that the energy is zero at the last step. This is so that
          * when we try to bisect the point along the track where the speed of
@@ -139,6 +159,8 @@ particle *particle_init(int id, double T0, size_t n)
          * range here using the dE/dx table instead of reading in the range. */
         p->T[n-1] = 0;
 
+        p->shower_photons = rad*PROTON_PHOTONS_PER_MEV;
+
         break;
     default:
         fprintf(stderr, "unknown particle id %i\n", id);
@@ -170,9 +192,9 @@ void particle_free(particle *p)
     free(p);
 }
 
-static double get_expected_charge_shower(double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double alpha, double beta, double rad)
+static double get_expected_charge_shower(particle *p, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double *q_delta_ray, double *q_indirect_delta_ray)
 {
-    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter;
+    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant;
 
     SUB(pmt_dir,pmt_pos,pos);
 
@@ -180,6 +202,9 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p
 
     cos_theta_pmt = DOT(pmt_dir,pmt_normal);
 
+    charge = 0.0;
+    *q_delta_ray = 0.0;
+    *q_indirect_delta_ray = 0.0;
     *reflected = 0.0;
     if (cos_theta_pmt > 0) return 0.0;
 
@@ -201,11 +226,26 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p
 
     l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
 
-    charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI);
+    constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*omega/(2*M_PI);
+
+    /* Note: We assume here that the peak of the angular distribution is at the
+     * Cerenkov angle for a particle with beta = 1. This seems to be an OK
+     * assumption for high energy showers, but is not exactly correct for
+     * showers from lower energy electrons or for delta rays from lower energy
+     * muons. It seems good enough, but in the future it would be nice to
+     * parameterize this. */
+    if (p->shower_photons > 0)
+        charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o);
+
+    if (p->delta_ray_photons > 0)
+        *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o);
 
     scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o);
 
     *reflected = (f_reflec + 1.0 - scatter)*charge;
+    *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray);
+
+    *q_delta_ray *= f*scatter;
 
     return f*charge*scatter;
 }
@@ -345,14 +385,14 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
         return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
 }
 
-static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad)
+static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
 {
     size_t i;
     static double qs[MAX_NPOINTS];
     static double q_indirects[MAX_NPOINTS];
     static double ts[MAX_NPOINTS];
     static double ts2[MAX_NPOINTS];
-    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o;
+    double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray;
 
     if (n > MAX_NPOINTS) {
         fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
@@ -364,6 +404,8 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
     n_d2o = get_index_snoman_d2o(wavelength0);
     n_h2o = get_index_snoman_h2o(wavelength0);
 
+    double range = x[n-1] - x[0];
+
     for (i = 0; i < n; i++) {
         pos[0] = pos0[0] + x[i]*dir0[0];
         pos[1] = pos0[1] + x[i]*dir0[1];
@@ -373,8 +415,9 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
 
         t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
 
-        qs[i] = get_expected_charge_shower(pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, p->a, p->b, rad)*pdf[i];
-        q_indirects[i] *= pdf[i];
+        qs[i] = get_expected_charge_shower(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray);
+        qs[i] = qs[i]*pdf[i] + q_delta_ray/range;
+        q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/range;
         ts[i] = t*qs[i];
         ts2[i] = t*t*qs[i];
     }
@@ -593,10 +636,11 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     charge = abs_prob*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI);
 
-    /* Add expected number of photons from electromagnetic shower for electrons
-     * and positrons. */
-    if (p->id == IDP_E_MINUS || p->id == IDP_E_PLUS)
-        charge += abs_prob*get_weighted_quantum_efficiency()*p->rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI);
+    /* Add expected number of photons from electromagnetic shower. */
+    if (p->shower_photons > 0)
+        charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI);
+    if (p->delta_ray_photons > 0)
+        charge += abs_prob*get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI);
 
     *mu_reflected = f_reflected*charge;
 
@@ -891,6 +935,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
         npoints = (size_t) (range/dx + 0.5);
 
         if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS;
+        if (npoints > MAX_NPOINTS) npoints = MAX_NPOINTS;
 
         /* Calculate total energy */
         E0 = v[j].T0 + p->mass;
@@ -903,32 +948,43 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
         if (!fast) {
             path = path_init(v[j].pos, v[j].dir, v[j].T0, range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass);
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
+            switch (v[j].id) {
+            case IDP_E_MINUS:
+            case IDP_E_PLUS:
                 electron_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b);
+                break;
+            case IDP_MU_MINUS:
+            case IDP_MU_PLUS:
+                muon_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b);
+                break;
+            default:
+                fprintf(stderr, "unknown particle id %i\n", v[j].id);
+                exit(1);
+            }
 
-                /* Lower and upper limits to integrate for the electromagnetic
-                 * shower. */
-                xlo = 0.0;
-                xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b);
+            /* Lower and upper limits to integrate for the electromagnetic
+             * shower. */
+            xlo = 0.0;
+            xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b);
 
-                /* Number of points to sample along the longitudinal direction for
-                 * the electromagnetic shower. */
-                npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5);
+            /* Number of points to sample along the longitudinal direction for
+             * the electromagnetic shower. */
+            npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5);
 
-                if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS;
+            if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS;
+            if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS;
 
-                x = malloc(sizeof(double)*npoints_shower);
-                pdf = malloc(sizeof(double)*npoints_shower);
-                for (i = 0; i < npoints_shower; i++) {
-                    x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1);
-                    pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
-                }
+            x = malloc(sizeof(double)*npoints_shower);
+            pdf = malloc(sizeof(double)*npoints_shower);
+            for (i = 0; i < npoints_shower; i++) {
+                x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1);
+                pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
+            }
 
-                /* Make sure the PDF is normalized. */
-                double norm = trapz(pdf,x[1]-x[0],npoints_shower);
-                for (i = 0; i < npoints_shower; i++) {
-                    pdf[i] /= norm;
-                }
+            /* Make sure the PDF is normalized. */
+            double norm = trapz(pdf,x[1]-x[0],npoints_shower);
+            for (i = 0; i < npoints_shower; i++) {
+                pdf[i] /= norm;
             }
         }
 
@@ -982,27 +1038,16 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
                 ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov);
             }
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
-                /* If we are fitting for an electron or positron, we need to
-                 * estimate the expected number of photons produced by the
-                 * electromagnetic shower.
-                 *
-                 * Technically we should also have a term for the electromagnetic
-                 * shower produced by muons, but for the energy range I'm
-                 * considering right now, it is very small. */
-                integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j],p->rad);
-                mu_indirect[j] += q_indirect;
-                ts_shower[i][j] = v[j].t0 + result;
-            }
+            integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]);
+            mu_indirect[j] += q_indirect;
+            ts_shower[i][j] = v[j].t0 + result;
         }
 
         if (!fast) {
             path_free(path);
 
-            if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) {
-                free(x);
-                free(pdf);
-            }
+            free(x);
+            free(pdf);
         }
 
         particle_free(p);
diff --git a/src/likelihood.h b/src/likelihood.h
index 84fa7ff..4b72742 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -45,13 +45,6 @@
 #define GTVALID 400.0
 #define BETA_MIN 0.8
 
-/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy
- * lost to radiation for electrons.
- *
- * FIXME: This is just a rough estimate, should use an energy dependent
- * quantity from simulation. */
-#define PHOTONS_PER_MEV 400.0
-
 /* Maximum number of vertices to fit. */
 #define MAX_VERTICES 10
 
@@ -70,12 +63,15 @@ typedef struct particle {
     int id;
     double mass;
     double range;
-    double rad;
     double *x;
     double *T;
     size_t n;
     double a;
     double b;
+    double shower_photons;
+    double delta_ray_a;
+    double delta_ray_b;
+    double delta_ray_photons;
 } particle;
 
 particle *particle_init(int id, double T0, size_t n);
diff --git a/src/muon.c b/src/muon.c
index 76e5223..6ab0aa7 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -40,6 +40,116 @@ static gsl_spline *spline_range;
 static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6;
 static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3;
 
+void muon_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 muons 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;
+     *     muon_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.
+     *
+     * FIXME: Double check that this is correct for muons. */
+    double tmax;
+
+    *b = RADIATION_LENGTH;
+    tmax = log(T0/MUON_CRITICAL_ENERGY_D2O) - 0.5;
+    *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1);
+}
+
+double muon_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 muons 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 muon in MeV and c0, c1, c2,
+     * and c3 are constants which I fit out. */
+    return 8.238633e-01 + 3.896665e-03/log(1.581060e-05*T0 + 9.991576e-01);
+}
+
+double muon_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 muons 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 muon in MeV and c0, c1, c2,
+     * and c3 are constants which I fit out. */
+    return 2.236058e-01 + 5.376336e-03/log(1.200215e-05*T0 + 1.002540e+00);
+}
+
+double muon_get_delta_ray_photons(double T0)
+{
+    /* Returns the number of photons in the range 200-800 nm produced by delta
+     * rays.
+     *
+     * This comes from a simple polynomial fit to the number of photons
+     * generated by simulating muons using RAT-PAC in heavy water. */
+    return fmax(0.0,-7532.39 + 39.4548*T0);
+}
+
+void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b)
+{
+    /* To describe the angular distribution of photons from delta rays I use
+     * the same heuristic form as used for electromagnetic showers:
+     *
+     * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
+     *
+     * I simulated high energy muons 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 muon in MeV and c0, c1, c2,
+     * and c3 are constants which I fit out. */
+    *a = 3.463093e-01 + 1.110835e-02/log(5.662948e-06*T0 + 1.009215e+00);
+    *b = 2.297358e-01 + 4.085721e-03/log(8.218774e-06*T0 + 1.007135e+00);
+
+    return;
+}
+
 static int init()
 {
     int i, j;
diff --git a/src/muon.h b/src/muon.h
index 586d67d..456c11e 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -5,6 +5,20 @@
 
 #define EULER_CONSTANT 0.57721
 
+/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy
+ * lost to radiation for electrons.
+ *
+ * FIXME: This is just a rough estimate, should use an energy dependent
+ * quantity from simulation.
+ *
+ * FIXME: Actually determine what this is. */
+#define MUON_PHOTONS_PER_MEV 7368.0
+
+void muon_get_position_distribution_parameters(double T0, double *a, double *b);
+double muon_get_angular_distribution_alpha(double T0);
+double muon_get_angular_distribution_beta(double T0);
+void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b);
+double muon_get_delta_ray_photons(double T0);
 double muon_get_range(double T, double rho);
 double muon_get_dEdx_rad(double T, double rho);
 double muon_get_dEdx(double T, double rho);
diff --git a/src/proton.h b/src/proton.h
index 008f835..cffeb9f 100644
--- a/src/proton.h
+++ b/src/proton.h
@@ -3,6 +3,15 @@
 
 #include <stddef.h> /* for size_t */
 
+/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy
+ * lost to radiation for electrons.
+ *
+ * FIXME: This is just a rough estimate, should use an energy dependent
+ * quantity from simulation.
+ *
+ * FIXME: Actually determine what this is. */
+#define PROTON_PHOTONS_PER_MEV 400.0
+
 double proton_get_range(double T, double rho);
 double proton_get_dEdx_rad(double T, double rho);
 double proton_get_dEdx(double T, double rho);