add code to expand the track of a particle using a KL expansion

To fit the path of muons and electrons I use the Karhunen-Loeve expansion of a
random 2D walk in the polar angle in x and y. This allows you to decompose the
path into a sum over sine functions whose coefficients become random variables.
The nice thing about fitting the path in this way is that you can capture
*most* of the variation in the path using a small number of variables by only
summing over the first N terms in the expansion and it is easy to calculate the
probability of the coefficients since they are all uncorrelated.

1 files changed, 82 insertions, 1 deletions

diff --git a/src/test.c b/src/test.c
index 747e664..b2faedc 100644
--- a/src/test.c
+++ b/src/test.c
@@ -8,6 +8,9 @@
 #include "misc.h"
 #include "sno_charge.h"
 #include <gsl/gsl_integration.h>
+#include "path.h"
+#include "random.h"
+#include "pdg.h"
 
 typedef int testFunction(char *err);
 
@@ -324,6 +327,83 @@ int test_logsumexp(char *err)
     return 0;
 }
 
+double getKineticEnergy(double x, double T0)
+{
+    return 1.0;
+}
+
+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 pos_expected[3], t_expected;
+    double pos[3], dir[3];
+    double E, mom, beta;
+    double s;
+
+    T0 = 1.0;
+    range = 1.0;
+    m = 1.0;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        pos0[0] = genrand_real2();
+        pos0[1] = genrand_real2();
+        pos0[2] = genrand_real2();
+
+        rand_sphere(dir0);
+
+        path *p = path_init(pos0,dir0,T0,range,0.1,getKineticEnergy,NULL,NULL,0,m);
+
+        for (j = 0; j < 100; j++) {
+            s = range*j/99;
+            pos_expected[0] = pos0[0] + dir0[0]*s;
+            pos_expected[1] = pos0[1] + dir0[1]*s;
+            pos_expected[2] = pos0[2] + dir0[2]*s;
+
+            /* Calculate total energy */
+            E = T0 + m;
+            mom = sqrt(E*E - m*m);
+            beta = mom/E;
+
+            t_expected = s/(beta*SPEED_OF_LIGHT);
+
+            path_eval(p,s,pos,dir,&T,&t);
+
+            for (k = 0; k < 3; k++) {
+                if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
+                    sprintf(err, "path_eval(%.2g) returned pos[%i] = %.5g, but expected %.5g", s, k, pos[k], pos_expected[k]);
+                    return 1;
+                }
+            }
+
+            for (k = 0; k < 3; k++) {
+                if (!isclose(dir[k], dir0[k], 1e-9, 1e-9)) {
+                    sprintf(err, "path_eval(%.2g) returned dir[%i] = %.5g, but expected %.5g", s, k, dir[k], dir0[k]);
+                    return 1;
+                }
+            }
+
+            if (!isclose(T, 1.0, 1e-9, 1e-9)) {
+                sprintf(err, "path_eval(%.2g) returned T = %.5g, but expected %.5g", s, T, 1.0);
+                return 1;
+            }
+
+            if (!isclose(t, t_expected, 1e-9, 1e-9)) {
+                sprintf(err, "path_eval(%.2g) returned t = %.5g, but expected %.5g", s, t, t_expected);
+                return 1;
+            }
+        }
+
+        path_free(p);
+    }
+
+    return 0;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -335,7 +415,8 @@ struct tests {
     {test_muon_get_range, "test_muon_get_range"},
     {test_muon_get_T, "test_muon_get_T"},
     {test_logsumexp, "test_logsumexp"},
-    {test_sno_charge_integral, "test_sno_charge_integral"}
+    {test_sno_charge_integral, "test_sno_charge_integral"},
+    {test_path, "test_path"}
 };
 
 int main(int argc, char **argv)