fix how multiple Coulomb scattering is treated

Previously I had been assuming that a particle undergoing many small angle
Coulomb scatters had a track direction whose polar angle was a Gaussian.
However, this was just due to a misunderstanding of the PDG section "Multiple
scattering through small angles" in the "Passage of particles through matter"
article. In fact, what is described by a Gaussian is the polar angle projected
onto a plane. Therefore the distribution of the polar angle is actually:

(1/(sqrt(2*pi)*theta0**2))*theta*exp(-theta**2/(2*theta0))

This commit updates the code in scattering.c to correctly calculate the
probability that a photon is emitted at a particular angle. I also updated
test-likelihood.c to simulate a track correctly.

1 files changed, 3 insertions, 6 deletions

diff --git a/src/scattering.c b/src/scattering.c
index a1a63df..66e8398 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -40,12 +40,9 @@ static double prob_scatter(double wavelength, void *params)
 
     index = get_index_snoman_d2o(wavelength);
 
-    delta = (1.0/index - beta_cos_theta)/(2*beta_sin_theta_theta0);
+    delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0;
 
-    /* FIXME: ignore GSL error for underflow here. */
-    if (delta*delta > 500) return 0.0;
-    else if (delta*delta == 0.0) return INFINITY;
-    return qe*exp(-delta*delta)*gsl_sf_bessel_K0(delta*delta)/pow(wavelength,2)*1e7/(4*M_PI*M_PI);
+    return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI);
 }
 
 void init_interpolation(void)
@@ -106,7 +103,7 @@ double get_probability(double beta, double cos_theta, double theta0)
      * we are going to square it everywhere. */
     sin_theta = fabs(sin(acos(cos_theta)));
 
-    return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/sin_theta;
+    return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta);
 }
 
 void free_interpolation(void)