add PMT response

This commit adds code to read in the PMT response from the PMTR bank from
SNOMAN. This file was used for the grey disk model in SNOMAN and was created
using a full 3D simulation of the PMT and concentrator. Since the PMT response
in SNOMAN included the quantum efficiency of the PMT, we have to divide that
out to get just the PMT response independent of the quantum efficiency.

I also updated the likelihood calculation to use the pmt response. Currently
the energy is being fit too high which I think will improve when we update the
solid angle calculation to use the radius of the concentrator instead of the
PMT.

1 files changed, 8 insertions, 3 deletions

diff --git a/src/muon.c b/src/muon.c
index 54a1aa1..d0eb03a 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -13,6 +13,7 @@
 #include "muon.h"
 #include "sno.h"
 #include "scattering.h"
+#include "pmt_response.h"
 
 static int initialized = 0;
 
@@ -226,14 +227,16 @@ double get_dEdx(double T, double rho)
 
 double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
 {
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R;
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt;
 
     z = 1.0;
 
     SUB(pmt_dir,pmt_pos,pos);
     normalize(pmt_dir);
 
-    if (DOT(pmt_dir,pmt_normal) > 0) return 0;
+    cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+
+    if (cos_theta_pmt > 0) return 0;
 
     /* Calculate the cosine of the angle between the track direction and the
      * vector to the PMT. */
@@ -259,6 +262,8 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl
 
     if (beta < 1/n) return 0;
 
+    f = get_weighted_pmt_response(acos(-cos_theta_pmt));
+
     /* FIXME: add angular response and scattering/absorption. */
-    return omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+    return f*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
 }