diff options
| -rw-r--r-- | Makefile | 6 | ||||
| -rw-r--r-- | muon.c | 36 | ||||
| -rw-r--r-- | muon.h | 1 | ||||
| -rw-r--r-- | test.c | 33 |
4 files changed, 73 insertions, 3 deletions
@@ -1,6 +1,8 @@ -CFLAGS=-O4 +CFLAGS=-O4 -Wall LDLIBS=-lm -lgsl -lgslcblas +all: test + calculate_limits: calculate_limits.c solid_angle.o: solid_angle.c @@ -8,6 +10,6 @@ solid_angle.o: solid_angle.c test: test.c solid_angle.c refractive_index.c muon.c clean: - rm calculate_limits + rm *.o calculate_limits test .PHONY: clean @@ -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 @@ -2,6 +2,7 @@ #define MUON_H double get_range(double T, double rho); +double get_E(double T, double x, double rho); double get_dEdx(double T, double rho); #endif @@ -109,6 +109,36 @@ int isclose(double a, double b, double rel_tol, double abs_tol) return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol); } +int test_muon_get_E(char *err) +{ + /* A very simple test to make sure the energy as a function of distance + * along the track makes sense. Should include more detailed tests later. */ + double T, E, range; + + /* Assume initial kinetic energy is 1 GeV. */ + T = 1000.0; + E = get_E(T,1e-9,1.0); + + /* At the beginning of the track we should have roughly the same energy. */ + if (!isclose(E, T, 1e-5, 0)) { + sprintf(err, "KE = %.5f, but expected %.5f", E, T); + return 1; + } + + /* Assume initial kinetic energy is 1 GeV. */ + T = 1000.0; + range = get_range(T,1.0); + E = get_E(T,range,1.0); + + /* At the end of the track we should have roughly the same energy. */ + if (!isclose(E, 0, 1e-5, 1e-5)) { + sprintf(err, "KE = %.5f, but expected %.5f", E, 0.0); + return 1; + } + + return 0; +} + int test_muon_get_range(char *err) { /* A very simple test to make sure we read in the PDG table correctly. */ @@ -198,7 +228,8 @@ struct tests { {test_solid_angle, "test_solid_angle"}, {test_refractive_index, "test_refractive_index"}, {test_muon_get_dEdx, "test_muon_get_dEdx"}, - {test_muon_get_range, "test_muon_get_range"} + {test_muon_get_range, "test_muon_get_range"}, + {test_muon_get_E, "test_muon_get_E"} }; int main(int argc, char **argv) |