fix how the RMS scattering angle is calculated

The RMS scattering angle calculation comes from Equation 33.15 in the PDG
article on the passage of particles through matter. It's not entirely obvious
if this equation is correct for a long track. It seems like it should be
integrated along the track to add up the contributions at different energies,
but it's not obvious how to do that with the log term.

In any case, the way I was previously calculating it (by using the momentum and
velocity at each point along the track) was definitely wrong.

I will try this out and perhaps try to integrate it later.

5 files changed, 13 insertions, 8 deletions

diff --git a/src/fit.c b/src/fit.c
index b6e441a..1f90f12 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -134,7 +134,7 @@ int main(int argc, char **argv)
                     /* New event, so we need to fit the old event. */
                     iter = 0;
 
-                    gsl_vector_set(x,0,1000.0);
+                    gsl_vector_set(x,0,500.0);
                     gsl_vector_set(x,1,100.0);
                     gsl_vector_set(x,2,100.0);
                     gsl_vector_set(x,3,100.0);
diff --git a/src/likelihood.c b/src/likelihood.c
index 98a4ad7..4e24a79 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -103,7 +103,7 @@ static double gsl_muon_time(double x, void *params)
 
     t = x/SPEED_OF_LIGHT + distance*n/SPEED_OF_LIGHT;
 
-    return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+    return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
 }
 
 static double gsl_muon_charge(double x, void *params)
@@ -124,7 +124,7 @@ static double gsl_muon_charge(double x, void *params)
     pos[1] = pos0[1] + dir[1]*x;
     pos[2] = pos0[2] + dir[2]*x;
 
-    return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+    return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
 }
 
 double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
diff --git a/src/likelihood.h b/src/likelihood.h
index ebfeb97..43b0823 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -4,7 +4,7 @@
 #include "event.h"
 
 #define MAX_PE 100
-#define CHARGE_FRACTION 150000.0
+#define CHARGE_FRACTION 60000.0
 #define DARK_RATE 500.0
 /* Single PE transit time spread (ns). */
 #define PMT_TTS 1.5
diff --git a/src/muon.c b/src/muon.c
index 0279984..9a539de 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -216,9 +216,9 @@ double get_dEdx(double T, double rho)
     return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
 }
 
-double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
+double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
 {
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R;
+    double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R, E0, p0, beta0;
 
     z = 1.0;
 
@@ -253,8 +253,13 @@ double get_expected_charge(double x, double T, double *pos, double *dir, double
 
     if (beta < 1/n) return 0;
 
+    /* Calculate total energy */
+    E0 = T0 + MUON_MASS;
+    p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+    beta0 = p0/E0;
+
     /* FIXME: is this formula valid for muons? */
-    theta0 = get_scattering_rms(x,p,beta,z,rho);
+    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);
diff --git a/src/muon.h b/src/muon.h
index a592591..e95cb3c 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -6,6 +6,6 @@
 double get_range(double T, double rho);
 double get_T(double T0, double x, double rho);
 double get_dEdx(double T, double rho);
-double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
+double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
 
 #endif