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, 43 insertions, 0 deletions

diff --git a/src/scattering.c b/src/scattering.c
index c4bdd57..68256dc 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -14,6 +14,8 @@
 #include <gsl/gsl_interp.h>
 #include <gsl/gsl_spline.h>
 
+static int initialized = 0;
+
 static double xlo = -1.0;
 static double xhi = 1.0;
 static size_t nx = 1000;
@@ -36,6 +38,9 @@ static size_t nx2 = 1000;
 static double *x2;
 static double *y2;
 
+/* Quantum efficiency weighted by the Cerenkov spectrum. */
+static double weighted_qe;
+
 static gsl_spline *spline2;
 static gsl_interp_accel *xacc2;
 
@@ -73,6 +78,30 @@ static double prob_scatter2(double wavelength, void *params)
     return 2*M_PI*FINE_STRUCTURE_CONSTANT*(1-(1/(beta*beta*index*index)))*qe/pow(wavelength,2)*1e7;
 }
 
+double get_weighted_quantum_efficiency(void)
+{
+    /* Returns the quantum efficiency weighted by the Cerenkov wavelength
+     * distribution, i.e. the probability that a photon randomly sampled from
+     * the Cerenkov wavelenght distribution between 200 and 800 nm is detected. */
+    if (!initialized) {
+        fprintf(stderr, "quantum efficiency hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return weighted_qe;
+}
+
+static double gsl_quantum_efficiency(double wavelength, void *params)
+{
+    /* Returns the quantum efficiency times the Cerenkov wavelength
+     * distribution. */
+    double qe;
+
+    qe = get_quantum_efficiency(wavelength);
+
+    return qe/pow(wavelength,2);
+}
+
 void init_interpolation(void)
 {
     size_t i, j;
@@ -149,6 +178,20 @@ void init_interpolation(void)
 
     gsl_spline_init(spline2, x2, y2, nx2);
 
+    F.function = &gsl_quantum_efficiency;
+    F.params = params;
+
+    status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+    if (status) {
+        fprintf(stderr, "error integrating quantum efficiency distribution: %s\n", gsl_strerror(status));
+        exit(1);
+    }
+
+    weighted_qe = result*800/3.0;
+
+    initialized = 1;
+
     gsl_integration_cquad_workspace_free(w);
 }