From 8447870e721dd738bce12b45e732c9cc01dc2595 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Sun, 11 Nov 2018 13:22:18 -0600 Subject: update likelihood function to fit electrons! To characterize the angular distribution of photons from an electromagnetic shower I came up with the following functional form: f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) and fit this to data simulated using RAT-PAC at several different energies. I then fit the alpha and beta coefficients as a function of energy to the functional form: alpha = c0 + c1/log(c2*T0 + c3) beta = c0 + c1/log(c2*T0 + c3). where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3 are parameters which I fit. The longitudinal distribution of the photons generated from an electromagnetic shower is described by a gamma distribution: f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a). This parameterization comes from the PDG "Passage of particles through matter" section 32.5. I also fit the data from my RAT-PAC simulation, but currently I am not using it, and instead using a simpler form to calculate the coefficients from the PDG (although I estimated the b parameter from the RAT-PAC data). I also sped up the calculation of the solid angle by making a lookup table since it was taking a significant fraction of the time to compute the likelihood function. --- src/test.c | 219 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 219 insertions(+) (limited to 'src/test.c') diff --git a/src/test.c b/src/test.c index f8cedb9..7179f6b 100644 --- a/src/test.c +++ b/src/test.c @@ -19,6 +19,8 @@ #include "likelihood.h" #include "id_particles.h" #include +#include /* for gsl_strerror() */ +#include 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) -- cgit