From 24c8bcfe7f76b20124e2862ea050f815c0f768e7 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 14 Aug 2018 10:08:27 -0500 Subject: move everything to src directory --- src/test.c | 360 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 360 insertions(+) create mode 100644 src/test.c (limited to 'src/test.c') 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 +#include +#include "optics.h" +#include "muon.h" +#include "mt19937ar.h" +#include +#include "misc.h" +#include "sno_charge.h" +#include + +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; +} -- cgit