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.

3 files changed, 26 insertions, 22 deletions

diff --git a/src/muon.c b/src/muon.c
index 9a539de..7012214 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -262,5 +262,5 @@ double get_expected_charge(double x, double T, double T0, double *pos, double *d
     theta0 = get_scattering_rms(x,p0,beta0,z,rho);
 
     /* FIXME: add angular response and scattering/absorption. */
-    return 2*omega*2*M_PI*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0)/(sqrt(2*M_PI)*theta0);
+    return omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
 }
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)
diff --git a/src/test-likelihood.c b/src/test-likelihood.c
index f0266b9..45bbfe2 100644
--- a/src/test-likelihood.c
+++ b/src/test-likelihood.c
@@ -21,7 +21,7 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t
      * distribution is simulated as a gaussian distribution with standard
      * deviation `theta0`. */
     int i;
-    double theta, phi, wavelength, u, qe, index, cerenkov_angle, dir[3], n[3], dest[3], E, p, beta, cos_theta;
+    double theta, phi, wavelength, u, qe, index, cerenkov_angle, dir[3], n[3], dest[3], E, p, beta, cos_theta, thetax, thetay;
 
     i = 0;
     while (i < N) {
@@ -35,7 +35,7 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t
         /* Check to see if the photon was detected. */
         if (genrand_real2() > qe) continue;
 
-        index = get_index(HEAVY_WATER_DENSITY, wavelength, 10.0);
+        index = get_index_snoman_d2o(wavelength);
 
         /* Calculate total energy */
         E = T + MUON_MASS;
@@ -45,21 +45,28 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t
         cerenkov_angle = acos(1/(index*beta));
 
         /* Assuming the muon track is dominated by small angle scattering, the
-         * angular distribution will be a Gaussian centered around 0 with a
-         * standard deviation of `theta0`. Here, we draw a random angle from
-         * this distribution. */
-        theta = randn()*theta0;
-
-        n[0] = sin(theta);
-        n[1] = 0;
+         * angular distribution looks like the product of two uncorrelated
+         * Gaussian distributions with a standard deviation of `theta0` in the
+         * plane perpendicular to the track direction. Here, we draw two random
+         * angles and then compute the polar and azimuthal angle for the track
+         * direction. */
+        thetax = randn()*theta0;
+        thetay = randn()*theta0;
+
+        theta = sqrt(thetax*thetax + thetay*thetay);
+        phi = atan2(thetay,thetax);
+
+        n[0] = sin(theta)*cos(phi);
+        n[1] = sin(theta)*sin(phi);
         n[2] = cos(theta);
 
-        /* To compute the direction of the photon, we start with a vector in
-         * the x-z plane which is offset from the track direction by the
-         * Cerenkov angle and then rotate it around the track direction by a
-         * random angle `phi`. */
-        dir[0] = sin(cerenkov_angle + theta);
-        dir[1] = 0;
+        /* To compute the direction of the photon, we start with a vector which
+         * has the same azimuthal angle as the track direction but is offset
+         * from the track direction in the polar angle by the Cerenkov angle
+         * and then rotate it around the track direction by a random angle
+         * `phi`. */
+        dir[0] = sin(cerenkov_angle + theta)*cos(phi);
+        dir[1] = sin(cerenkov_angle + theta)*sin(phi);
         dir[2] = cos(cerenkov_angle + theta);
 
         phi = genrand_real2()*2*M_PI;