update likelihood function to speed it up

This commit makes a few small changes to try and reduce the number of divisions
and multiplications done in get_expected_charge() to speed up the likelihood
function.

1 files changed, 9 insertions, 9 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 17cde5d..950b646 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -547,18 +547,18 @@ 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, z, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o;
-
-    z = 1.0;
+    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;
 
     R = NORM(pos);
 
     n = (R <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
 
+    cos_theta_cerenkov = 1/(beta*n);
+
     *q = 0.0;
     *t = 0.0;
     *reflected = 0.0;
-    if (beta < 1/n) return;
+    if (cos_theta_cerenkov > 1) return;
 
     SUB(pmt_dir,pmts[pmt].pos,pos);
 
@@ -569,17 +569,17 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
     cos_theta = DOT(dir,pmt_dir);
     sin_theta = fast_sqrt(1-cos_theta*cos_theta);
 
-    cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
+    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));
 
-    if (fabs(cos_theta-1.0/(n*beta))/(sin_theta*theta0) > 5) return;
+    if (fabs(cos_theta-cos_theta_cerenkov)/(sin_theta*theta0) > 5) return;
 
     omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
 
-    f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
-    f = get_weighted_pmt_response(-cos_theta_pmt);
+    f_reflec = get_weighted_pmt_reflectivity(cos_theta_pmt);
+    f = get_weighted_pmt_response(cos_theta_pmt);
 
     get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
@@ -603,7 +603,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
      * pretty low, this is expected to be a very small effect. */
     prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
-    charge = omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, sin_theta, theta0);
+    charge = omega*(1-cos_theta_cerenkov*cos_theta_cerenkov)*get_probability(beta, cos_theta, sin_theta, theta0);
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;