add functions to get the PMT reflectivity from the PMTR bank

2 files changed, 103 insertions, 7 deletions

diff --git a/src/pmt_response.c b/src/pmt_response.c
index 21e8377..4a61652 100644
--- a/src/pmt_response.c
+++ b/src/pmt_response.c
@@ -14,7 +14,8 @@ static int initialized = 0;
 
 char pmtr_err[256];
 
-static gsl_spline2d *spline;
+static gsl_spline2d *spline_resp;
+static gsl_spline2d *spline_reflec;
 static gsl_interp_accel *xacc;
 static gsl_interp_accel *yacc;
 
@@ -170,6 +171,35 @@ static double resp[NUM_ANGLES*NUM_WAVELENGTHS];
  * 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
+ * and the quantum efficiency. */
+static double weighted_reflec[NUM_ANGLES];
+
+double get_weighted_pmt_reflectivity(double 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);
+}
+
 double get_weighted_pmt_response(double theta)
 {
     /* Returns the probability that a photon is detected as a function of the
@@ -195,6 +225,38 @@ double get_weighted_pmt_response(double theta)
     return interp1d(deg, thetas, weighted_resp, NUM_ANGLES);
 }
 
+double get_pmt_reflectivity(double wavelength, double theta)
+{
+    /* Returns the probability that a photon is reflected as a function of the
+     * photon angle with respect to the PMT normal (in radians) and the
+     * wavelength (in nm). */
+    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;
+
+    if (deg > thetas[NUM_ANGLES-1])
+        deg = thetas[NUM_ANGLES-1];
+
+    if (wavelength < wavelengths[0])
+        wavelength = wavelengths[0];
+
+    if (wavelength > wavelengths[NUM_WAVELENGTHS-1])
+        wavelength = wavelengths[NUM_WAVELENGTHS-1];
+
+    return gsl_spline2d_eval(spline_reflec, deg, wavelength, xacc, xacc);
+}
+
 double get_pmt_response(double wavelength, double theta)
 {
     /* Returns the probability that a photon is detected as a function of the
@@ -238,7 +300,21 @@ double get_pmt_response(double wavelength, double theta)
      * wavelengths the photon won't be detected, it doesn't really matter. */
     if (qe == 0.0) return 0.0;
 
-    return gsl_spline2d_eval(spline, deg, wavelength, xacc, yacc)/qe;
+    return gsl_spline2d_eval(spline_resp, deg, wavelength, xacc, yacc)/qe;
+}
+
+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;
+
+    theta = ((double *) params)[0];
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe*get_pmt_reflectivity(wavelength,theta)/pow(wavelength,2);
 }
 
 static double gsl_pmt_response(double wavelength, void *params)
@@ -278,7 +354,8 @@ int pmt_response_init(dict *db)
     gsl_function F;
     double params[1];
 
-    spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS);
+    spline_resp = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS);
+    spline_reflec = gsl_spline2d_alloc(gsl_interp2d_bilinear, NUM_ANGLES, NUM_WAVELENGTHS);
     xacc = gsl_interp_accel_alloc();
     yacc = gsl_interp_accel_alloc();
 
@@ -291,11 +368,13 @@ int pmt_response_init(dict *db)
 
     for (i = 0; i < NUM_ANGLES; i++) {
         for (j = 0; j < NUM_WAVELENGTHS; j++) {
-            gsl_spline2d_set(spline, resp, i, j, pmtr[30+KPMTR_RESP+i+j*NUM_ANGLES]);
+            gsl_spline2d_set(spline_resp, resp, i, j, pmtr[30+KPMTR_RESP+i+j*NUM_ANGLES]);
+            gsl_spline2d_set(spline_reflec, reflec, i, j, pmtr[30+KPMTR_REFLEC+i+j*NUM_ANGLES]);
         }
     }
 
-    gsl_spline2d_init(spline, thetas, wavelengths, resp, NUM_ANGLES, NUM_WAVELENGTHS);
+    gsl_spline2d_init(spline_resp, thetas, wavelengths, resp, NUM_ANGLES, NUM_WAVELENGTHS);
+    gsl_spline2d_init(spline_reflec, thetas, wavelengths, reflec, NUM_ANGLES, NUM_WAVELENGTHS);
 
     initialized = 1;
 
@@ -312,7 +391,6 @@ int pmt_response_init(dict *db)
     }
 
     F.function = &gsl_pmt_response;
-    F.params = params;
 
     for (i = 0; i < NUM_ANGLES; i++) {
         params[0] = thetas[i]*M_PI/180.0;
@@ -327,6 +405,21 @@ 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;
+
+        status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+        weighted_reflec[i] = result/norm;
+
+        if (status) {
+            fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status));
+            exit(1);
+        }
+    }
+
     gsl_integration_cquad_workspace_free(w);
 
     return 0;
@@ -334,7 +427,8 @@ int pmt_response_init(dict *db)
 
 void pmt_response_free(void)
 {
-    gsl_spline2d_free(spline);
+    gsl_spline2d_free(spline_resp);
+    gsl_spline2d_free(spline_reflec);
     gsl_interp_accel_free(xacc);
     gsl_interp_accel_free(yacc);
 }
diff --git a/src/pmt_response.h b/src/pmt_response.h
index 300f6b2..694ab0a 100644
--- a/src/pmt_response.h
+++ b/src/pmt_response.h
@@ -15,7 +15,9 @@ extern char pmtr_err[256];
 #define KPMTR_RESP 5
 #define KPMTR_REFLEC 4506
 
+double get_weighted_pmt_reflectivity(double theta);
 double get_weighted_pmt_response(double theta);
+double get_pmt_reflectivity(double wavelength, double theta);
 double get_pmt_response(double wavelength, double theta);
 int pmt_response_init(dict *db);
 void pmt_response_free(void);