diff options
Diffstat (limited to 'src/optics.c')
| -rw-r--r-- | src/optics.c | 200 |
1 files changed, 119 insertions, 81 deletions
diff --git a/src/optics.c b/src/optics.c index 6b6d761..89bee66 100644 --- a/src/optics.c +++ b/src/optics.c @@ -24,29 +24,35 @@ #include <gsl/gsl_spline.h> #include "misc.h" #include "db.h" -#include "dict.h" #include "sno.h" char optics_err[256]; -double avg_index_h2o, avg_index_d2o; +/* Index of refraction in water and heavy water averaged over the Cerenkov + * spectrum and the PMT quantum efficiency. + * + * Set these to 0.0 so that it's obvious if someone forgets to call + * optics_init() before using them. */ +double avg_index_h2o = 0.0, avg_index_d2o = 0.0; /* Absorption coefficients for H2O, D2O, and acrylic as a function of * wavelength from SNOMAN. * * Note: These numbers come from mcprod/media_qoca_d2o_20060110.cmd. */ -static double absorption_coefficient_h2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; -static double absorption_coefficient_h2o[6] = {8.410e-05, 2.880e-06, 1.431e-04, 1.034e-04, 5.239e-04, 2.482e-03}; - -static double rayleigh_coefficient_wavelengths[51]; -static double rayleigh_coefficient_h2o[51]; -static double rayleigh_coefficient_d2o[51]; - static double absorption_coefficient_d2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_d2o[6] = {6.057e-05, 2.680e-05, 3.333e-05, 3.006e-05, 2.773e-05, 2.736e-05}; +static double rayleigh_scl_fac_d2o = 1.000; +static double isothermal_comp_d2o = 4.92e-10; + +static double absorption_coefficient_h2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; +static double absorption_coefficient_h2o[6] = {8.410e-05, 2.880e-06, 1.431e-04, 1.034e-04, 5.239e-04, 2.482e-03}; +static double rayleigh_scl_fac_acrylic = 0.950; +static double isothermal_comp_acrylic = 3.55e-10; static double absorption_coefficient_acrylic_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_acrylic[6] = {5.610e-02, 2.279e-02, 1.204e-02, 7.587e-03, 7.036e-03, 7.068e-03}; +static double rayleigh_scl_fac_h2o = 0.870; +static double isothermal_comp_h2o = 4.78e-10; static int initialized = 0; @@ -76,6 +82,7 @@ static double fabs_h2o[N]; static double fabs_acrylic[N]; static double fsct_d2o[N]; static double fsct_h2o[N]; +static double fsct_acrylic[N]; /* From Table 4 in the paper. */ static double A0 = 0.243905091; @@ -97,6 +104,10 @@ static double RIND_D2O_C1 = 1.302; static double RIND_D2O_C2 = 0.01333; static double RIND_D2O_C3 = 0.32; +static double RIND_ACRYLIC_C1 = 1.452; +static double RIND_ACRYLIC_C2 = 0.02; +static double RIND_ACRYLIC_C3 = 0.32; + /* Wavelength of 1 eV particle in nm. * * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */ @@ -153,6 +164,61 @@ double get_index_snoman_d2o(double wavelength) return RIND_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E); } +double get_index_snoman_acrylic(double wavelength) +{ + /* Returns the index of refraction of SNO acrylic for a given wavelength. + * The wavelength should be in units of nm. + * + * The coefficients come from media.dat in SNOMAN version 5.0294 and the + * formula used to compute the index of refraction comes from the SNOMAN + * companion page for titles MEDA. */ + double E; + + /* Calculate the energy of the photon in eV. */ + E = HC/wavelength; + + return RIND_ACRYLIC_C1 + RIND_ACRYLIC_C2*exp(RIND_ACRYLIC_C3*E); +} + +/* Returns the probability per unit length (1/cm) for Rayleigh scattering in a + * material with an isothermal compressibility of `isothermal_comp` at a given + * wavelength in nm. + * + * The isothermal compressibility should be in units of N^-1 m^2. + * + * This function is mostly copied from snoman's rayint_prob.for. */ +double rayint_prob(double wavelength, double n, double isothermal_comp) +{ + double E; + double confac = 1.53E26; + + /* Calculate the energy of the photon in MeV. */ + E = 1e-6*HC/wavelength; + + return isothermal_comp*confac*pow(E,4)*pow(n*n - 1.0,2)*pow(n*n + 2.0,2); +} + +/* Returns the probability per unit length (1/cm) for Rayleigh scattering in + * D2O at a given wavelength in nm. */ +double rayint_prob_d2o(double wavelength) +{ + return rayleigh_scl_fac_d2o*rayint_prob(wavelength, get_index_snoman_d2o(wavelength), isothermal_comp_d2o); +} + +/* Returns the probability per unit length (1/cm) for Rayleigh scattering in + * H2O at a given wavelength in nm. */ +double rayint_prob_h2o(double wavelength) +{ + return rayleigh_scl_fac_h2o*rayint_prob(wavelength, get_index_snoman_h2o(wavelength), isothermal_comp_h2o); +} + +/* Returns the probability per unit length (1/cm) for Rayleigh scattering in + * acrylic at a given wavelength in nm. */ +double rayint_prob_acrylic(double wavelength) +{ + return rayleigh_scl_fac_acrylic*rayint_prob(wavelength, get_index_snoman_acrylic(wavelength), isothermal_comp_acrylic); +} + /* Returns the average probability that a photon is not absorbed after a * distance `x` in cm in D2O. * @@ -223,6 +289,19 @@ double get_fsct_h2o(double x) return interp1d(x,fabs_x,fsct_h2o,LEN(fabs_x)); } +/* Returns the average probability that a photon is not scattered after a + * distance `x` in cm in SNO acrylic. See get_fabs_d2o() for how this is + * calculated. */ +double get_fsct_acrylic(double x) +{ + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fsct_acrylic,LEN(fabs_x)); +} + /* Returns the absorption length as a function of wavelength in D2O. * * The wavelength should be in nm and the absorption length is in cm. */ @@ -307,85 +386,43 @@ double get_absorption_length_snoman_acrylic(double wavelength) return 1.0/gsl_spline_eval(spline, wavelength, acc); } -/* Returns the Rayleigh scattering length at a particular wavelength in D2O. - * - * The wavelength should be in nm and the returned scattering length is in cm. */ -double get_rayleigh_scattering_length_snoman_d2o(double wavelength) -{ - /* Note: We use GSL splines here instead of interp1d() since we aren't - * guaranteed that the values in the lookup table are evenly spaced. */ - static gsl_spline *spline; - static gsl_interp_accel *acc; - - if (!spline) { - spline = gsl_spline_alloc(gsl_interp_linear, LEN(rayleigh_coefficient_wavelengths)); - gsl_spline_init(spline, rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o, LEN(rayleigh_coefficient_wavelengths)); - acc = gsl_interp_accel_alloc(); - } - - /* If the wavelength is outside the range of the lookup table we just - * return the value at the beginning or end. We do this because by default - * gsl_spline_eval will exit if the wavelength is out of bounds. */ - if (wavelength < rayleigh_coefficient_wavelengths[0]) - wavelength = rayleigh_coefficient_wavelengths[0]; - - if (wavelength > rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1]) - wavelength = rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1]; - - return 1.0/gsl_spline_eval(spline, wavelength, acc); -} - -/* Returns the Rayleigh scattering length at a particular wavelength in H2O. - * - * The wavelength should be in nm and the returned scattering length is in cm. */ -double get_rayleigh_scattering_length_snoman_h2o(double wavelength) +static double gsl_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params) { - /* Returns the Rayleigh scattering length in H2O in cm. */ - static gsl_spline *spline; - static gsl_interp_accel *acc; - - if (!spline) { - spline = gsl_spline_alloc(gsl_interp_linear, LEN(rayleigh_coefficient_wavelengths)); - gsl_spline_init(spline, rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o, LEN(rayleigh_coefficient_wavelengths)); - acc = gsl_interp_accel_alloc(); - } + /* Returns the Rayleigh scattering length in D2O weighted by the quantum + * efficiency and the Cerenkov spectrum. */ + double qe, x; - /* If the wavelength is outside the range of the lookup table we just - * return the value at the beginning or end. We do this because by default - * gsl_spline_eval will exit if the wavelength is out of bounds. */ - if (wavelength < rayleigh_coefficient_wavelengths[0]) - wavelength = rayleigh_coefficient_wavelengths[0]; + x = *((double *) params); - if (wavelength > rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1]) - wavelength = rayleigh_coefficient_wavelengths[LEN(rayleigh_coefficient_wavelengths)-1]; + qe = get_quantum_efficiency(wavelength); - return 1.0/gsl_spline_eval(spline, wavelength, acc); + return qe*exp(-x*rayint_prob_d2o(wavelength))/pow(wavelength,2); } -static double gsl_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params) +static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params) { - /* Returns the Rayleigh scattering length in D2O weighted by the quantum efficiency - * and the Cerenkov spectrum. */ + /* Returns the Rayleigh scattering length in H2O weighted by the quantum + * efficiency and the Cerenkov spectrum. */ double qe, x; x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*exp(-x/get_rayleigh_scattering_length_snoman_d2o(wavelength))/pow(wavelength,2); + return qe*exp(-x*rayint_prob_h2o(wavelength))/pow(wavelength,2); } -static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params) +static double gsl_rayleigh_scattering_length_snoman_acrylic(double wavelength, void *params) { - /* Returns the Rayleigh scattering length in H2O weighted by the quantum efficiency - * and the Cerenkov spectrum. */ + /* Returns the Rayleigh scattering length in acrylic weighted by the + * quantum efficiency and the Cerenkov spectrum. */ double qe, x; x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*exp(-x/get_rayleigh_scattering_length_snoman_h2o(wavelength))/pow(wavelength,2); + return qe*exp(-x*rayint_prob_acrylic(wavelength))/pow(wavelength,2); } static double gsl_absorption_length_snoman_d2o(double wavelength, void *params) @@ -459,7 +496,7 @@ static double gsl_cerenkov(double wavelength, void *params) return qe/pow(wavelength,2); } -int optics_init(dict *db) +int optics_init(void) { int i; double norm; @@ -468,7 +505,6 @@ int optics_init(dict *db) int status; gsl_integration_cquad_workspace *w; gsl_function F; - dbval *rspr; w = gsl_integration_cquad_workspace_alloc(100); @@ -555,19 +591,6 @@ int optics_init(dict *db) fabs_acrylic[i] = result/norm; } - rspr = get_bank(db, "RSPR", 1); - - if (!rspr) { - sprintf(optics_err, "failed to load RSPR bank\n"); - return -1; - } - - for (i = 0; i < 51; i++) { - rayleigh_coefficient_wavelengths[i] = rspr[30+i*3].u32; - rayleigh_coefficient_d2o[i] = rspr[30+i*3+1].f; - rayleigh_coefficient_h2o[i] = rspr[30+i*3+2].f; - } - F.function = &gsl_rayleigh_scattering_length_snoman_d2o; for (i = 0; i < LEN(fabs_x); i++) { @@ -598,6 +621,21 @@ int optics_init(dict *db) fsct_h2o[i] = result/norm; } + F.function = &gsl_rayleigh_scattering_length_snoman_acrylic; + + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fsct_acrylic[i] = result/norm; + } + gsl_integration_cquad_workspace_free(w); initialized = 1; |