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, 79 insertions, 1 deletions

diff --git a/src/solid_angle.c b/src/solid_angle.c
index 7201533..354547e 100644
--- a/src/solid_angle.c
+++ b/src/solid_angle.c
@@ -4,6 +4,7 @@
 #include <gsl/gsl_interp2d.h>
 #include <gsl/gsl_spline2d.h>
 #include "vector.h"
+#include "misc.h"
 
 static double hd[11] = {0.1,0.125,0.150,0.175,0.2,0.3,0.4,0.5,0.6,0.8,1.0};
 static double Rd[13] = {0.1,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0};
@@ -40,7 +41,8 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r)
      * regimes, the approximation is not very good and so it is modified by a
      * correction factor which comes from a lookup table.
      *
-     * This formula comes from "The Solid Angle Subtended at a Point by a * Circular Disk." Gardener et al. 1969. */
+     * This formula comes from "The Solid Angle Subtended at a Point by a
+     * Circular Disk." Gardener et al. 1969. */
     double dir[3];
     double h, r0, D, a, u1, u2, F;
     static gsl_spline2d *spline;
@@ -78,6 +80,82 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r)
     }
 }
 
+double get_solid_angle2(double L2, double r2)
+{
+    double k = 2*sqrt(r2/(pow(L2,2)+pow(1+r2,2)));
+    double a2 = 4*r2/pow(1+r2,2);
+    double L_Rmax = L2/sqrt(pow(L2,2)+pow(1+r2,2));
+    double r0_r = (r2-1)/(r2+1);
+    double K, P;
+
+    if (fabs(r2-1) < 1e-5) {
+        /* If r0 is very close to r, the GSL routines below will occasionally
+         * throw a domain error. */
+        K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+        return M_PI - 2*L_Rmax*K;
+    } else if (r2 <= 1) {
+        K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+        P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE);
+        return 2*M_PI - 2*L_Rmax*K + (2*L_Rmax)*r0_r*P;
+    } else {
+        K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+        P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE);
+        return -2*L_Rmax*K + (2*L_Rmax)*r0_r*P;
+    }
+}
+
+double get_solid_angle_lookup(double L2, double r2)
+{
+    size_t i, j;
+    static int initialized = 0;
+    static double xs[1000];
+    static double ys[1000];
+    static double zs[1000*1000];
+    static double xlo = 0.0;
+    static double xhi = 1000;
+    static double ylo = 0.0;
+    static double yhi = 1000;
+
+    if (!initialized) {
+        for (i = 0; i < LEN(xs); i++) {
+            xs[i] = xlo + (xhi-xlo)*i/(LEN(xs)-1);
+        }
+        for (i = 0; i < LEN(ys); i++) {
+            ys[i] = ylo + (yhi-ylo)*i/(LEN(ys)-1);
+        }
+        for (i = 0; i < LEN(xs); i++) {
+            for (j = 0; j < LEN(ys); j++) {
+                zs[i*LEN(ys) + j] = get_solid_angle2(xs[i],ys[j]);
+            }
+        }
+        initialized = 1;
+    }
+
+    if (L2 < xlo || L2 > xhi || r2 < ylo || r2 > yhi) {
+        /* If the arguments are out of bounds of the lookup table, just call
+         * get_solid_angle2(). */
+        return get_solid_angle2(L2,r2);
+    }
+
+    return interp2d(L2,r2,xs,ys,zs,LEN(xs),LEN(ys));
+}
+
+double get_solid_angle_fast(double *pos, double *pmt, double *n, double r)
+{
+    double dir[3];
+    double L, r0, R;
+
+    dir[0] = pos[0] - pmt[0];
+    dir[1] = pos[1] - pmt[1];
+    dir[2] = pos[2] - pmt[2];
+
+    L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+    R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+    r0 = sqrt(R*R - L*L);
+
+    return get_solid_angle_lookup(L/r,r0/r);
+}
+
 double get_solid_angle(double *pos, double *pmt, double *n, double r)
 {
     /* Returns the solid angle subtended by a circular disk of radius r at a