1 files changed, 51 insertions, 0 deletions
diff --git a/src/muon.c b/src/muon.c
index 575158b..1cccad6 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -256,6 +256,57 @@ double get_dEdx(double T, double rho)
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
 }
 
+double get_expected_reflected_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, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+
+    z = 1.0;
+
+    SUB(pmt_dir,pmt_pos,pos);
+
+    normalize(pmt_dir);
+
+    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. */
+    cos_theta = DOT(dir,pmt_dir);
+
+    /* Calculate total energy */
+    E = T + MUON_MASS;
+    p = sqrt(E*E - MUON_MASS*MUON_MASS);
+    beta = p/E;
+
+    omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
+
+    R = NORM(pos);
+
+    /* FIXME: I just calculate delta assuming 400 nm light. */
+    wavelength0 = 400.0;
+
+    if (R <= AV_RADIUS) {
+        n = get_index_snoman_d2o(wavelength0);
+    } else {
+        n = get_index_snoman_h2o(wavelength0);
+    }
+
+    if (beta < 1/n) return 0;
+
+    f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
+
+    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
+    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
+    absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
+
+    get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+    l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
+
+    return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+}
+
 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, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;