speed up likelihood function by not calling trapz()

This commit speeds up the likelihood function by integrating the charge along
the track inline instead of creating an array and then calling trapz(). It also
introduces two global variables avg_index_d2o and avg_index_h2o which are the
average indices of refraction for D2O and H2O weighted by the PMT quantum
efficiency and the Cerenkov spectrum.

1 files changed, 48 insertions, 0 deletions

diff --git a/src/optics.c b/src/optics.c
index 830ada9..6b6d761 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -29,6 +29,8 @@
 
 char optics_err[256];
 
+double avg_index_h2o, avg_index_d2o;
+
 /* Absorption coefficients for H2O, D2O, and acrylic as a function of
  * wavelength from SNOMAN.
  *
@@ -425,6 +427,28 @@ static double gsl_absorption_length_snoman_acrylic(double wavelength, void *para
     return qe*exp(-x/get_absorption_length_snoman_acrylic(wavelength))/pow(wavelength,2);
 }
 
+/* Returns the product of the index of refraction in D2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_d2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_d2o(wavelength)/pow(wavelength,2);
+}
+
+/* Returns the product of the index of refraction in H2O multiplied by the
+ * quantum efficiency and the Cerenkov spectrum. */
+static double gsl_index_h2o(double wavelength, void *params)
+{
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_index_snoman_h2o(wavelength)/pow(wavelength,2);
+}
+
 static double gsl_cerenkov(double wavelength, void *params)
 {
     /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */
@@ -458,6 +482,30 @@ int optics_init(dict *db)
         return -1;
     }
 
+    F.function = &gsl_index_d2o;
+    F.params = NULL;
+
+    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;
+    }
+
+    avg_index_d2o = result/norm;
+
+    F.function = &gsl_index_h2o;
+    F.params = NULL;
+
+    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;
+    }
+
+    avg_index_h2o = result/norm;
+
     for (i = 0; i < LEN(fabs_x); i++) {
         fabs_x[i] = xlo + (xhi-xlo)*i/(LEN(fabs_x)-1);
     }