fix a bug in the absorption and scattering probabilities

Previously I was computing the fraction of light absorbed and scattered by
calculating an average absorption and scattering length weighted by the
Cerenkov spectrum and the PMT quantum efficiency, which isn't correct since we
should be averaging the absorption and scattering probabilities, not the
absorption and scattering lengths.

This commit fixes this by instead computing the average probability that a
photon is absorbed or scattered as a function of the distance travelled by
integrating the absorption and scattering probabilities over all wavelengths
weighted by the PMT quantum efficiency and the Cerenkov spectrum.

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
diff --git a/src/sno.h b/src/sno.h
index 0a605c9..bf0a1eb 100644
--- a/src/sno.h
+++ b/src/sno.h
@@ -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