diff options
Diffstat (limited to 'test.c')
| -rw-r--r-- | test.c | 162 |
1 files changed, 134 insertions, 28 deletions
@@ -3,6 +3,11 @@ #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); @@ -19,30 +24,30 @@ struct refractive_index_results { double wavelength; double n; } refractive_index_results[] = { - {1000.0, 0 , 226.50, 1.39468}, - {1000.0, 10, 226.50, 1.39439}, - {1000.0, 20, 226.50, 1.39353}, - {1000.0, 30, 226.50, 1.39224}, - {1000.0, 0 , 404.41, 1.34431}, - {1000.0, 10, 404.41, 1.34404}, - {1000.0, 20, 404.41, 1.34329}, - {1000.0, 30, 404.41, 1.34218}, - {1000.0, 0 , 589.00, 1.33447}, - {1000.0, 10, 589.00, 1.33422}, - {1000.0, 20, 589.00, 1.33350}, - {1000.0, 30, 589.00, 1.33243}, - {1000.0, 0 , 632.80, 1.33321}, - {1000.0, 10, 632.80, 1.33296}, - {1000.0, 20, 632.80, 1.33224}, - {1000.0, 30, 632.80, 1.33118}, - {1000.0, 0 , 1013.98, 1.32626}, - {1000.0, 10, 1013.98, 1.32604}, - {1000.0, 20, 1013.98, 1.32537}, - {1000.0, 30, 1013.98, 1.32437}, - {1000.0, 0 , 2325.42, 1.27663}, - {1000.0, 10, 2325.42, 1.27663}, - {1000.0, 20, 2325.42, 1.27627}, - {1000.0, 30, 2325.42, 1.27563}, + {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. @@ -109,7 +114,7 @@ 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) +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. */ @@ -117,7 +122,7 @@ int test_muon_get_E(char *err) /* Assume initial kinetic energy is 1 GeV. */ T = 1000.0; - E = get_E(T,1e-9,1.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)) { @@ -128,7 +133,7 @@ int test_muon_get_E(char *err) /* Assume initial kinetic energy is 1 GeV. */ T = 1000.0; range = get_range(T,1.0); - E = get_E(T,range,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)) { @@ -189,6 +194,38 @@ int test_refractive_index(char *err) 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. */ @@ -221,15 +258,84 @@ int test_solid_angle(char *err) 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_E, "test_muon_get_E"} + {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) |