move everything to src directory

1 files changed, 360 insertions, 0 deletions

diff --git a/src/test.c b/src/test.c
new file mode 100644
index 0000000..747e664
--- /dev/null
+++ b/src/test.c
@@ -0,0 +1,360 @@
+#include "solid_angle.h"
+#include <math.h>
+#include <stdio.h>
+#include "optics.h"
+#include "muon.h"
+#include "mt19937ar.h"
+#include <gsl/gsl_sf_gamma.h>
+#include "misc.h"
+#include "sno_charge.h"
+#include <gsl/gsl_integration.h>
+
+typedef int testFunction(char *err);
+
+/* Table of some of the tabulated values of the refractive index of water as a
+ * function of wavelength and temperature. In all cases, I just used the values
+ * for standard atmospheric pressure and assume this corresponds approximately
+ * to a density of 1000 kg/m^3.
+ *
+ * See Table 7 in https://aip.scitation.org/doi/pdf/10.1063/1.555859. */
+
+struct refractive_index_results {
+    double p;
+    double T;
+    double wavelength;
+    double n;
+} refractive_index_results[] = {
+    {1.0, 0 ,  226.50, 1.39468},
+    {1.0, 10,  226.50, 1.39439},
+    {1.0, 20,  226.50, 1.39353},
+    {1.0, 30,  226.50, 1.39224},
+    {1.0, 0 ,  404.41, 1.34431},
+    {1.0, 10,  404.41, 1.34404},
+    {1.0, 20,  404.41, 1.34329},
+    {1.0, 30,  404.41, 1.34218},
+    {1.0, 0 ,  589.00, 1.33447},
+    {1.0, 10,  589.00, 1.33422},
+    {1.0, 20,  589.00, 1.33350},
+    {1.0, 30,  589.00, 1.33243},
+    {1.0, 0 ,  632.80, 1.33321},
+    {1.0, 10,  632.80, 1.33296},
+    {1.0, 20,  632.80, 1.33224},
+    {1.0, 30,  632.80, 1.33118},
+    {1.0, 0 , 1013.98, 1.32626},
+    {1.0, 10, 1013.98, 1.32604},
+    {1.0, 20, 1013.98, 1.32537},
+    {1.0, 30, 1013.98, 1.32437},
+    {1.0, 0 , 2325.42, 1.27663},
+    {1.0, 10, 2325.42, 1.27663},
+    {1.0, 20, 2325.42, 1.27627},
+    {1.0, 30, 2325.42, 1.27563},
+};
+
+/* Table of the values of solid angle for various values of r0/r and L/r.
+ *
+ * See Table 1 in http://www.umich.edu/~ners312/CourseLibrary/SolidAngleOfADiskOffAxis.pdf. */
+
+struct solid_angle_results {
+    double L;
+    double r0;
+    double omega;
+} solid_angle_results[] = {
+    {0.5,0.0,3.4732594},
+    {0.5,0.2,3.4184435},
+    {0.5,0.4,3.2435434},
+    {0.5,0.6,2.9185178},
+    {0.5,0.8,2.4122535},
+    {0.5,1.0,1.7687239},
+    {0.5,1.2,1.1661307},
+    {0.5,1.4,0.7428889},
+    {0.5,1.6,0.4841273},
+    {0.5,1.8,0.3287007},
+    {0.5,2.0,0.2324189},
+    {1.0,0.0,1.8403024},
+    {1.0,0.2,1.8070933},
+    {1.0,0.4,1.7089486},
+    {1.0,0.6,1.5517370},
+    {1.0,0.8,1.3488367},
+    {1.0,1.0,1.1226876},
+    {1.0,1.2,0.9003572},
+    {1.0,1.4,0.7039130},
+    {1.0,1.6,0.5436956},
+    {1.0,1.8,0.4195415},
+    {1.0,2.0,0.3257993},
+    {1.5,0.0,1.0552591},
+    {1.5,0.2,1.0405177},
+    {1.5,0.4,0.9975504},
+    {1.5,0.6,0.9301028},
+    {1.5,0.8,0.8441578},
+    {1.5,1.0,0.7472299},
+    {1.5,1.2,0.6472056},
+    {1.5,1.4,0.5509617},
+    {1.5,1.6,0.4632819},
+    {1.5,1.8,0.3866757},
+    {1.5,2.0,0.3217142},
+    {2.0,0.0,0.6633335},
+    {2.0,0.2,0.6566352},
+    {2.0,0.4,0.6370508},
+    {2.0,0.6,0.6060694},
+    {2.0,0.8,0.5659755},
+    {2.0,1.0,0.5195359},
+    {2.0,1.2,0.4696858},
+    {2.0,1.4,0.4191714},
+    {2.0,1.6,0.3702014},
+    {2.0,1.8,0.3243908},
+    {2.0,2.0,0.282707}
+};
+
+int isclose(double a, double b, double rel_tol, double abs_tol)
+{
+    /* Returns 1 if a and b are "close". This algorithm is taken from Python's
+     * math.isclose() function.
+     *
+     * See https://www.python.org/dev/peps/pep-0485/. */
+    return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
+}
+
+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;
+
+    /* Assume initial kinetic energy is 1 GeV. */
+    T = 1000.0;
+    E = get_T(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_T(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. */
+    double value;
+
+    value = get_range(1.0,1.0);
+
+    if (!isclose(value, 1.863e-3, 1e-5, 0)) {
+        sprintf(err, "range = %.5f, but expected %.5f", value, 1.863e-3);
+        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);
+
+    if (!isclose(value, 6.097, 1e-5, 0)) {
+        sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 6.097);
+        return 1;
+    }
+
+    return 0;
+}
+
+int test_refractive_index(char *err)
+{
+    /* Tests the get_index() function. */
+    int i;
+    double n;
+    struct refractive_index_results result;
+
+    for (i = 0; i < sizeof(refractive_index_results)/sizeof(struct refractive_index_results); i++) {
+        result = refractive_index_results[i];
+        n = get_index(result.p, result.wavelength, result.T);
+
+        if (!isclose(n, result.n, 1e-2, 0)) {
+            sprintf(err, "n = %.5f, but expected %.5f", n, result.n);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+int test_solid_angle_approx(char *err)
+{
+    /* Tests the get_solid_angle_approx() function. */
+    int i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+
+    solid_angle = get_solid_angle_approx(pos,pmt,n,r);
+
+    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-2, 0)) {
+        sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
+        return 1;
+    }
+
+    for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
+        pos[0] = solid_angle_results[i].r0*r;
+        pos[2] = solid_angle_results[i].L*r;
+
+        solid_angle = get_solid_angle_approx(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-2, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+int test_solid_angle(char *err)
+{
+    /* Tests the get_solid_angle() function. */
+    int i;
+    double pmt[3] = {0,0,0};
+    double pos[3] = {0,0,1};
+    double n[3] = {0,0,1};
+    double r = 1.0;
+    double solid_angle;
+
+    solid_angle = get_solid_angle(pos,pmt,n,r);
+
+    if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-9, 0)) {
+        sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
+        return 1;
+    }
+
+    for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
+        pos[0] = solid_angle_results[i].r0*r;
+        pos[2] = solid_angle_results[i].L*r;
+
+        solid_angle = get_solid_angle(pos,pmt,n,r);
+
+        if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
+            sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+static double sno_charge(double q, void *params)
+{
+    int n = ((int *) params)[0];
+    return pq(q,n);
+}
+
+int test_sno_charge_integral(char *err)
+{
+    /* Test that the single PE charge distribution integrates to one. */
+    double result, error;
+    gsl_function F;
+    size_t nevals;
+    gsl_integration_cquad_workspace *w;
+    int params[1];
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &sno_charge;
+    params[0] = 1;
+    F.params = params;
+
+    init_charge();
+
+    gsl_integration_cquad(&F, -10.0, 1000.0, 0, 1e-9, w, &result, &error, &nevals);
+
+    if (!isclose(result, 1.0, 1e-9, 1e-9)) {
+        sprintf(err, "integral of single PE charge distribution is %.2f", result);
+        goto err;
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    return 0;
+
+err:
+    gsl_integration_cquad_workspace_free(w);
+    return 1;
+}
+
+int test_logsumexp(char *err)
+{
+    /* Tests the logsumexp() function. */
+    int i, j;
+    double logp[100], mu, result, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        mu = genrand_real2();
+
+        for (j = 0; j < sizeof(logp)/sizeof(double); j++) {
+            logp[j] = -mu + j*log(mu) - gsl_sf_lnfact(j);
+        }
+        result = logsumexp(logp, sizeof(logp)/sizeof(double));
+
+        expected = 0.0;
+
+        if (!isclose(result, expected, 1e-9, 1e-9)) {
+            sprintf(err, "logsumexp(logp) = %.5g, but expected %.5g", result, expected);
+            return 1;
+        }
+    }
+
+    return 0;
+}
+
+struct tests {
+    testFunction *test;
+    char *name;
+} tests[] = {
+    {test_solid_angle, "test_solid_angle"},
+    {test_solid_angle_approx, "test_solid_angle_approx"},
+    {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_T, "test_muon_get_T"},
+    {test_logsumexp, "test_logsumexp"},
+    {test_sno_charge_integral, "test_sno_charge_integral"}
+};
+
+int main(int argc, char **argv)
+{
+    int i;
+    char err[256];
+    int retval = 0;
+    struct tests test;
+
+    for (i = 0; i < sizeof(tests)/sizeof(struct tests); i++) {
+        test = tests[i];
+
+        if (!test.test(err)) {
+            printf("[\033[92mok\033[0m] %s\n", test.name);
+        } else {
+            printf("[\033[91mfail\033[0m] %s: %s\n", test.name, err);
+            retval = 1;
+        }
+    }
+
+    return retval;
+}