1 files changed, 60 insertions, 48 deletions

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);