update likelihood function to fit electrons!

To characterize the angular distribution of photons from an electromagnetic
shower I came up with the following functional form:

    f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)

and fit this to data simulated using RAT-PAC at several different energies. I
then fit the alpha and beta coefficients as a function of energy to the
functional form:

    alpha = c0 + c1/log(c2*T0 + c3)
    beta  = c0 + c1/log(c2*T0 + c3).

where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3
are parameters which I fit.

The longitudinal distribution of the photons generated from an electromagnetic
shower is described by a gamma distribution:

    f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a).

This parameterization comes from the PDG "Passage of particles through matter"
section 32.5. I also fit the data from my RAT-PAC simulation, but currently I
am not using it, and instead using a simpler form to calculate the coefficients
from the PDG (although I estimated the b parameter from the RAT-PAC data).

I also sped up the calculation of the solid angle by making a lookup table
since it was taking a significant fraction of the time to compute the
likelihood function.

1 files changed, 58 insertions, 27 deletions

diff --git a/src/fit.c b/src/fit.c
index eba125e..82586f9 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -22,6 +22,7 @@
 #include <signal.h> /* for signal() */
 #include "release.h"
 #include "id_particles.h"
+#include "misc.h"
 
 char *GitSHA1(void);
 char *GitDirty(void);
@@ -40,8 +41,14 @@ typedef struct fitParams {
     int n;
     int fast;
     int print;
+    int id;
 } fitParams;
 
+int particles[] = {
+    IDP_E_MINUS,
+    IDP_MU_MINUS,
+};
+
 /* In order to start the fitter close to the minimum, we first do a series of
  * "quick" minimizations starting at the following points. We keep track of the
  * parameters with the best likelihood value and then start the "real"
@@ -4996,7 +5003,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
     z2[0] = x[8];
 
     gettimeofday(&tv_start, NULL);
-    fval = nll_muon(fpars->ev, IDP_MU_MINUS, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast);
+    fval = nll_muon(fpars->ev, fpars->id, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast);
     gettimeofday(&tv_stop, NULL);
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -5099,11 +5106,11 @@ void guess_direction(event *ev, double *pos, double *theta, double *phi)
     *phi = atan2(dir[1],dir[0]);
 }
 
-int fit_event(event *ev, double *xopt, double *fmin)
+int fit_event(event *ev, double *xopt, double *fmin, int id)
 {
     size_t i;
     fitParams fpars;
-    double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin;
+    double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin, mass;
     struct timeval tv_start, tv_stop;
 
     opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
@@ -5124,12 +5131,29 @@ int fit_event(event *ev, double *xopt, double *fmin)
 
     n = get_index_snoman_d2o(400.0);
 
+    switch (id) {
+    case IDP_E_MINUS:
+    case IDP_E_PLUS:
+        mass = ELECTRON_MASS;
+        break;
+    case IDP_MU_MINUS:
+    case IDP_MU_PLUS:
+        mass = MUON_MASS;
+        break;
+    case IDP_PROTON:
+        mass = PROTON_MASS;
+        break;
+    default:
+        fprintf(stderr, "unknown particle id %i\n", id);
+        exit(1);
+    }
+
     /* Calculate the Cerenkov threshold for a muon. */
-    Tmin = MUON_MASS/sqrt(1.0-1.0/(n*n)) - MUON_MASS;
+    Tmin = mass/sqrt(1.0-1.0/(n*n)) - mass;
 
     /* If our guess is below the Cerenkov threshold, start at the Cerenkov
      * threshold. */
-    double Tmin2 = sqrt(MUON_MASS*MUON_MASS/(1-pow(0.9,2)))-MUON_MASS;
+    double Tmin2 = sqrt(mass*mass/(1-pow(0.9,2)))-mass;
     if (T0 < Tmin2) T0 = Tmin2;
 
     x0[0] = T0;
@@ -5178,6 +5202,7 @@ int fit_event(event *ev, double *xopt, double *fmin)
     nlopt_set_initial_step(opt, ss);
 
     fpars.ev = ev;
+    fpars.id = id;
 
     iter = 0;
 
@@ -5538,32 +5563,38 @@ int main(int argc, char **argv)
 
             rv = get_event(f,&ev,&b);
 
-            gettimeofday(&tv_start, NULL);
-            if (fit_event(&ev,xopt,&fmin) == NLOPT_FORCED_STOP) {
-                printf("ctrl-c caught. quitting...\n");
-                goto end;
-            }
-            gettimeofday(&tv_stop, NULL);
-
-            long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
-
             if (fout) {
                 if (first_ev) fprintf(fout, "    ev:\n");
                 fprintf(fout, "      - gtid: %i\n", ev.gtid);
                 fprintf(fout, "        fit:\n");
-                fprintf(fout, "          -\n");
-                fprintf(fout, "            energy: %.2f\n", xopt[0]);
-                fprintf(fout, "            posx: %.2f\n", xopt[1]);
-                fprintf(fout, "            posy: %.2f\n", xopt[2]);
-                fprintf(fout, "            posz: %.2f\n", xopt[3]);
-                fprintf(fout, "            theta: %.4f\n", xopt[4]);
-                fprintf(fout, "            phi: %.4f\n", xopt[5]);
-                fprintf(fout, "            t0: %.2f\n", xopt[6]);
-                fprintf(fout, "            z1: %.2f\n", xopt[7]);
-                fprintf(fout, "            z2: %.2f\n", xopt[8]);
-                fprintf(fout, "            fmin: %.2f\n", fmin);
-                fprintf(fout, "            time: %lld\n", elapsed);
-                fflush(fout);
+            }
+
+            for (i = 0; i < LEN(particles); i++) {
+                gettimeofday(&tv_start, NULL);
+                if (fit_event(&ev,xopt,&fmin,particles[i]) == NLOPT_FORCED_STOP) {
+                    printf("ctrl-c caught. quitting...\n");
+                    goto end;
+                }
+                gettimeofday(&tv_stop, NULL);
+
+                long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+
+                if (fout) {
+                    fprintf(fout, "          -\n");
+                    fprintf(fout, "            id: %i\n", particles[i]);
+                    fprintf(fout, "            energy: %.2f\n", xopt[0]);
+                    fprintf(fout, "            posx: %.2f\n", xopt[1]);
+                    fprintf(fout, "            posy: %.2f\n", xopt[2]);
+                    fprintf(fout, "            posz: %.2f\n", xopt[3]);
+                    fprintf(fout, "            theta: %.4f\n", xopt[4]);
+                    fprintf(fout, "            phi: %.4f\n", xopt[5]);
+                    fprintf(fout, "            t0: %.2f\n", xopt[6]);
+                    fprintf(fout, "            z1: %.2f\n", xopt[7]);
+                    fprintf(fout, "            z2: %.2f\n", xopt[8]);
+                    fprintf(fout, "            fmin: %.2f\n", fmin);
+                    fprintf(fout, "            time: %lld\n", elapsed);
+                    fflush(fout);
+                }
             }
 
             first_ev = 0;