update shower position distribution parameters for muons

This commit updates the a and b parameters for the gamma distribution used to
describe the position distribution of shower photons produced along the
direction of the muon. Previously I had been assuming b was equal to the
radiation length and using a formula from the PDG to calculate a from that.
However, this formula doesn't seem to be valid for muons (the formula comes
from a section describing the shower profile of electrons and gammas, so it's
not surprising). Therefore, now we don't assume any relationship between a and
b.

Now, the value of a is approximated by a constant since I couldn't find any
functional relationship as a function of energy to describe a very well (and
it's approximately constant), and b is approximated by a single degree
polynomial fit to the values I got from simulating muons in RAT-PAC as a
function of energy.

Note that looking at the simulation data it seems like the position
distribution of shower photons from muons isn't even very well described by a
gamma distribution, so in the future it might be a good idea to come up with a
better parameterization.

Even if I stick with the gamma distribution, it would be good to revisit this
in the future and fit for a and b over a wider range of energies.

1 files changed, 14 insertions, 12 deletions

diff --git a/src/muon.c b/src/muon.c
index 9d6c91a..14574d3 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -95,15 +95,20 @@ void muon_get_position_distribution_parameters(double T0, double *a, double *b)
      *
      *     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.
+     * I determined a and b by simulating high energy muons using
+     * RAT-PAC and fitting the histogram of the position of all photons as a
+     * function of the distance along the track length.
      *
-     *     tmax = (a-1)/b = ln(E/E_C) - 0.5
+     * Note: Unlike the case of a shower produced by an electron, the
+     * distribution of photons from high energy muons does not seem to follow a
+     * gamma distribution very well. In addition, the numbers I use here are
+     * really approximate. The b parameter was obtained by a single degree
+     * polynomial fit because it looked pretty good, but for the a parameter, I
+     * couldn't find any functional form that would describe it well as a
+     * function of energy and so I decided to just approximate it by a
+     * constant.
      *
-     * Therefore, we calculate a as:
-     *
-     *     a = tmax*b+1.
+     * FIXME: Should update this in the future.
      *
      * `T` should be in units of MeV.
      *
@@ -116,11 +121,8 @@ void muon_get_position_distribution_parameters(double T0, double *a, double *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);
+    *b = -7.8 + 0.118928*T0;
+    *a = 1.5;
 }
 
 double muon_get_angular_distribution_alpha(double T0)