update fit to fit for electrons and protons

1 files changed, 157 insertions, 17 deletions

diff --git a/src/test.c b/src/test.c
index 1a30976..2645982 100644
--- a/src/test.c
+++ b/src/test.c
@@ -14,6 +14,10 @@
 #include <gsl/gsl_spline.h>
 #include "vector.h"
 #include <gsl/gsl_statistics.h>
+#include "electron.h"
+#include "proton.h"
+#include "likelihood.h"
+#include "id_particles.h"
 
 typedef int testFunction(char *err);
 
@@ -111,17 +115,89 @@ struct solid_angle_results {
     {2.0,2.0,0.282707}
 };
 
+int test_proton_get_energy(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 T0, T;
+    particle *p;
+
+    /* Assume initial kinetic energy is 1 GeV. */
+    T0 = 1000.0;
+    p = particle_init(IDP_PROTON, T0,1.0,10000);
+    T = particle_get_energy(1e-9,p);
+
+    /* At the beginning of the track we should have roughly the same energy. */
+    if (!isclose(T, T0, 1e-5, 0)) {
+        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
+        goto err;
+    }
+
+    T = particle_get_energy(p->range,p);
+
+    /* 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);
+        goto err;
+    }
+
+    particle_free(p);
+
+    return 0;
+
+err:
+    particle_free(p);
+
+    return 1;
+}
+
+int test_electron_get_energy(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 T0, T;
+    particle *p;
+
+    /* Assume initial kinetic energy is 1 GeV. */
+    T0 = 1000.0;
+    p = particle_init(IDP_E_MINUS, T0,1.0,10000);
+    T = particle_get_energy(1e-9,p);
+
+    /* At the beginning of the track we should have roughly the same energy. */
+    if (!isclose(T, T0, 1e-5, 0)) {
+        sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
+        goto err;
+    }
+
+    T = particle_get_energy(p->range,p);
+
+    /* 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);
+        goto err;
+    }
+
+    particle_free(p);
+
+    return 0;
+
+err:
+    particle_free(p);
+
+    return 1;
+}
+
 int test_muon_get_energy(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 T0, T, range;
-    muon_energy *m;
+    double T0, T;
+    particle *p;
 
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
-    m = muon_init_energy(T0,1.0,10000);
-    T = muon_get_energy(1e-9,m);
+    p = particle_init(IDP_MU_MINUS, T0,1.0,10000);
+    T = particle_get_energy(1e-9,p);
 
     /* At the beginning of the track we should have roughly the same energy. */
     if (!isclose(T, T0, 1e-5, 0)) {
@@ -129,8 +205,7 @@ int test_muon_get_energy(char *err)
         goto err;
     }
 
-    range = get_range(T0,1.0);
-    T = muon_get_energy(range,m);
+    T = particle_get_energy(p->range,p);
 
     /* At the end of the track the energy should be approximately 0. */
     if (!isclose(T, 0, 1e-5, 1e-5)) {
@@ -138,15 +213,44 @@ int test_muon_get_energy(char *err)
         goto err;
     }
 
-    muon_free_energy(m);
+    particle_free(p);
 
     return 0;
 
 err:
-    muon_free_energy(m);
+    particle_free(p);
 
     return 1;
+}
 
+int test_proton_get_range(char *err)
+{
+    /* A very simple test to make sure we read in the PDG table correctly. */
+    double value;
+
+    value = proton_get_range(1.0,1.0);
+
+    if (!isclose(value, 2.458E-03, 1e-5, 0)) {
+        sprintf(err, "range = %.5f, but expected %.5f", value, 2.458E-03);
+        return 1;
+    }
+
+    return 0;
+}
+
+int test_electron_get_range(char *err)
+{
+    /* A very simple test to make sure we read in the PDG table correctly. */
+    double value;
+
+    value = electron_get_range(1.0,1.0);
+
+    if (!isclose(value, 4.367e-01, 1e-5, 0)) {
+        sprintf(err, "range = %.5f, but expected %.5f", value, 4.367e-01);
+        return 1;
+    }
+
+    return 0;
 }
 
 int test_muon_get_range(char *err)
