diff options
| -rw-r--r-- | src/electron.c | 23 | ||||
| -rw-r--r-- | src/electron.h | 8 | ||||
| -rw-r--r-- | src/likelihood.c | 2 |
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: |