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.

1 files changed, 119 insertions, 81 deletions

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