speed up the likelihood calculation by avoiding calls to acos()

This commit speeds up the likelihood calculation by eliminating most calls to
acos(). This is done by updating the PMT response lookup tables to be as a
function of the cosine of the angle between the photon and the PMT normal
instead of the angle itself.

3 files changed, 72 insertions, 62 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 379ac9a..d29be5e 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, theta_pmt, charge, constant, prob_abs, prob_sct;
+    double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct;
 
     SUB(pmt_dir,pmt_pos,pos);
 
@@ -230,9 +230,8 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir,
 
     omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
 
-    theta_pmt = acos(-cos_theta_pmt);
-    f_reflec = get_weighted_pmt_reflectivity(theta_pmt);
-    f = get_weighted_pmt_response(theta_pmt);
+    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
+    f = get_weighted_pmt_response(-cos_theta_pmt);
 
     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);
@@ -261,7 +260,7 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir,
 
 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, theta_pmt, charge, prob_abs, prob_sct;
+    double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct;
 
     z = 1.0;
 
@@ -289,9 +288,8 @@ static double get_expected_charge(double x, double beta, double theta0, double *
 
     omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
 
-    theta_pmt = acos(-cos_theta_pmt);
-    f_reflec = get_weighted_pmt_reflectivity(theta_pmt);
-    f = get_weighted_pmt_response(theta_pmt);
+    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
+    f = get_weighted_pmt_response(-cos_theta_pmt);
 
     /* Probability that a photon is absorbed. We calculate this by computing:
      *
@@ -545,7 +543,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, l_h2o, l_d2o, f_reflec, charge, prob, frac, prob_abs, prob_sct;
+    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, 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. */
@@ -642,9 +640,8 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     frac = sqrt(2)*n_d2o*x*beta0*theta0;
 
-    theta_pmt = acos(-cos_theta_pmt);
-    f = get_weighted_pmt_response(theta_pmt);
-    f_reflec = get_weighted_pmt_reflectivity(theta_pmt);
+    f = get_weighted_pmt_response(-cos_theta_pmt);
+    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
 
     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);
diff --git a/src/pmt_response.c b/src/pmt_response.c
index cedc341..7d42a39 100644
--- a/src/pmt_response.c
+++ b/src/pmt_response.c
@@ -28,13 +28,18 @@
 
 static int initialized = 0;
 
+/* Global error string set when pmt_response_init() returns -1. */
 char pmtr_err[256];
 
+/* 2D lookup table for the PMT response as a function of angle and wavelength. */
 static gsl_spline2d *spline_resp;
 static gsl_spline2d *spline_reflec;
 static gsl_interp_accel *xacc;
 static gsl_interp_accel *yacc;
 
