update theta0 calculation to use the radiation length in light water

Previously I was using the radiation length in light water but scaling it by
the density of heavy water, which isn't correct. Since the radiation length in
heavy and light water is almost identical, we just use the radiation length in
light water.

3 files changed, 11 insertions, 5 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 7da81bb..7c3d10f 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -530,7 +530,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
     beta0 = p0/E0;
 
     /* FIXME: is this formula valid for muons? */
-    theta0 = get_scattering_rms(range/2,p0,beta0,1.0,HEAVY_WATER_DENSITY)/sqrt(range/2);
+    theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2);
 
     params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
 
diff --git a/src/pdg.c b/src/pdg.c
index 780fbc8..a3eb2d5 100644
--- a/src/pdg.c
+++ b/src/pdg.c
@@ -1,7 +1,8 @@
 #include "pdg.h"
 #include "math.h"
+#include "sno.h"
 
-double get_scattering_rms(double x, double p, double beta, double z, double rho)
+double get_scattering_rms(double x, double p, double beta, double z)
 {
     /* Returns the RMS width of the scattering angle for a particle deflected
      * by many small-angle scatters after a distance `x`. `p` is the momentum
@@ -22,5 +23,5 @@ double get_scattering_rms(double x, double p, double beta, double z, double rho)
      * See Equation 33.15 in
      * http://pdg.lbl.gov/2018/reviews/rpp2018-rev-passage-particles-matter.pdf. */
     if (x == 0.0) return 0.0;
-    return (13.6/(beta*p))*z*sqrt(x*rho/RADIATION_LENGTH)*(1+0.038*log((x*z*z)/(RADIATION_LENGTH*beta*beta/rho)));
+    return (13.6/(beta*p))*z*sqrt(x*WATER_DENSITY/RADIATION_LENGTH)*(1+0.038*log((x*z*z)/(RADIATION_LENGTH*beta*beta/WATER_DENSITY)));
 }
diff --git a/src/pdg.h b/src/pdg.h
index 98046f5..df73380 100644
--- a/src/pdg.h
+++ b/src/pdg.h
@@ -2,13 +2,18 @@
 #define PDG_H
 
 #define SPEED_OF_LIGHT 29.9792458 /* cm/ns */
-/* From http://pdg.lbl.gov/2017/AtomicNuclearProperties/HTML/water_liquid.html */
+/* From http://pdg.lbl.gov/2017/AtomicNuclearProperties/HTML/water_liquid.html.
+ *
+ * Technically there is a different radiation length for D2O and H2O, and each
+ * should be divided by the density of the water. However, the results for
+ * water and heavy water are very close so instead of having to keep track of
+ * which medium we are in we just use the radiation length in light water. */
 #define RADIATION_LENGTH 36.08 /* g/cm^2 */
 #define ELECTRON_MASS 0.5109989461 /* MeV */
 #define MUON_MASS 105.6583745 /* MeV */
 #define PROTON_MASS 938.272081 /* MeV */
 #define FINE_STRUCTURE_CONSTANT 7.297352566417e-3
 
-double get_scattering_rms(double x, double p, double beta, double z, double rho);
+double get_scattering_rms(double x, double p, double beta, double z);
 
 #endif