fix delta ray charge calculation

Previously I was calculating the expected number of delta ray photons when
integrating over the shower path, but since the delta rays are produced along
the particle path and not further out like the shower photons, this wasn't
correct. The normalization of the probability distribution for the photons
produced along the path was also not handled correctly.

This commit adds a new function called integrate_path_delta_ray() to compute
the expected number of photons from delta rays hitting each PMT. Currently this
means that the likelihood function for muons will be significantly slower than
previously, but hopefully I can speed it up again in the future (for example by
skipping the shower calculation which is negligible for lower energy muons).

1 files changed, 10 insertions, 9 deletions

diff --git a/src/test.c b/src/test.c
index 2507546..936ccd5 100644
--- a/src/test.c
+++ b/src/test.c
@@ -1257,21 +1257,22 @@ int test_norm_cdf(char *err)
 
 static double gsl_time_pdf(double x, void *params)
 {
-    double mu_noise, mu_indirect, mu_direct, mu_shower, ts, ts_shower, tmean, tsigma, ts_sigma;
+    double mu_noise, mu_indirect, tmean;
+    double mu[2], ts[2], ts_sigma[2];
 
     double *pars = (double *) params;
 
     mu_noise = pars[0];
     mu_indirect = pars[1];
-    mu_direct = pars[2];
-    mu_shower = pars[3];
-    ts = pars[4];
-    ts_shower = pars[5];
+    mu[0] = pars[2];
+    mu[1] = pars[3];
+    ts[0] = pars[4];
+    ts[1] = pars[5];
     tmean = pars[6];
-    tsigma = pars[7];
-    ts_sigma = pars[8];
+    ts_sigma[0] = pars[7];
+    ts_sigma[1] = pars[8];
 
-    return time_pdf(x,mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma);
+    return time_pdf(x,mu_noise,mu_indirect,mu,2,ts,tmean,ts_sigma);
 }
 
 int test_time_pdf_norm(char *err)
@@ -1378,7 +1379,7 @@ int test_time_cdf(char *err)
             goto err;
         }
 
-        expected = time_cdf(t,params[0],params[1],&params[2],&params[3],1,&params[4],&params[5],params[6],params[7],&params[8]);
+        expected = time_cdf(t,params[0],params[1],&params[2],2,&params[4],params[6],&params[7]);
 
         if (!isclose(result, expected, 1e-2, 0)) {
             sprintf(err, "integral of time_pdf = %.5f, but expected %.5f", result, expected);