+/* Angles (in degrees) that the PMT response lookup table is defined for.
+ * Ideally these would have been stored with the table itself, but these are
+ * hardcoded in SNOMAN. */
 static double thetas[NUM_ANGLES] = {
     0.0,
     1.0,
@@ -128,6 +133,9 @@ static double thetas[NUM_ANGLES] = {
     89.0,
 };
 
+/* Wavelengths (in nm) that the PMT response lookup table is defined for.
+ * Ideally these would have been stored with the table itself, but these are
+ * hardcoded in SNOMAN. */
 static double wavelengths[NUM_WAVELENGTHS] = {
     220.0,
     230.0,
@@ -183,62 +191,62 @@ static double wavelengths[NUM_WAVELENGTHS] = {
 
 /* PMT response as a function of angle and wavelength. */
 static double resp[NUM_ANGLES*NUM_WAVELENGTHS];
-/* PMT response as a function of angle weighted by the Cerenkov spectrum and
- * the quantum efficiency. */
-static double weighted_resp[NUM_ANGLES];
-
 /* PMT reflectivity as a function of angle and wavelength. */
 static double reflec[NUM_ANGLES*NUM_WAVELENGTHS];
-/* PMT reflectivity as a function of angle weighted by the Cerenkov spectrum
+
+#define N 1000
+
+/* Lower bound for cos(theta). */
+static double xlo = 0.0;
+/* Upper bound for cos(theta). */
+static double xhi = 1.0;
+
+/* Array holding the values of cos(theta). */
+static double cos_thetas[N];
+
+/* PMT response as a function of cos(theta) weighted by the Cerenkov spectrum
  * and the quantum efficiency. */
-static double weighted_reflec[NUM_ANGLES];
+static double weighted_resp[N];
+/* PMT reflectivity as a function of cos(theta) weighted by the Cerenkov
+ * spectrum and the quantum efficiency. */
+static double weighted_reflec[N];
 
-double get_weighted_pmt_reflectivity(double theta)
+/* Returns the probability that a photon is reflected as a function of the
+ * cosine of the photon angle with respect to the PMT normal.
+ *
+ * Note: The angle should be relative to the negative of the PMT normal, i.e. a
+ * photon incident perpendicular to the front of the PMT should have
+ *
+ *     cos(theta) = 1.
+ *
+ */
+double get_weighted_pmt_reflectivity(double cos_theta)
 {
-    /* Returns the probability that a photon is reflected as a function of the
-     * photon angle with respect to the PMT normal. `theta` should be in
-     * radians. */
-    double deg;
-
     if (!initialized) {
         fprintf(stderr, "pmt response hasn't been initialized!\n");
         exit(1);
     }
 
-    /* Convert to degrees since the PMTR table is in degrees. */
-    deg = theta*180.0/M_PI;
-
-    deg = fmod(deg,180.0);
-
-    if (deg < 0.0)
-        deg += 180.0;
-
-    return interp1d(deg, thetas, weighted_reflec, NUM_ANGLES);
+    return interp1d(cos_theta, cos_thetas, weighted_reflec, LEN(cos_thetas));
 }
 
-double get_weighted_pmt_response(double theta)
+/* Returns the probability that a photon is detected as a function of the
+ * cosine of the photon angle with respect to the PMT normal.
+ *
+ * Note: The angle should be relative to the negative of the PMT normal, i.e. a
+ * photon incident perpendicular to the front of the PMT should have
+ *
+ *     cos(theta) = 1.
+ *
+ * Note: This does *not* include the PMT quantum efficiency. */
+double get_weighted_pmt_response(double cos_theta)
 {
-    /* Returns the probability that a photon is detected as a function of the
-     * photon angle with respect to the PMT normal. `theta` should be in
-     * radians.
-     *
-     * Note: This does *not* include the PMT quantum efficiency. */
-    double deg;
-
     if (!initialized) {
         fprintf(stderr, "pmt response hasn't been initialized!\n");
         exit(1);
     }
 
-    /* Convert to degrees since the PMTR table is in degrees. */
-    deg = theta*180.0/M_PI;
-
-    deg = fmod(deg,180.0);
-
-    if (deg < 0.0)
-        deg += 180.0;
-
-    return interp1d(deg, thetas, weighted_resp, NUM_ANGLES);
+    return interp1d(cos_theta, cos_thetas, weighted_resp, LEN(cos_thetas));
 }
 
 double get_pmt_reflectivity(double wavelength, double theta)
@@ -324,13 +332,13 @@ static double gsl_pmt_reflectivity(double wavelength, void *params)
     /* Returns the probability that a photon is reflected as a function of
      * wavelength for a specific angle weighted by the quantum efficiency and
      * the Cerenkov spectrum. */
-    double qe, theta;
+    double qe, cos_theta;
 
-    theta = ((double *) params)[0];
+    cos_theta = ((double *) params)[0];
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_pmt_reflectivity(wavelength,theta)/pow(wavelength,2);
+    return qe*get_pmt_reflectivity(wavelength,acos(cos_theta))/pow(wavelength,2);
 }
 
 static double gsl_pmt_response(double wavelength, void *params)
@@ -338,13 +346,13 @@ static double gsl_pmt_response(double wavelength, void *params)
     /* Returns the probability that a photon is detected as a function of
      * wavelength for a specific angle weighted by the quantum efficiency and
      * the Cerenkov spectrum. */
-    double qe, theta;
+    double qe, cos_theta;
 
-    theta = ((double *) params)[0];
+    cos_theta = ((double *) params)[0];
 
     qe = get_quantum_efficiency(wavelength);
 
-    return qe*get_pmt_response(wavelength,theta)/pow(wavelength,2);
+    return qe*get_pmt_response(wavelength,acos(cos_theta))/pow(wavelength,2);
 }
 
 static double gsl_cerenkov(double wavelength, void *params)
@@ -406,10 +414,14 @@ int pmt_response_init(dict *db)
         exit(1);
     }
 
+    for (i = 0; i < LEN(cos_thetas); i++) {
+        cos_thetas[i] = xlo + (xhi-xlo)*i/(LEN(cos_thetas)-1);
+    }
+
     F.function = &gsl_pmt_response;
 
-    for (i = 0; i < NUM_ANGLES; i++) {
-        params[0] = thetas[i]*M_PI/180.0;
+    for (i = 0; i < LEN(cos_thetas); i++) {
+        params[0] = cos_thetas[i];
 
         status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
@@ -423,8 +435,8 @@ int pmt_response_init(dict *db)
 
     F.function = &gsl_pmt_reflectivity;
 
-    for (i = 0; i < NUM_ANGLES; i++) {
-        params[0] = thetas[i]*M_PI/180.0;
+    for (i = 0; i < LEN(cos_thetas); i++) {
+        params[0] = cos_thetas[i];
 
         status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
 
diff --git a/src/pmt_response.h b/src/pmt_response.h
index fd23ab2..c87ce0d 100644
--- a/src/pmt_response.h
+++ b/src/pmt_response.h
@@ -19,6 +19,7 @@
 
 #include "dict.h"
 
+/* Global error string set when pmt_response_init() returns -1. */
 extern char pmtr_err[256];
 
 #define NUM_ANGLES 90
@@ -31,8 +32,8 @@ extern char pmtr_err[256];
 #define KPMTR_RESP 4
 #define KPMTR_REFLEC 4504
 
-double get_weighted_pmt_reflectivity(double theta);
-double get_weighted_pmt_response(double theta);
+double get_weighted_pmt_reflectivity(double cos_theta);
+double get_weighted_pmt_response(double cos_theta);
 double get_pmt_reflectivity(double wavelength, double theta);
 double get_pmt_response(double wavelength, double theta);
 int pmt_response_init(dict *db);