add absorption length

This commit adds the absorption length to the likelihood calculation. For now
I'm just using a single number independent of wavelength. I should update this
in the future to actually use the absorption lengths as measured by SNO and
then calculate an overall absorption length weighted by the Cerenkov spectrum
and the PMT quantum efficiency.

3 files changed, 14 insertions, 2 deletions

diff --git a/src/muon.c b/src/muon.c
index d0eb03a..186363e 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -227,11 +227,14 @@ 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, f, cos_theta_pmt;
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length, distance;
 
     z = 1.0;
 
     SUB(pmt_dir,pmt_pos,pos);
+
+    distance = NORM(pmt_dir);
+
     normalize(pmt_dir);
 
     cos_theta_pmt = DOT(pmt_dir,pmt_normal);
@@ -264,6 +267,8 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl
 
     f = get_weighted_pmt_response(acos(-cos_theta_pmt));
 
+    absorption_length = get_absorption_length_snoman_d2o(wavelength0);
+
     /* FIXME: add angular response and scattering/absorption. */
-    return f*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+    return f*exp(-distance/absorption_length)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
 }
diff --git a/src/optics.c b/src/optics.c
index dd9b82d..99c2a51 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -26,6 +26,12 @@ static double RIND_D2O_C3 = 0.32;
  * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
 static double HC = 1.239841973976e3;
 
+double get_absorption_length_snoman_d2o(double wavelength)
+{
+    /* Returns the absorption length in D2O in cm. */
+    return 1e4;
+}
+
 double get_index(double p, double wavelength, double T)
 {
     /* Returns the index of refraction of pure water for a given density,
diff --git a/src/optics.h b/src/optics.h
index e2df659..5a6dec9 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -1,6 +1,7 @@
 #ifndef OPTICS_H
 #define OPTICS_H
 
+double get_absorption_length_snoman_d2o(double wavelength);
 double get_index(double p, double wavelength, double T);
 double get_index_snoman_h2o(double wavelength);
 double get_index_snoman_d2o(double wavelength);