diff options
Diffstat (limited to 'muon.c')
| -rw-r--r-- | muon.c | 36 |
1 files changed, 36 insertions, 0 deletions
@@ -4,6 +4,7 @@ #include <stdlib.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_spline.h> +#include <math.h> static int initialized = 0; @@ -16,6 +17,8 @@ static gsl_spline *spline_dEdx; static gsl_interp_accel *acc_range; static gsl_spline *spline_range; +static const double MUON_CRITICAL_ENERGY = 1.029e6; + static int init() { int i, j; @@ -153,6 +156,39 @@ double get_range(double T, double rho) return gsl_spline_eval(spline_range, T, acc_range)/rho; } +double get_E(double T, double x, double rho) +{ + /* Returns the approximate energy of a muon in water after travelling `x` + * cm with an initial kinetic energy `T`. + * + * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3. + * + * Return value is in MeV. + * + * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */ + double a, b, range, E; + + if (!initialized) { + if (init()) { + exit(1); + } + } + + range = gsl_spline_eval(spline_range, T, acc_range)/rho; + /* FIXME: Need to double check if it's kosher to be using kinetic energies + * here instead of the total energy. Equation 1 in the document uses the + * total energy, but here I'm using the critical energy in kinetic energy, + * so I should check to see if I need to convert both. */ + b = log(1 + T/MUON_CRITICAL_ENERGY)/range; + a = MUON_CRITICAL_ENERGY*b; + + E = T + a*(1-exp(b*x))/b; + + if (E < 0) return 0; + + return E; +} + double get_dEdx(double T, double rho) { /* Returns the approximate dE/dx for a muon in water with kinetic energy |