add energy dependent number of shower photons

This commit updates the code to calculate the number of Cerenkov photons from
secondary particles produced in an electromagnetic shower from electrons to use
an energy dependent formula I fit to data simulated with RAT-PAC.

3 files changed, 25 insertions, 8 deletions

diff --git a/src/electron.c b/src/electron.c
index 957b278..5b4133a 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -59,6 +59,29 @@ static gsl_spline *spline_range;
 static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33;
 static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33;
 
+/* Returns the average number of Cerenkov photons in the range 200-800 nm
+ * produced by secondary particles in an electron shower.
+ *
+ * This comes from fitting the ratio # shower photons/rad loss to the function:
+ *
+ *     c0 + c1/log(T0*c2 + c3)
+ *
+ * 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
+ * ratio turned out to vary from approximately 520 at low energies (10 MeV) to
+ * a roughly constant value of 420 at higher energies (> 1 GeV).
+ *
+ * This functional form just happened to fit the ratio as a function of energy
+ * pretty well.
+ *
+ * `T0` is the initial kinetic energy of the electron in MeV and `rad` is the
+ * energy lost due to radiation in MeV. */
+double electron_get_shower_photons(double T0, double rad)
+{
+    return rad*(406.745 + 0.271816/log(5.31309e-05*T0+1.00174));
+}
+
 void electron_get_position_distribution_parameters(double T0, double *a, double *b)
 {
     /* Computes the gamma distribution parameters describing the longitudinal
diff --git a/src/electron.h b/src/electron.h
index 98201b8..cc550a7 100644
--- a/src/electron.h
+++ b/src/electron.h
@@ -19,13 +19,7 @@
 
 #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
-
+double electron_get_shower_photons(double T0, double rad);
 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);
diff --git a/src/likelihood.c b/src/likelihood.c
index f0430cf..c53e1fc 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -107,7 +107,7 @@ 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;
+        p->shower_photons = electron_get_shower_photons(T0, rad);
 
         break;
     case IDP_MU_MINUS: