diff options
Diffstat (limited to 'src/pmt_response.c')
| -rw-r--r-- | src/pmt_response.c | 340 |
1 files changed, 340 insertions, 0 deletions
diff --git a/src/pmt_response.c b/src/pmt_response.c new file mode 100644 index 0000000..21e8377 --- /dev/null +++ b/src/pmt_response.c @@ -0,0 +1,340 @@ +#include "pmt_response.h" +#include "dict.h" +#include <stdint.h> /* for uint32_t */ +#include <gsl/gsl_interp2d.h> +#include <gsl/gsl_spline2d.h> +#include <gsl/gsl_integration.h> +#include <stdio.h> /* for sprintf() */ +#include "misc.h" +#include "quantum_efficiency.h" +#include "db.h" +#include <gsl/gsl_errno.h> /* for gsl_strerror() */ + +static int initialized = 0; + +char pmtr_err[256]; + +static gsl_spline2d *spline; +static gsl_interp_accel *xacc; +static gsl_interp_accel *yacc; + +static double thetas[NUM_ANGLES] = { + 0.0, + 1.0, + 2.0, + 3.0, + 4.0, + 5.0, + 6.0, + 7.0, + 8.0, + 9.0, + 10.0, + 11.0, + 12.0, + 13.0, + 14.0, + 15.0, + 16.0, + 17.0, + 18.0, + 19.0, + 20.0, + 21.0, + 22.0, + 23.0, + 24.0, + 25.0, + 26.0, + 27.0, + 28.0, + 29.0, + 30.0, + 31.0, + 32.0, + 33.0, + 34.0, + 35.0, + 36.0, + 37.0, + 38.0, + 39.0, + 40.0, + 41.0, + 42.0, + 43.0, + 44.0, + 45.0, + 46.0, + 47.0, + 48.0, + 49.0, + 50.0, + 51.0, + 52.0, + 53.0, + 54.0, + 55.0, + 56.0, + 57.0, + 58.0, + 59.0, + 60.0, + 61.0, + 62.0, + 63.0, + 64.0, + 65.0, + 66.0, + 67.0, + 68.0, + 69.0, + 70.0, + 71.0, + 72.0, + 73.0, + 74.0, + 75.0, + 76.0, + 77.0, + 78.0, + 79.0, + 80.0, + 81.0, + 82.0, + 83.0, + 84.0, + 85.0, + 86.0, + 87.0, + 88.0, + 89.0, +}; + +static double wavelengths[NUM_WAVELENGTHS] = { + 220.0, + 230.0, + 240.0, + 250.0, + 260.0, + 270.0, + 280.0, + 290.0, + 300.0, + 310.0, + 320.0, + 330.0, + 340.0, + 350.0, + 360.0, + 370.0, + 380.0, + 390.0, + 400.0, + 410.0, + 420.0, + 430.0, + 440.0, + 450.0, + 460.0, + 470.0, + 480.0, + 490.0, + 500.0, + 510.0, + 520.0, + 530.0, + 540.0, + 550.0, + 560.0, + 570.0, + 580.0, + 590.0, + 600.0, + 610.0, + 620.0, + 630.0, + 640.0, + 650.0, + 660.0, + 670.0, + 680.0, + 690.0, + 700.0, + 710.0, +}; + +/* 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]; + +double get_weighted_pmt_response(double 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); +} + +double get_pmt_response(double wavelength, double theta) +{ + /* Returns the probability that a photon is detected as a function of the + * photon angle with respect to the PMT normal (in radians) and the + * wavelength (in nm). + * + * Note: This does *not* include the PMT quantum efficiency. */ + double deg, qe; + + 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; + + if (deg > thetas[NUM_ANGLES-1]) + deg = thetas[NUM_ANGLES-1]; + + if (wavelength < wavelengths[0]) + wavelength = wavelengths[0]; + + if (wavelength > wavelengths[NUM_WAVELENGTHS-1]) + wavelength = wavelengths[NUM_WAVELENGTHS-1]; + + /* The PMTR bank in SNOMAN included the effect of the PMT quantum + * efficiency since it was used for the "grey disk" model. Therefore, + * since we already account for the quantum efficiency it is necessary to + * divide the response by the quantum efficiency. */ + + qe = get_quantum_efficiency(wavelength); + + /* If the quantum efficiency is zero, we have no way of knowing what the + * response should be since it was multiplied by zero. Since for these + * wavelengths the photon won't be detected, it doesn't really matter. */ + if (qe == 0.0) return 0.0; + + return gsl_spline2d_eval(spline, deg, wavelength, xacc, yacc)/qe; +} + +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; + + theta = ((double *) params)[0]; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_pmt_response(wavelength,theta)/pow(wavelength,2); +} + +static double gsl_cerenkov(double wavelength, void *params) +{ + /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe/pow(wavelength,2); +} + +int pmt_response_init(dict *db) +{ + int i, j; + float *pmtr; + double norm; + + double result, error; + size_t nevals; + int status; + gsl_integration_cquad_workspace *w; + gsl_function F; + double params[1]; + + spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS); + xacc = gsl_interp_accel_alloc(); + yacc = gsl_interp_accel_alloc(); + + pmtr = (float *) get_bank(db, "PMTR", 1); + + if (!pmtr) { + sprintf(pmtr_err, "failed to load PMTR bank\n"); + return -1; + } + + for (i = 0; i < NUM_ANGLES; i++) { + for (j = 0; j < NUM_WAVELENGTHS; j++) { + gsl_spline2d_set(spline, resp, i, j, pmtr[30+KPMTR_RESP+i+j*NUM_ANGLES]); + } + } + + gsl_spline2d_init(spline, thetas, wavelengths, resp, NUM_ANGLES, NUM_WAVELENGTHS); + + initialized = 1; + + w = gsl_integration_cquad_workspace_alloc(100); + + F.function = &gsl_cerenkov; + F.params = params; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &norm, &error, &nevals); + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + F.function = &gsl_pmt_response; + F.params = params; + + for (i = 0; i < NUM_ANGLES; i++) { + params[0] = thetas[i]*M_PI/180.0; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + weighted_resp[i] = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + } + + gsl_integration_cquad_workspace_free(w); + + return 0; +} + +void pmt_response_free(void) +{ + gsl_spline2d_free(spline); + gsl_interp_accel_free(xacc); + gsl_interp_accel_free(yacc); +} |