@@ -154,7 +258,7 @@ int test_muon_get_range(char *err)
     /* A very simple test to make sure we read in the PDG table correctly. */
     double value;
 
-    value = get_range(1.0,1.0);
+    value = muon_get_range(1.0,1.0);
 
     if (!isclose(value, 2.071e-3, 1e-5, 0)) {
         sprintf(err, "range = %.5f, but expected %.5f", value, 1.863e-3);
@@ -164,12 +268,42 @@ int test_muon_get_range(char *err)
     return 0;
 }
 
+int test_proton_get_dEdx(char *err)
+{
+    /* A very simple test to make sure we read in the PDG table correctly. */
+    double value;
+
+    value = proton_get_dEdx(1.0,1.0);
+
+    if (!isclose(value, 2.608E+02, 1e-5, 0)) {
+        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 2.608E+02);
+        return 1;
+    }
+
+    return 0;
+}
+
+int test_electron_get_dEdx(char *err)
+{
+    /* A very simple test to make sure we read in the PDG table correctly. */
+    double value;
+
+    value = electron_get_dEdx(1.0,1.0);
+
+    if (!isclose(value, 1.862, 1e-5, 0)) {
+        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 1.862);
+        return 1;
+    }
+
+    return 0;
+}
+
 int test_muon_get_dEdx(char *err)
 {
     /* A very simple test to make sure we read in the PDG table correctly. */
     double value;
 
-    value = get_dEdx(1.0,1.0);
+    value = muon_get_dEdx(1.0,1.0);
 
     if (!isclose(value, 5.471, 1e-5, 0)) {
         sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 6.097);
@@ -341,10 +475,10 @@ int test_path(char *err)
     /* Tests the logsumexp() function. */
     int i, j, k;
     double pos0[3], dir0[3], T0, range, m;
-    double T, t;
+    double beta, t;
     double pos_expected[3], t_expected;
     double pos[3], dir[3];
-    double E, mom, beta;
+    double E, mom, beta0;
     double s;
     double theta0;
 
@@ -372,11 +506,11 @@ int test_path(char *err)
             /* Calculate total energy */
             E = T0 + m;
             mom = sqrt(E*E - m*m);
-            beta = mom/E;
+            beta0 = mom/E;
 
-            t_expected = s/(beta*SPEED_OF_LIGHT);
+            t_expected = s/(beta0*SPEED_OF_LIGHT);
 
-            path_eval(p,s,pos,dir,&T,&t,&theta0);
+            path_eval(p,s,pos,dir,&beta,&t,&theta0);
 
             for (k = 0; k < 3; k++) {
                 if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
@@ -392,8 +526,8 @@ int test_path(char *err)
                 }
             }
 
-            if (!isclose(T, 1.0, 1e-9, 1e-9)) {
-                sprintf(err, "path_eval(%.2g) returned T = %.5g, but expected %.5g", s, T, 1.0);
+            if (!isclose(beta, beta0, 1e-9, 1e-9)) {
+                sprintf(err, "path_eval(%.2g) returned beta = %.5g, but expected %.5g", s, beta, beta0);
                 return 1;
             }
 
@@ -737,6 +871,12 @@ struct tests {
     {test_get_path_length, "test_get_path_length"},
     {test_mean, "test_mean"},
     {test_std, "test_std"},
+    {test_electron_get_dEdx, "test_electron_get_dEdx"},
+    {test_electron_get_range, "test_electron_get_range"},
+    {test_electron_get_energy, "test_electron_get_energy"},
+    {test_proton_get_dEdx, "test_proton_get_dEdx"},
+    {test_proton_get_range, "test_proton_get_range"},
+    {test_proton_get_energy, "test_proton_get_energy"},
 };
 
 int main(int argc, char **argv)