diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 21 | ||||
| -rw-r--r-- | src/pmt_response.c | 108 | ||||
| -rw-r--r-- | src/pmt_response.h | 5 |
3 files changed, 72 insertions, 62 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 379ac9a..d29be5e 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -210,7 +210,7 @@ void particle_free(particle *p) static double get_expected_charge_shower(particle *p, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double *q_delta_ray, double *q_indirect_delta_ray) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, theta_pmt, charge, constant, prob_abs, prob_sct; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct; SUB(pmt_dir,pmt_pos,pos); @@ -230,9 +230,8 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir, omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); - theta_pmt = acos(-cos_theta_pmt); - f_reflec = get_weighted_pmt_reflectivity(theta_pmt); - f = get_weighted_pmt_response(theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); + f = get_weighted_pmt_response(-cos_theta_pmt); prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o); @@ -261,7 +260,7 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir, static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o) { - double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, theta_pmt, charge, prob_abs, prob_sct; + double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct; z = 1.0; @@ -289,9 +288,8 @@ static double get_expected_charge(double x, double beta, double theta0, double * omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); - theta_pmt = acos(-cos_theta_pmt); - f_reflec = get_weighted_pmt_reflectivity(theta_pmt); - f = get_weighted_pmt_response(theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); + f = get_weighted_pmt_response(-cos_theta_pmt); /* Probability that a photon is absorbed. We calculate this by computing: * @@ -545,7 +543,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa * * `smax` is currently calculated as the point where the particle velocity * drops to 0.8 times the speed of light. */ - double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, f, cos_theta_pmt, theta_pmt, l_h2o, l_d2o, f_reflec, charge, prob, frac, prob_abs, prob_sct; + double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, f, cos_theta_pmt, l_h2o, l_d2o, f_reflec, charge, prob, frac, prob_abs, prob_sct; /* First, we find the point along the track where the PMT is at the * Cerenkov angle. */ @@ -642,9 +640,8 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa frac = sqrt(2)*n_d2o*x*beta0*theta0; - theta_pmt = acos(-cos_theta_pmt); - f = get_weighted_pmt_response(theta_pmt); - f_reflec = get_weighted_pmt_reflectivity(theta_pmt); + f = get_weighted_pmt_response(-cos_theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o); diff --git a/src/pmt_response.c b/src/pmt_response.c index cedc341..7d42a39 100644 --- a/src/pmt_response.c +++ b/src/pmt_response.c @@ -28,13 +28,18 @@ static int initialized = 0; +/* Global error string set when pmt_response_init() returns -1. */ char pmtr_err[256]; +/* 2D lookup table for the PMT response as a function of angle and wavelength. */ static gsl_spline2d *spline_resp; static gsl_spline2d *spline_reflec; static gsl_interp_accel *xacc; static gsl_interp_accel *yacc; +/* Angles (in degrees) that the PMT response lookup table is defined for. + * Ideally these would have been stored with the table itself, but these are + * hardcoded in SNOMAN. */ static double thetas[NUM_ANGLES] = { 0.0, 1.0, @@ -128,6 +133,9 @@ static double thetas[NUM_ANGLES] = { 89.0, }; +/* Wavelengths (in nm) that the PMT response lookup table is defined for. + * Ideally these would have been stored with the table itself, but these are + * hardcoded in SNOMAN. */ static double wavelengths[NUM_WAVELENGTHS] = { 220.0, 230.0, @@ -183,62 +191,62 @@ static double wavelengths[NUM_WAVELENGTHS] = { /* PMT response as a function of angle and wavelength. */ static double resp[NUM_ANGLES*NUM_WAVELENGTHS]; -/* PMT response as a function of angle weighted by the Cerenkov spectrum and - * the quantum efficiency. */ -static double weighted_resp[NUM_ANGLES]; - /* PMT reflectivity as a function of angle and wavelength. */ static double reflec[NUM_ANGLES*NUM_WAVELENGTHS]; -/* PMT reflectivity as a function of angle weighted by the Cerenkov spectrum + +#define N 1000 + +/* Lower bound for cos(theta). */ +static double xlo = 0.0; +/* Upper bound for cos(theta). */ +static double xhi = 1.0; + +/* Array holding the values of cos(theta). */ +static double cos_thetas[N]; + +/* PMT response as a function of cos(theta) weighted by the Cerenkov spectrum * and the quantum efficiency. */ -static double weighted_reflec[NUM_ANGLES]; +static double weighted_resp[N]; +/* PMT reflectivity as a function of cos(theta) weighted by the Cerenkov + * spectrum and the quantum efficiency. */ +static double weighted_reflec[N]; -double get_weighted_pmt_reflectivity(double theta) +/* Returns the probability that a photon is reflected as a function of the + * cosine of the photon angle with respect to the PMT normal. + * + * Note: The angle should be relative to the negative of the PMT normal, i.e. a + * photon incident perpendicular to the front of the PMT should have + * + * cos(theta) = 1. + * + */ +double get_weighted_pmt_reflectivity(double cos_theta) { - /* Returns the probability that a photon is reflected as a function of the - * photon angle with respect to the PMT normal. `theta` should be in - * radians. */ - double deg; - if (!initialized) { fprintf(stderr, "pmt response hasn't been initialized!\n"); exit(1); } - /* Convert to degrees since the PMTR table is in degrees. */ - deg = theta*180.0/M_PI; - - deg = fmod(deg,180.0); - - if (deg < 0.0) - deg += 180.0; - - return interp1d(deg, thetas, weighted_reflec, NUM_ANGLES); + return interp1d(cos_theta, cos_thetas, weighted_reflec, LEN(cos_thetas)); } -double get_weighted_pmt_response(double theta) +/* Returns the probability that a photon is detected as a function of the + * cosine of the photon angle with respect to the PMT normal. + * + * Note: The angle should be relative to the negative of the PMT normal, i.e. a + * photon incident perpendicular to the front of the PMT should have + * + * cos(theta) = 1. + * + * Note: This does *not* include the PMT quantum efficiency. */ +double get_weighted_pmt_response(double cos_theta) { - /* Returns the probability that a photon is detected as a function of the - * photon angle with respect to the PMT normal. `theta` should be in - * radians. - * - * Note: This does *not* include the PMT quantum efficiency. */ - double deg; - if (!initialized) { fprintf(stderr, "pmt response hasn't been initialized!\n"); exit(1); } - /* Convert to degrees since the PMTR table is in degrees. */ - deg = theta*180.0/M_PI; - - deg = fmod(deg,180.0); - - if (deg < 0.0) - deg += 180.0; - - return interp1d(deg, thetas, weighted_resp, NUM_ANGLES); + return interp1d(cos_theta, cos_thetas, weighted_resp, LEN(cos_thetas)); } double get_pmt_reflectivity(double wavelength, double theta) @@ -324,13 +332,13 @@ static double gsl_pmt_reflectivity(double wavelength, void *params) /* Returns the probability that a photon is reflected as a function of * wavelength for a specific angle weighted by the quantum efficiency and * the Cerenkov spectrum. */ - double qe, theta; + double qe, cos_theta; - theta = ((double *) params)[0]; + cos_theta = ((double *) params)[0]; qe = get_quantum_efficiency(wavelength); - return qe*get_pmt_reflectivity(wavelength,theta)/pow(wavelength,2); + return qe*get_pmt_reflectivity(wavelength,acos(cos_theta))/pow(wavelength,2); } static double gsl_pmt_response(double wavelength, void *params) @@ -338,13 +346,13 @@ static double gsl_pmt_response(double wavelength, void *params) /* Returns the probability that a photon is detected as a function of * wavelength for a specific angle weighted by the quantum efficiency and * the Cerenkov spectrum. */ - double qe, theta; + double qe, cos_theta; - theta = ((double *) params)[0]; + cos_theta = ((double *) params)[0]; qe = get_quantum_efficiency(wavelength); - return qe*get_pmt_response(wavelength,theta)/pow(wavelength,2); + return qe*get_pmt_response(wavelength,acos(cos_theta))/pow(wavelength,2); } static double gsl_cerenkov(double wavelength, void *params) @@ -406,10 +414,14 @@ int pmt_response_init(dict *db) exit(1); } + for (i = 0; i < LEN(cos_thetas); i++) { + cos_thetas[i] = xlo + (xhi-xlo)*i/(LEN(cos_thetas)-1); + } + F.function = &gsl_pmt_response; - for (i = 0; i < NUM_ANGLES; i++) { - params[0] = thetas[i]*M_PI/180.0; + for (i = 0; i < LEN(cos_thetas); i++) { + params[0] = cos_thetas[i]; status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); @@ -423,8 +435,8 @@ int pmt_response_init(dict *db) F.function = &gsl_pmt_reflectivity; - for (i = 0; i < NUM_ANGLES; i++) { - params[0] = thetas[i]*M_PI/180.0; + for (i = 0; i < LEN(cos_thetas); i++) { + params[0] = cos_thetas[i]; status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); diff --git a/src/pmt_response.h b/src/pmt_response.h index fd23ab2..c87ce0d 100644 --- a/src/pmt_response.h +++ b/src/pmt_response.h @@ -19,6 +19,7 @@ #include "dict.h" +/* Global error string set when pmt_response_init() returns -1. */ extern char pmtr_err[256]; #define NUM_ANGLES 90 @@ -31,8 +32,8 @@ extern char pmtr_err[256]; #define KPMTR_RESP 4 #define KPMTR_REFLEC 4504 -double get_weighted_pmt_reflectivity(double theta); -double get_weighted_pmt_response(double theta); +double get_weighted_pmt_reflectivity(double cos_theta); +double get_weighted_pmt_response(double cos_theta); double get_pmt_reflectivity(double wavelength, double theta); double get_pmt_response(double wavelength, double theta); int pmt_response_init(dict *db); |