switch back to calling cquad just once to speed things up

I found when simulating high energy muons that the expected charge for some
PMTs which should be getting hit was zero. The reason for this is that the
integrand was very sharply peaked at the Cerenkov angle which makes it
difficult to integrate for numerical integration routines like cquad. To solve
this I split up the integral at the point when the track was at the Cerenkov
angle from the PMT to make sure that cquad didn't miss the peak. However,
calling cquad twice takes a lot of time so it's not necessarily good to do this
for all fits. Also, it's not obvious if it is necessary any more now that the
angular distribution calculation was fixed.

I think the real reason that cquad was missing those integrals was that for a
high energy muon the range is going to be very large (approximately 40 meters
for a 10 GeV muon). In this case, I should really only integrate up to the edge
of the cavity or PSUP and hopefully cquad picks enough points in there to get a
non zero value.

I also added a check to only compute tmean when at least one PMT has a valid
time. This prevents a divide by zero which causes the likelihood function to
return nan.

1 files changed, 7 insertions, 22 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index f06fade..c90914b 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -201,33 +201,17 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
          * Cerenkov angle. */
         x = R*sin(theta_cerenkov-theta)/sin(theta_cerenkov);
 
-        if (x > 0 && x < range) {
-            /* Split up the integral at the point where the PMT is at the
-             * Cerenkov angle. */
-            F.function = &gsl_muon_charge;
-            gsl_integration_cquad(&F, 0, x, 0, 1e-2, w, &result, &error, &nevals);
-            mu_direct[i] = result;
-            gsl_integration_cquad(&F, x, range, 0, 1e-2, w, &result, &error, &nevals);
-            mu_direct[i] += result;
+        F.function = &gsl_muon_charge;
+        gsl_integration_cquad(&F, 0, range, 0, 1e-2, w, &result, &error, &nevals);
+        mu_direct[i] = result;
 
-            F.function = &gsl_muon_time;
-            gsl_integration_cquad(&F, 0, x, 0, 1e-2, w, &result, &error, &nevals);
-            ts[i] = result;
-            gsl_integration_cquad(&F, x, range, 0, 1e-2, w, &result, &error, &nevals);
-            ts[i] += result;
-        } else {
-            F.function = &gsl_muon_charge;
-            gsl_integration_cquad(&F, 0, range, 0, 1e-2, w, &result, &error, &nevals);
-            mu_direct[i] = result;
+        total_charge += mu_direct[i];
 
+        if (mu_direct[i] > 0.001) {
             F.function = &gsl_muon_time;
             gsl_integration_cquad(&F, 0, range, 0, 1e-2, w, &result, &error, &nevals);
             ts[i] = result;
-        }
 
-        total_charge += mu_direct[i];
-
-        if (mu_direct[i] > 0.001) {
             ts[i] /= mu_direct[i];
             ts[i] += t0;
             tmean += ts[i];
@@ -239,7 +223,8 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
 
     path_free(params.p);
 
-    tmean /= npmt;
+    if (npmt)
+        tmean /= npmt;
 
     gsl_integration_cquad_workspace_free(w);