add a function to compute the number of shower photons for muons

Similarly to electrons, I fit an analytic form to the ratio of the number of
photons produced via shower particles over the radiative energy loss. In this
case, I chose the functional form:

    ratio = a*(1-exp(-T/b))

since the ratio seemed to reach a constant value after a certain energy. I then
simulated a 10 GeV muon and it appears that the ratio might actually decrease
after that, so for higher energies I may have to come up with a different fit
function.

1 files changed, 23 insertions, 0 deletions

diff --git a/src/muon.c b/src/muon.c
index dabf9bc..aea2309 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -57,6 +57,29 @@ static gsl_spline *spline_range;
 static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6;
 static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3;
 
+/* Returns the average number of Cerenkov photons in the range 200-800 nm
+ * produced by secondary particles in a muon shower.
+ *
+ * This comes from fitting the ratio # shower photons/rad loss to the function:
+ *
+ *     c0*(1-exp(-T/c1))
+ *
+ * I don't really have any good theoretical motivation for this. My initial
+ * thought was that the number of photons should be roughly proportional to the
+ * energy lost due to radiation which is why I chose to fit the ratio.
+ *
+ * This functional form just happened to fit the ratio as a function of energy
+ * pretty well from 300 MeV to 10 GeV. At 10 GeV, it looks like the ratio is
+ * starting to decrease so a different form for energies past that is probably
+ * needed.
+ *
+ * `T0` is the initial kinetic energy of the electron in MeV and `rad` is the
+ * energy lost due to radiation in MeV. */
+double muon_get_shower_photons(double T0, double rad)
+{
+    return rad*(9.288929e+03*(1 - exp(-T0/8.403863e+02)));
+}
+
 void muon_get_position_distribution_parameters(double T0, double *a, double *b)
 {
     /* Computes the gamma distribution parameters describing the longitudinal