update muon kinetic energy calculation

This commit updates the calculation of the muon kinetic energy as a function of
distance along the track. Previously I was using an approximation from the PDG,
but it doesn't seem to be very accurate and won't generalize to the case of
electrons. The kinetic energy is now calculated using the tabulated values of
dE/dx as a function of energy.

1 files changed, 14 insertions, 14 deletions

diff --git a/src/test.c b/src/test.c
index 6d7d934..b74728d 100644
--- a/src/test.c
+++ b/src/test.c
@@ -113,26 +113,26 @@ int test_muon_get_T(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;
+    double T0, T, range;
+    muon_energy *m;
 
     /* Assume initial kinetic energy is 1 GeV. */
-    T = 1000.0;
-    E = get_T(T,1e-9,1.0);
+    T0 = 1000.0;
+    m = muon_init_energy(T0,1.0,10000);
+    T = muon_get_energy(1e-9,m);
 
     /* 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);
+    if (!isclose(T, T0, 1e-5, 0)) {
+        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
         return 1;
     }
 
-    /* Assume initial kinetic energy is 1 GeV. */
-    T = 1000.0;
-    range = get_range(T,1.0);
-    E = get_T(T,range,1.0);
+    range = get_range(T0,1.0);
+    T = muon_get_energy(range,m);
 
-    /* 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);
+    /* At the end of the track the energy should be approximately 0. */
+    if (!isclose(T, 0, 1e-5, 1e-5)) {
+        sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
         return 1;
     }
 
@@ -319,7 +319,7 @@ int test_logsumexp(char *err)
     return 0;
 }
 
-double getKineticEnergy(double x, double T0)
+double getKineticEnergy(double x, void *params)
 {
     return 1.0;
 }
@@ -349,7 +349,7 @@ int test_path(char *err)
 
         rand_sphere(dir0);
 
-        path *p = path_init(pos0,dir0,T0,range,0.1,getKineticEnergy,NULL,NULL,0,m);
+        path *p = path_init(pos0,dir0,T0,range,0.1,getKineticEnergy,NULL,NULL,NULL,0,m);
 
         for (j = 0; j < 100; j++) {
             s = range*j/99;