update rayleigh scattering calculation

This commit updates the optics code to calculate the rayleigh scattering length
using the Einstein-Smoluchowski formula instead of using the effective rayleigh
scattering lengths from the RSPR bank.

4 files changed, 144 insertions, 95 deletions

diff --git a/src/fit.c b/src/fit.c
index 4674e3d..92584aa 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -6044,7 +6044,7 @@ int main(int argc, char **argv)
         goto err;
     }
 
-    if (optics_init(db)) {
+    if (optics_init()) {
         fprintf(stderr, "failed to initialize optics: %s\n", optics_err);
         goto err;
     }
diff --git a/src/likelihood.c b/src/likelihood.c
index b0f6b94..442a8b3 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -235,7 +235,7 @@ static void get_expected_charge_delta_ray(particle *p, double *pos, double *dir,
     get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_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);
+    prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
     constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
 
@@ -281,7 +281,7 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
     get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_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);
+    prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
     constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
 
@@ -730,7 +730,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
     get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_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);
+    prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
     /* Assume the particle is travelling at the speed of light. */
     *t = s/SPEED_OF_LIGHT + l_d2o*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/SPEED_OF_LIGHT;
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;
diff --git a/src/optics.h b/src/optics.h
index 4b41c6b..0e7740b 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -1,4 +1,7 @@
-/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
+/* Library for computing optical properites of the heavy water, light water,
+ * and acrylic in the SNO detector.
+ *
+ * Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
  *
  * This program is free software: you can redistribute it and/or modify it
  * under the terms of the GNU General Public License as published by the Free
@@ -17,21 +20,31 @@
 #ifndef OPTICS_H
 #define OPTICS_H
 
-#include "dict.h"
-
 /* Global error string when optics_init() returns -1. */
 extern char optics_err[256];
 
+/* Index of refraction in water and heavy water averaged over the Cerenkov
+ * spectrum and the PMT quantum efficiency.
+ *
+ * Note: You must call optics_init() before using these! */
 extern double avg_index_h2o, avg_index_d2o;
 
-/* Initialize the optics data by reading in the RSPR bank and precomputing the
- * average absorption and scattering tables. */
-int optics_init(dict *db);
+/* Initialize the optics data by precomputing the average absorption and
+ * scattering tables. */
+int optics_init(void);
 
 /* 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);
+double get_index_snoman_acrylic(double wavelength);
+
+/* Functions for computing the scattering probability as a funcion of
+ * wavelength. */
+double rayint_prob(double wavelength, double n, double isothermal_comp);
+double rayint_prob_d2o(double wavelength);
+double rayint_prob_h2o(double wavelength);
+double rayint_prob_acrylic(double wavelength);
 
 /* Functions for computing the probability that a photon is not absorbed or
  * scattered after a certain distance. */
@@ -40,13 +53,11 @@ double get_fabs_h2o(double x);
 double get_fabs_acrylic(double x);
 double get_fsct_d2o(double x);
 double get_fsct_h2o(double x);
+double get_fsct_acrylic(double x);
 
-/* Functions for computing the absorption and scattering length as a funcion of
- * wavelength. */
+/* Functions for computing the absorption 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