fix a bug in theta0 min calculation

Previously I was accidentally passing the absolute position of the particle
instead of the distance to the PMT to get_theta0_min().

1 files changed, 4 insertions, 2 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 344f320..9037cb9 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -550,7 +550,7 @@ double get_theta0_min(double R, double r, double sin_theta_pmt)
 
 static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
 {
-    double pmt_dir[3], cos_theta, sin_theta, n, omega, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov;
+    double pmt_dir[3], cos_theta, sin_theta, n, omega, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov, distance_to_pmt;
 
     R = NORM(pos);
 
@@ -565,6 +565,8 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
 
     SUB(pmt_dir,pmts[pmt].pos,pos);
 
+    distance_to_pmt = NORM(pmt_dir);
+
     normalize(pmt_dir);
 
     /* Calculate the cosine of the angle between the track direction and the
@@ -575,7 +577,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
     cos_theta_pmt = -DOT(pmt_dir,pmts[pmt].normal);
     sin_theta_pmt = fast_sqrt(1-cos_theta_pmt*cos_theta_pmt);
 
-    theta0 = fmax(theta0,get_theta0_min(R,PMT_RADIUS,sin_theta_pmt));
+    theta0 = fmax(theta0,get_theta0_min(distance_to_pmt,PMT_RADIUS,sin_theta_pmt));
 
     if (fabs(cos_theta-cos_theta_cerenkov)/(sin_theta*theta0) > 5) return;