diff options
| -rw-r--r-- | src/fit.c | 1 | ||||
| -rw-r--r-- | src/muon.c | 4 | ||||
| -rw-r--r-- | src/optics.c | 119 | ||||
| -rw-r--r-- | src/optics.h | 3 |
4 files changed, 125 insertions, 2 deletions
@@ -5432,6 +5432,7 @@ int main(int argc, char **argv) init_interpolation(); init_charge(); + optics_init(); dict *db = db_init(); if (load_file(db, "DQXX_0000010000.dat")) { @@ -298,8 +298,8 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl f = get_weighted_pmt_response(acos(-cos_theta_pmt)); - absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0); - absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0); + absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); + absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); diff --git a/src/optics.c b/src/optics.c index af89802..b01ce7d 100644 --- a/src/optics.c +++ b/src/optics.c @@ -1,6 +1,10 @@ #include <math.h> #include "optics.h" #include "misc.h" +#include "quantum_efficiency.h" +#include <stdio.h> /* for fprintf() */ +#include <gsl/gsl_integration.h> +#include <gsl/gsl_errno.h> /* for gsl_strerror() */ /* Absorption coefficients for H2O and D2O as a function of wavelength from * SNOMAN. @@ -14,6 +18,15 @@ static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e- static double absorption_coefficient_d2o_wavelengths[7] = {254.0, 313.0, 366.0, 406.0, 436.0, 548.0, 578.0}; static double absorption_coefficient_d2o[7] = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5}; +static int initialized = 0; + +/* D2O absorption length weighted by the Cerenkov spectrum and the quantum + * efficiency. */ +static double absorption_length_d2o_weighted = 0.0; +/* H2O absorption length weighted by the Cerenkov spectrum and the quantum + * efficiency. */ +static double absorption_length_h2o_weighted = 0.0; + /* From Table 4 in the paper. */ static double A0 = 0.243905091; static double A1 = 9.53518094e-3; @@ -39,12 +52,118 @@ static double RIND_D2O_C3 = 0.32; * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */ static double HC = 1.239841973976e3; +static double gsl_absorption_length_snoman_d2o(double wavelength, void *params) +{ + /* Returns the absorption length in D2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_absorption_length_snoman_d2o(wavelength)/pow(wavelength,2); +} + +static double gsl_absorption_length_snoman_h2o(double wavelength, void *params) +{ + /* Returns the absorption length in H2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_absorption_length_snoman_h2o(wavelength)/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 optics_init(void) +{ + double norm; + double result, error; + size_t nevals; + int status; + gsl_integration_cquad_workspace *w; + gsl_function F; + + w = gsl_integration_cquad_workspace_alloc(100); + + F.function = &gsl_cerenkov; + F.params = NULL; + + 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_absorption_length_snoman_d2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + absorption_length_d2o_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + F.function = &gsl_absorption_length_snoman_h2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + absorption_length_h2o_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + gsl_integration_cquad_workspace_free(w); + + initialized = 1; + + return 0; +} + +double get_weighted_absorption_length_snoman_h2o(void) +{ + /* Returns the weighted absorption length in H2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); + exit(1); + } + + return absorption_length_h2o_weighted; +} + double get_absorption_length_snoman_h2o(double wavelength) { /* Returns the absorption length in H2O in cm. */ return 1.0/interp1d(wavelength, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, sizeof(absorption_coefficient_h2o_wavelengths)/sizeof(absorption_coefficient_h2o_wavelengths[0])); } +double get_weighted_absorption_length_snoman_d2o(void) +{ + /* Returns the weighted absorption length in D2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); + exit(1); + } + + return absorption_length_d2o_weighted; +} + double get_absorption_length_snoman_d2o(double wavelength) { /* Returns the absorption length in D2O in cm. */ diff --git a/src/optics.h b/src/optics.h index 13b3e0b..fdeefdd 100644 --- a/src/optics.h +++ b/src/optics.h @@ -1,7 +1,10 @@ #ifndef OPTICS_H #define OPTICS_H +int optics_init(void); +double get_weighted_absorption_length_snoman_h2o(void); double get_absorption_length_snoman_h2o(double wavelength); +double get_weighted_absorption_length_snoman_d2o(void); double get_absorption_length_snoman_d2o(double wavelength); double get_index(double p, double wavelength, double T); double get_index_snoman_h2o(double wavelength); |