add muon critical energy for D2O

1 files changed, 10 insertions, 2 deletions

diff --git a/src/muon.c b/src/muon.c
index 8db6237..54a1aa1 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -25,7 +25,15 @@ static gsl_spline *spline_dEdx;
 static gsl_interp_accel *acc_range;
 static gsl_spline *spline_range;
 
-static const double MUON_CRITICAL_ENERGY = 1.029e6;
+/* Muon critical energy in H2O and D2O. These values are used in computing the
+ * kinetic energy of the muon as a function of distance in get_T().
+ *
+ * These values come from
+ * http://pdgprod.lbl.gov/~deg/AtomicNuclearProperties/HTML/deuterium_oxide_liquid.html for D2O and
+ * http://pdg.lbl.gov/2018/AtomicNuclearProperties/HTML/water_liquid.html for H2O.
+ */
+static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6;
+static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3;
 
 static int init()
 {
@@ -184,7 +192,7 @@ double get_T(double T0, double x, double rho)
 
     /* This comes from Equation 33.42 in the PDG Passage of Particles Through
      * Matter article. */
-    b = log(1 + T0/MUON_CRITICAL_ENERGY)/range;
+    b = log(1 + T0/MUON_CRITICAL_ENERGY_D2O)/range;
     /* Now we compute the ionization energy loss from the known range and b. */
     a = b*T0/(exp(b*range)-1.0);