add absorption lengths for D2O and H2O weighted by the Cerenkov spectrum and the quantum efficiency

4 files changed, 125 insertions, 2 deletions

diff --git a/src/fit.c b/src/fit.c
index 293c19c..09da3a0 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5432,6 +5432,7 @@ int main(int argc, char **argv)
 
     init_interpolation();
     init_charge();
+    optics_init();
     dict *db = db_init();
 
     if (load_file(db, "DQXX_0000010000.dat")) {
diff --git a/src/muon.c b/src/muon.c
index 9114dda..ae9c08c 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -298,8 +298,8 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl
 
     f = get_weighted_pmt_response(acos(-cos_theta_pmt));
 
-    absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
-    absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
+    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
+    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
 
     get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
 
diff --git a/src/optics.c b/src/optics.c
index af89802..b01ce7d 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -1,6 +1,10 @@
 #include <math.h>
 #include "optics.h"
 #include "misc.h"
+#include "quantum_efficiency.h"
+#include <stdio.h> /* for fprintf() */
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_errno.h> /* for gsl_strerror() */
 
 /* Absorption coefficients for H2O and D2O as a function of wavelength from
  * SNOMAN.
@@ -14,6 +18,15 @@ static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e-
 static double absorption_coefficient_d2o_wavelengths[7] = {254.0, 313.0, 366.0, 406.0, 436.0, 548.0, 578.0};
 static double absorption_coefficient_d2o[7] = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5};
 
+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;
+
 /* From Table 4 in the paper. */
 static double A0 = 0.243905091;
 static double A1 = 9.53518094e-3;
@@ -39,12 +52,118 @@ static double RIND_D2O_C3 = 0.32;
  * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
 static double HC = 1.239841973976e3;
 
+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;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*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;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2);
+}
+
+static double gsl_cerenkov(double wavelength, void *params)
+{
+    /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe/pow(wavelength,2);
+}
+
+int optics_init(void)
+{
+    double norm;
+    double result, error;
+    size_t nevals;
+    int status;
+    gsl_integration_cquad_workspace *w;
+    gsl_function F;
+
+    w = gsl_integration_cquad_workspace_alloc(100);
+
+    F.function = &gsl_cerenkov;
+    F.params = NULL;
+
+    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);
+    }
+
+    F.function = &gsl_absorption_length_snoman_d2o;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    absorption_length_d2o_weighted = result/norm;
+
+    if (status) {
+        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    F.function = &gsl_absorption_length_snoman_h2o;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    absorption_length_h2o_weighted = result/norm;
+
+    if (status) {
+        fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    gsl_integration_cquad_workspace_free(w);
+
+    initialized = 1;
+
+    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_absorption_length_snoman_h2o(double wavelength)
 {
     /* Returns the absorption length in H2O in cm. */
     return 1.0/interp1d(wavelength, absorption_coefficient_h2o_wavelengths, absorption_coefficient_h2o, sizeof(absorption_coefficient_h2o_wavelengths)/sizeof(absorption_coefficient_h2o_wavelengths[0]));
 }
 
+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. */
diff --git a/src/optics.h b/src/optics.h
index 13b3e0b..fdeefdd 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -1,7 +1,10 @@
 #ifndef OPTICS_H
 #define OPTICS_H
 
+int optics_init(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_index(double p, double wavelength, double T);
 double get_index_snoman_h2o(double wavelength);