add a function to compute the energy along a muon track

4 files changed, 73 insertions, 3 deletions

diff --git a/Makefile b/Makefile
index 0e68776..3e0c0ee 100644
--- a/Makefile
+++ b/Makefile
@@ -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
diff --git a/muon.c b/muon.c
index 98a74ea..12a68b2 100644
--- a/muon.c
+++ b/muon.c
@@ -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
diff --git a/muon.h b/muon.h
index 461339b..d8efcb9 100644
--- a/muon.h
+++ b/muon.h
@@ -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
diff --git a/test.c b/test.c
index 1491711..3448565 100644
--- a/test.c
+++ b/test.c
@@ -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)