diff options
Diffstat (limited to 'src/test.c')
| -rw-r--r-- | src/test.c | 219 |
1 files changed, 219 insertions, 0 deletions
@@ -19,6 +19,8 @@ #include "likelihood.h" #include "id_particles.h" #include <gsl/gsl_spline2d.h> +#include <gsl/gsl_errno.h> /* for gsl_strerror() */ +#include <gsl/gsl_randist.h> typedef int testFunction(char *err); @@ -970,6 +972,218 @@ err: return 1; } +static double gsl_electron_get_angular_pdf(double cos_theta, void *params) +{ + double alpha = ((double *) params)[0]; + double beta = ((double *) params)[1]; + double mu = ((double *) params)[2]; + return electron_get_angular_pdf(cos_theta,alpha,beta,mu); +} + +int test_electron_get_angular_pdf(char *err) +{ + /* Tests that the electron_get_angular_pdf() function integrates to 1. */ + size_t i; + double params[3]; + gsl_integration_cquad_workspace *w; + gsl_function F; + double result, error; + int status; + size_t nevals; + double T0; + + w = gsl_integration_cquad_workspace_alloc(100); + + F.function = &gsl_electron_get_angular_pdf; + F.params = params; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + T0 = genrand_real2()*1000; + + params[0] = electron_get_angular_distribution_alpha(T0); + params[1] = electron_get_angular_distribution_beta(T0); + params[2] = genrand_real2(); + + status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals); + + if (status) { + sprintf(err, "error integrating electron angular distribution: %s\n", gsl_strerror(status)); + goto err; + } + + if (!isclose(result, 1.0, 0, 1e-3)) { + sprintf(err, "integral of electron_get_angular_pdf() returned %.5g, but expected %.5g", result, 1.0); + goto err; + } + } + + gsl_integration_cquad_workspace_free(w); + + return 0; + +err: + gsl_integration_cquad_workspace_free(w); + + return 1; +} + +int test_gamma_pdf(char *err) +{ + /* Tests the gamma_pdf() function. */ + size_t i; + double x, k, theta, expected; + double result; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + x = genrand_real2(); + k = genrand_real2(); + theta = genrand_real2(); + + result = gamma_pdf(x,k,theta); + + expected = gsl_ran_gamma_pdf(x,k,theta); + + if (!isclose(result, expected, 0, 1e-9)) { + sprintf(err, "gamma_pdf() returned %.5g, but expected %.5g", result, expected); + goto err; + } + } + + return 0; + +err: + return 1; +} + +int test_solid_angle2(char *err) +{ + /* Tests the get_solid_angle()2 function. */ + size_t 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; + double dir[3]; + double L, r0, R; + double expected; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + rand_sphere(pmt); + rand_sphere(n); + r = genrand_real2(); + + dir[0] = pos[0] - pmt[0]; + dir[1] = pos[1] - pmt[1]; + dir[2] = pos[2] - pmt[2]; + + L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]); + R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); + r0 = sqrt(R*R - L*L); + + solid_angle = get_solid_angle2(L/r,r0/r); + + expected = get_solid_angle(pos,pmt,n,r); + + if (!isclose(solid_angle, expected, 1e-4, 0)) { + sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected); + return 1; + } + } + + return 0; +} + +int test_solid_angle_lookup(char *err) +{ + /* Tests the get_solid_angle_lookup() function. */ + size_t 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; + double dir[3]; + double L, r0, R; + double expected; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + rand_sphere(pmt); + MUL(pmt,800); + rand_sphere(n); + r = genrand_real2(); + + dir[0] = pos[0] - pmt[0]; + dir[1] = pos[1] - pmt[1]; + dir[2] = pos[2] - pmt[2]; + + L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]); + R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); + r0 = sqrt(R*R - L*L); + + solid_angle = get_solid_angle_lookup(L/r,r0/r); + + expected = get_solid_angle(pos,pmt,n,r); + + if (!isclose(solid_angle, expected, 1e-2, 0)) { + sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected); + return 1; + } + } + + return 0; +} + +int test_solid_angle_fast(char *err) +{ + /* Tests the get_solid_angle_lookup() function. */ + size_t 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; + double dir[3]; + double L, r0, R; + double expected; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + rand_sphere(pmt); + MUL(pmt,800); + rand_sphere(n); + r = genrand_real2(); + + dir[0] = pos[0] - pmt[0]; + dir[1] = pos[1] - pmt[1]; + dir[2] = pos[2] - pmt[2]; + + L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]); + R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); + r0 = sqrt(R*R - L*L); + + solid_angle = get_solid_angle_fast(pos,pmt,n,r); + + expected = get_solid_angle(pos,pmt,n,r); + + if (!isclose(solid_angle, expected, 1e-2, 0)) { + sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected); + return 1; + } + } + + return 0; +} + struct tests { testFunction *test; char *name; @@ -999,6 +1213,11 @@ struct tests { {test_log_norm, "test_log_norm"}, {test_trapz, "test_trapz"}, {test_interp2d, "test_interp2d"}, + {test_electron_get_angular_pdf, "test_electron_get_angular_pdf"}, + {test_gamma_pdf, "test_gamma_pdf"}, + {test_solid_angle2, "test_solid_angle2"}, + {test_solid_angle_lookup, "test_solid_angle_lookup"}, + {test_solid_angle_fast, "test_solid_angle_fast"}, }; int main(int argc, char **argv) |