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

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;