From b0c6382a80ef9704858c3a77511fa58f7949d7a3 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 22 May 2019 12:44:15 -0400 Subject: 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. --- src/likelihood.c | 3 +-- src/muon.c | 23 +++++++++++++++++++++++ src/muon.h | 1 + 3 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/likelihood.c b/src/likelihood.c index 245c20c..3b37bd5 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -188,8 +188,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; - /* FIXME: add shower photons for muons. */ - p->shower_photons = 0.0; + p->shower_photons = muon_get_shower_photons(T0, rad); /* Calculate the normalization constant for the angular distribution. */ norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o); 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 diff --git a/src/muon.h b/src/muon.h index e259999..a8d2ea1 100644 --- a/src/muon.h +++ b/src/muon.h @@ -30,6 +30,7 @@ * FIXME: Actually determine what this is. */ #define MUON_PHOTONS_PER_MEV 7368.0 +double muon_get_shower_photons(double T0, double rad); void muon_get_position_distribution_parameters(double T0, double *a, double *b); double muon_get_angular_distribution_alpha(double T0); double muon_get_angular_distribution_beta(double T0); -- cgit