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, 15 insertions, 15 deletions

diff --git a/src/test-time-pdf.c b/src/test-time-pdf.c
index f754b49..6be404f 100644
--- a/src/test-time-pdf.c
+++ b/src/test-time-pdf.c
@@ -41,18 +41,18 @@ int main(int argc, char **argv)
 {
     size_t i, n;
     double *x;
-    double mu_noise, mu_indirect, mu_direct, mu_shower, ts, ts_shower, tmean, tsigma, ts_sigma;
+    double mu[2], ts[2], ts_sigma[2], mu_noise, mu_indirect, tmean;
 
     n = 1000;
     mu_noise = 0.1;
     mu_indirect = 0.5;
-    mu_direct = 1.0;
-    mu_shower = 1.0;
-    ts = 100.0;
-    ts_shower = 120.0;
+    mu[0] = 1.0;
+    mu[1] = 1.0;
+    ts[0] = 100.0;
+    ts[1] = 120.0;
     tmean = 100.0;
-    tsigma = PMT_TTS;
-    ts_sigma = 10.0;
+    ts_sigma[0] = PMT_TTS;
+    ts_sigma[1] = 10.0;
 
     for (i = 1; i < argc; i++) {
         if (!strncmp(argv[i], "--", 2)) {
@@ -63,25 +63,25 @@ int main(int argc, char **argv)
                 mu_indirect = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"mu-direct")) {
-                mu_direct = strtod(argv[++i],NULL);
+                mu[0] = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"mu-shower")) {
-                mu_shower = strtod(argv[++i],NULL);
+                mu[1] = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"ts")) {
-                ts = strtod(argv[++i],NULL);
+                ts[0] = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"ts-shower")) {
-                ts_shower = strtod(argv[++i],NULL);
+                ts[1] = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"t-mean")) {
                 tmean = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"t-sigma")) {
-                tsigma = strtod(argv[++i],NULL);
+                ts_sigma[0] = strtod(argv[++i],NULL);
                 continue;
             } else if (!strcmp(argv[i]+2,"ts-sigma")) {
-                ts_sigma = strtod(argv[++i],NULL);
+                ts_sigma[1] = strtod(argv[++i],NULL);
                 continue;
             }
         } else if (argv[i][0] == '-') {
@@ -112,12 +112,12 @@ int main(int argc, char **argv)
     }
 
     for (i = 0; i < n; i++) {
-        fprintf(pipe, "%.10f %.10f\n", x[i], time_pdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma));
+        fprintf(pipe, "%.10f %.10f\n", x[i], time_pdf(x[i],mu_noise,mu_indirect,mu,2,ts,tmean,ts_sigma));
     }
     fprintf(pipe, "\n\n");
 
     for (i = 0; i < n; i++) {
-        fprintf(pipe, "%.10f %.10f\n", x[i], time_cdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma));
+        fprintf(pipe, "%.10f %.10f\n", x[i], time_cdf(x[i],mu_noise,mu_indirect,mu,2,ts,tmean,ts_sigma));
     }
     fprintf(pipe, "\n\n");