diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 80 | ||||
| -rw-r--r-- | src/optics.c | 600 | ||||
| -rw-r--r-- | src/optics.h | 31 | ||||
| -rw-r--r-- | src/sno.h | 1 |
4 files changed, 420 insertions, 292 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 45274b7..379ac9a 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, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, theta_pmt, charge, constant, prob_abs, prob_sct; SUB(pmt_dir,pmt_pos,pos); @@ -234,15 +234,10 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir, f_reflec = get_weighted_pmt_reflectivity(theta_pmt); f = get_weighted_pmt_response(theta_pmt); - absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); - absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); - absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); - scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); - scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); + 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); - l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; - - constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*omega/(2*M_PI); + constant = get_weighted_quantum_efficiency()*omega/(2*M_PI); /* Note: We assume here that the peak of the angular distribution is at the * Cerenkov angle for a particle with beta = 1. This seems to be an OK @@ -256,19 +251,17 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir, if (p->delta_ray_photons > 0) *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o); - scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); - - *reflected = (f_reflec + 1.0 - scatter)*charge; - *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray); + *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; + *q_indirect_delta_ray = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*(*q_delta_ray) + prob_sct*(*q_delta_ray); - *q_delta_ray *= f*scatter; + *q_delta_ray *= (1.0-prob_abs)*(1.0-prob_sct)*f; - return f*charge*scatter; + return (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } 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, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; + double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, theta_pmt, charge, prob_abs, prob_sct; z = 1.0; @@ -300,21 +293,29 @@ static double get_expected_charge(double x, double beta, double theta0, double * f_reflec = get_weighted_pmt_reflectivity(theta_pmt); f = get_weighted_pmt_response(theta_pmt); - absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); - absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); - absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); - scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); - scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); - - l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; - - charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); + /* Probability that a photon is absorbed. We calculate this by computing: + * + * 1.0 - P(not absorbed in D2O)*P(not absorbed in H2O)*P(not absorbed in acrylic) + * + * since if we worked with the absorption probabilities directly it would + * be more complicated, i.e. + * + * P(absorbed in D2O) + P(absorbed in acrylic|not absorbed in D2O)*P(not absorbed in D2O) + ... + * + */ + prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); + /* Similiar calculation for the probability that a photon is scattered. + * + * Technically we should compute this conditionally on the probability that + * a photon is not absorbed, but since the probability of scattering is + * pretty low, this is expected to be a very small effect. */ + prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o); - scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); + charge = omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); - *reflected = (f_reflec + 1.0 - scatter)*charge; + *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; - return f*charge*scatter; + return (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) @@ -544,7 +545,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, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, f_reflected, charge, prob, frac; + 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; /* First, we find the point along the track where the PMT is at the * Cerenkov angle. */ @@ -643,32 +644,27 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa theta_pmt = acos(-cos_theta_pmt); f = get_weighted_pmt_response(theta_pmt); - f_reflected = get_weighted_pmt_reflectivity(theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(theta_pmt); - absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); - absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); - absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); + 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); get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); - l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; - /* Assume the particle is travelling at the speed of light. */ *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; - double abs_prob = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic); - - charge = abs_prob*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); + charge = n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); /* Add expected number of photons from electromagnetic shower. */ if (p->shower_photons > 0) - charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); + charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); if (p->delta_ray_photons > 0) - charge += abs_prob*get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI); + charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI); - *mu_reflected = f_reflected*charge; + *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; - return f*charge; + return (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } typedef struct betaRootParams { diff --git a/src/optics.c b/src/optics.c index d626bf7..830ada9 100644 --- a/src/optics.c +++ b/src/optics.c @@ -25,11 +25,12 @@ #include "misc.h" #include "db.h" #include "dict.h" +#include "sno.h" char optics_err[256]; -/* Absorption coefficients for H2O and D2O as a function of wavelength from - * SNOMAN. +/* 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}; @@ -47,21 +48,32 @@ static double absorption_coefficient_acrylic[6] = {5.610e-02, 2.279e-02, 1.204e- 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; -/* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum - * efficiency. */ -static double absorption_length_acrylic_weighted = 0.0; -/* D2O Rayleigh scattering length weighted by the Cerenkov spectrum and the - * quantum efficiency. */ -static double rayleigh_scattering_length_d2o_weighted = 0.0; -/* H2O Rayleigh scattering length weighted by the Cerenkov spectrum and the - * quantum efficiency. */ -static double rayleigh_scattering_length_h2o_weighted = 0.0; +/* Number of points to sample when computing the average absorption and + * scattering probabilities. + * + * Note that we use this array to sample for very small distances when + * computing the absorption in the acrylic so (xhi-xlo)/N should be less than + * 10 cm. */ +#define N 10000 + +/* Lower and upper bounds for the distances used to calculate the average + * absorption and scattering probabilities. */ +static const double xlo = 0.0; +static const double xhi = 2*PSUP_RADIUS; + +/* Array of lengths used in the lookup tables for the average absorption and + * scattering probabilities. */ +static double fabs_x[N]; + +/* Average fraction of light which is not absorbed by D2O as a function of + * length weighted by the Cerenkov spectrum and the quantum efficiency. + * + * We calculate the fraction not absorbed instead of absorbed */ +static double fabs_d2o[N]; +static double fabs_h2o[N]; +static double fabs_acrylic[N]; +static double fsct_d2o[N]; +static double fsct_h2o[N]; /* From Table 4 in the paper. */ static double A0 = 0.243905091; @@ -88,59 +100,329 @@ static double RIND_D2O_C3 = 0.32; * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */ static double HC = 1.239841973976e3; +double get_index(double p, double wavelength, double T) +{ + /* Returns the index of refraction of pure water for a given density, + * wavelength, and temperature. The density should be in units of g/cm^3, + * the wavelength in nm, and the temperature in Celsius. + * + * See "Refractive Index of Water and Steam as a function of Wavelength, + * Temperature, and Density" by Schiebener et al. 1990. */ + /* normalize the temperature and pressure */ + wavelength = wavelength/589.0; + T = (T+273.15)/273.15; + + /* first we compute the right hand side of Equation 7 */ + double c = A0 + A1*p + A2*T + A3*pow(wavelength,2)*T + A4/pow(wavelength,2) + A5/(pow(wavelength,2)-pow(UV,2)) + A6/(pow(wavelength,2)-pow(IR,2)) + A7*pow(p,2); + c *= p; + + return sqrt((2*c+1)/(1-c)); +} + +double get_index_snoman_h2o(double wavelength) +{ + /* Returns the index of refraction of light water that SNOMAN uses 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_H2O_C1 + RIND_H2O_C2*exp(RIND_H2O_C3*E); +} + +double get_index_snoman_d2o(double wavelength) +{ + /* Returns the index of refraction of heavy water that SNOMAN uses 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_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E); +} + +/* Returns the average probability that a photon is not absorbed after a + * distance `x` in cm in D2O. + * + * This average probability is found by averaging + * + * e^(-x/l(wavelength)) + * + * over the wavelength range 200-800 nm weighted by the Cerenkov wavelenght + * distribution and the quantum efficiency. */ +double get_fabs_d2o(double x) +{ + if (!initialized) { + fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fabs_d2o,LEN(fabs_x)); +} + +/* Returns the average probability that a photon is not absorbed after a + * distance `x` in cm in H2O. See get_fabs_d2o() for how this is calculated. */ +double get_fabs_h2o(double x) +{ + if (!initialized) { + fprintf(stderr, "average probability hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fabs_h2o,LEN(fabs_x)); +} + +/* Returns the average probability that a photon is not absorbed after a + * distance `x` in cm in the SNO acrylic. See get_fabs_d2o() for how this is + * calculated. + * + * FIXME: Should include a term for the "dark" acrylic. */ +double get_fabs_acrylic(double x) +{ + if (!initialized) { + fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fabs_acrylic,LEN(fabs_x)); +} + +/* Returns the average probability that a photon is not scattered after a + * distance `x` in cm in D2O. See get_fabs_d2o() for how this is calculated. */ +double get_fsct_d2o(double x) +{ + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fsct_d2o,LEN(fabs_x)); +} + +/* Returns the average probability that a photon is not scattered after a + * distance `x` in cm in H2O. See get_fabs_d2o() for how this is calculated. */ +double get_fsct_h2o(double x) +{ + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return interp1d(x,fabs_x,fsct_h2o,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. */ +double get_absorption_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(absorption_coefficient_d2o_wavelengths)); + gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_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 < absorption_coefficient_d2o_wavelengths[0]) + wavelength = absorption_coefficient_d2o_wavelengths[0]; + + if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1]) + wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1]; + + return 1.0/gsl_spline_eval(spline, wavelength, acc); +} + +/* Returns the absorption length as a function of wavelength in H2O. + * + * The wavelength should be in nm and the absorption length is in cm. */ +double get_absorption_length_snoman_h2o(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(absorption_coefficient_h2o_wavelengths)); + gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_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 < absorption_coefficient_h2o_wavelengths[0]) + wavelength = absorption_coefficient_h2o_wavelengths[0]; + + if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1]) + wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1]; + + return 1.0/gsl_spline_eval(spline, wavelength, acc); +} + +/* Returns the absorption length as a function of wavelength in acrylic. + * + * The wavelength should be in nm and the absorption length is in cm. */ +double get_absorption_length_snoman_acrylic(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(absorption_coefficient_acrylic_wavelengths)); + gsl_spline_init(spline, absorption_coefficient_acrylic_wavelengths, absorption_coefficient_acrylic, LEN(absorption_coefficient_acrylic_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 < absorption_coefficient_acrylic_wavelengths[0]) + wavelength = absorption_coefficient_acrylic_wavelengths[0]; + + if (wavelength > absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]) + wavelength = absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]; + + 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) +{ + /* 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(); + } + + /* 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); +} + static double gsl_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params) { /* Returns the Rayleigh scattering length in D2O weighted by the quantum efficiency * and the Cerenkov spectrum. */ - double qe; + double qe, x; + + x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*get_rayleigh_scattering_length_snoman_d2o(wavelength)/pow(wavelength,2); + return qe*exp(-x/get_rayleigh_scattering_length_snoman_d2o(wavelength))/pow(wavelength,2); } static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params) { /* Returns the Rayleigh scattering length in H2O weighted by the quantum efficiency * and the Cerenkov spectrum. */ - double qe; + double qe, x; + + x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*get_rayleigh_scattering_length_snoman_h2o(wavelength)/pow(wavelength,2); + return qe*exp(-x/get_rayleigh_scattering_length_snoman_h2o(wavelength))/pow(wavelength,2); } 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; + double qe, x; + + x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*get_absorption_length_snoman_d2o(wavelength)/pow(wavelength,2); + return qe*exp(-x/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; + double qe, x; + + x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2); + return qe*exp(-x/get_absorption_length_snoman_h2o(wavelength))/pow(wavelength,2); } static double gsl_absorption_length_snoman_acrylic(double wavelength, void *params) { /* Returns the absorption length in acrylic weighted by the quantum * efficiency and the Cerenkov spectrum. */ - double qe; + double qe, x; + + x = *((double *) params); qe = get_quantum_efficiency(wavelength); - return qe*get_absorption_length_snoman_acrylic(wavelength)/pow(wavelength,2); + return qe*exp(-x/get_absorption_length_snoman_acrylic(wavelength))/pow(wavelength,2); } static double gsl_cerenkov(double wavelength, void *params) @@ -172,41 +454,57 @@ int optics_init(dict *db) 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); + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + for (i = 0; i < LEN(fabs_x); i++) { + fabs_x[i] = xlo + (xhi-xlo)*i/(LEN(fabs_x)-1); } F.function = &gsl_absorption_length_snoman_d2o; - status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; - absorption_length_d2o_weighted = result/norm; + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); - if (status) { - fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); - exit(1); + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fabs_d2o[i] = result/norm; } F.function = &gsl_absorption_length_snoman_h2o; - status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; - absorption_length_h2o_weighted = result/norm; + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); - if (status) { - fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); - exit(1); + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fabs_h2o[i] = result/norm; } F.function = &gsl_absorption_length_snoman_acrylic; - status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; - absorption_length_acrylic_weighted = result/norm; + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); - if (status) { - fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); - exit(1); + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fabs_acrylic[i] = result/norm; } rspr = get_bank(db, "RSPR", 1); @@ -224,24 +522,32 @@ int optics_init(dict *db) F.function = &gsl_rayleigh_scattering_length_snoman_d2o; - status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; - rayleigh_scattering_length_d2o_weighted = result/norm; + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); - if (status) { - fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); - exit(1); + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fsct_d2o[i] = result/norm; } F.function = &gsl_rayleigh_scattering_length_snoman_h2o; - status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + for (i = 0; i < LEN(fabs_x); i++) { + F.params = fabs_x+i; - rayleigh_scattering_length_h2o_weighted = result/norm; + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); - if (status) { - fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); - exit(1); + if (status) { + sprintf(optics_err, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + return -1; + } + + fsct_h2o[i] = result/norm; } gsl_integration_cquad_workspace_free(w); @@ -250,189 +556,3 @@ int optics_init(dict *db) 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_rayleigh_scattering_length_snoman_d2o(double wavelength) -{ - /* Returns the Rayleigh scattering length in D2O in cm. */ - return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o,LEN(rayleigh_coefficient_wavelengths)); -} - -double get_rayleigh_scattering_length_snoman_h2o(double wavelength) -{ - /* Returns the Rayleigh scattering length in H2O in cm. */ - return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o,LEN(rayleigh_coefficient_wavelengths)); -} - -double get_weighted_rayleigh_scattering_length_snoman_d2o(void) -{ - /* Returns the weighted Rayleigh scattering length in D2O in cm. */ - - if (!initialized) { - fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); - exit(1); - } - - return rayleigh_scattering_length_d2o_weighted; -} - -double get_weighted_rayleigh_scattering_length_snoman_h2o(void) -{ - /* Returns the weighted Rayleigh scattering length in H2O in cm. */ - - if (!initialized) { - fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); - exit(1); - } - - return rayleigh_scattering_length_h2o_weighted; -} - -double get_absorption_length_snoman_h2o(double wavelength) -{ - /* Returns the absorption length in H2O in cm. */ - static gsl_spline *spline; - static gsl_interp_accel *acc; - - if (!spline) { - spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_h2o_wavelengths)); - gsl_spline_init(spline, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, LEN(absorption_coefficient_h2o_wavelengths)); - acc = gsl_interp_accel_alloc(); - } - - if (wavelength < absorption_coefficient_h2o_wavelengths[0]) - wavelength = absorption_coefficient_h2o_wavelengths[0]; - - if (wavelength > absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1]) - wavelength = absorption_coefficient_h2o_wavelengths[LEN(absorption_coefficient_h2o_wavelengths)-1]; - - return 1.0/gsl_spline_eval(spline, wavelength, acc); -} - -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. */ - static gsl_spline *spline; - static gsl_interp_accel *acc; - - if (!spline) { - spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_d2o_wavelengths)); - gsl_spline_init(spline, absorption_coefficient_d2o_wavelengths, absorption_coefficient_d2o, LEN(absorption_coefficient_d2o_wavelengths)); - acc = gsl_interp_accel_alloc(); - } - - if (wavelength < absorption_coefficient_d2o_wavelengths[0]) - wavelength = absorption_coefficient_d2o_wavelengths[0]; - - if (wavelength > absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1]) - wavelength = absorption_coefficient_d2o_wavelengths[LEN(absorption_coefficient_d2o_wavelengths)-1]; - - return 1.0/gsl_spline_eval(spline, wavelength, acc); -} - -double get_weighted_absorption_length_snoman_acrylic(void) -{ - /* Returns the weighted absorption length in acrylic in cm. */ - - if (!initialized) { - fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); - exit(1); - } - - return absorption_length_acrylic_weighted; -} - -double get_absorption_length_snoman_acrylic(double wavelength) -{ - /* Returns the absorption length in acrylic in cm. */ - static gsl_spline *spline; - static gsl_interp_accel *acc; - - if (!spline) { - spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_acrylic_wavelengths)); - gsl_spline_init(spline, absorption_coefficient_acrylic_wavelengths, absorption_coefficient_acrylic, LEN(absorption_coefficient_acrylic_wavelengths)); - acc = gsl_interp_accel_alloc(); - } - - if (wavelength < absorption_coefficient_acrylic_wavelengths[0]) - wavelength = absorption_coefficient_acrylic_wavelengths[0]; - - if (wavelength > absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]) - wavelength = absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]; - - return 1.0/gsl_spline_eval(spline, wavelength, acc); -} - -double get_index(double p, double wavelength, double T) -{ - /* Returns the index of refraction of pure water for a given density, - * wavelength, and temperature. The density should be in units of g/cm^3, - * the wavelength in nm, and the temperature in Celsius. - * - * See "Refractive Index of Water and Steam as a function of Wavelength, - * Temperature, and Density" by Schiebener et al. 1990. */ - /* normalize the temperature and pressure */ - wavelength = wavelength/589.0; - T = (T+273.15)/273.15; - - /* first we compute the right hand side of Equation 7 */ - double c = A0 + A1*p + A2*T + A3*pow(wavelength,2)*T + A4/pow(wavelength,2) + A5/(pow(wavelength,2)-pow(UV,2)) + A6/(pow(wavelength,2)-pow(IR,2)) + A7*pow(p,2); - c *= p; - - return sqrt((2*c+1)/(1-c)); -} - -double get_index_snoman_h2o(double wavelength) -{ - /* Returns the index of refraction of light water that SNOMAN uses 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_H2O_C1 + RIND_H2O_C2*exp(RIND_H2O_C3*E); -} - -double get_index_snoman_d2o(double wavelength) -{ - /* Returns the index of refraction of heavy water that SNOMAN uses 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_D2O_C1 + RIND_D2O_C2*exp(RIND_D2O_C3*E); -} diff --git a/src/optics.h b/src/optics.h index 19db2b9..70ff53b 100644 --- a/src/optics.h +++ b/src/optics.h @@ -19,21 +19,32 @@ #include "dict.h" +/* Global error string when optics_init() returns -1. */ extern char optics_err[256]; +/* Initialize the optics data by reading in the RSPR bank and precomputing the + * average absorption and scattering tables. */ int optics_init(dict *db); -double get_rayleigh_scattering_length_snoman_d2o(double wavelength); -double get_rayleigh_scattering_length_snoman_h2o(double wavelength); -double get_weighted_rayleigh_scattering_length_snoman_d2o(void); -double get_weighted_rayleigh_scattering_length_snoman_h2o(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_weighted_absorption_length_snoman_acrylic(void); -double get_absorption_length_snoman_acrylic(double wavelength); + +/* Functions for computing the index of refraction. */ double get_index(double p, double wavelength, double T); double get_index_snoman_h2o(double wavelength); double get_index_snoman_d2o(double wavelength); +/* Functions for computing the probability that a photon is not absorbed or + * scattered after a certain distance. */ +double get_fabs_d2o(double x); +double get_fabs_h2o(double x); +double get_fabs_acrylic(double x); +double get_fsct_d2o(double x); +double get_fsct_h2o(double x); + +/* Functions for computing the absorption and scattering length as a funcion of + * wavelength. */ +double get_absorption_length_snoman_d2o(double wavelength); +double get_absorption_length_snoman_h2o(double wavelength); +double get_absorption_length_snoman_acrylic(double wavelength); +double get_rayleigh_scattering_length_snoman_d2o(double wavelength); +double get_rayleigh_scattering_length_snoman_h2o(double wavelength); + #endif @@ -31,6 +31,7 @@ #define AV_RADIUS_INNER 600.50 #define AV_RADIUS_OUTER 606.00 #define AV_RADIUS 603.25 +#define AV_THICKNESS 6.5 /* FIXME: is this right? */ #define PSUP_RADIUS 800.0 |