add the integration over the gamma production point

1 files changed, 85 insertions, 7 deletions

diff --git a/calculate_limits.c b/calculate_limits.c
index c4919d0..07fe052 100644
--- a/calculate_limits.c
+++ b/calculate_limits.c
@@ -124,9 +124,27 @@ double f3(double phi, void *params)
     double result, error;
     gsl_function F;
     double *data = (double *) params;
-    data[2] = phi;
+    data[5] = phi;
+    double x[3];
+    double r[3];
+    double distance;
 
-    return 1;
+    /* Compute cartesian position in local SNO coordinates. */
+    sphere2cartesian(data[3], data[4], data[5], &x[0], &x[1], &x[2]);
+
+    /* Cartesian coordinates of gamma production offset in earth centered
+     * coordinates .*/
+    double *gamma_offset = data+6;
+
+    /* Vector distance between integration in local coordinates and gamma
+     * production point .*/
+    r[0] = x_sno[0] + x[0] - gamma_offset[0];
+    r[1] = x_sno[1] + x[1] - gamma_offset[1];
+    r[2] = x_sno[2] + x[2] - gamma_offset[2];
+
+    distance = norm(r);
+
+    return exp(-distance/decay_length)/(4*M_PI*distance*distance*decay_length)*data[3]*data[3]*sin(data[4]);
 }
 
 /* Integral over theta. */
@@ -135,7 +153,7 @@ double f2(double theta, void *params)
     double result, error;
     gsl_function F;
     double *data = (double *) params;
-    data[1] = theta;
+    data[4] = theta;
 
     gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
 
@@ -155,7 +173,7 @@ double f1(double r, void *params)
     double result, error;
     gsl_function F;
     double *data = (double *) params;
-    data[0] = r;
+    data[3] = r;
 
     gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
 
@@ -169,12 +187,17 @@ double f1(double r, void *params)
     return result;
 }
 
-/* Integral over r. */
-double get_event_rate()
+double f4_earth(double phi_earth, void *params)
 {
     double result, error;
-    double params[3];
     gsl_function F;
+    double *data = (double *) params;
+    data[2] = phi_earth;
+    double gamma_offset[3];
+
+    /* Compute the cartesian coordinates of the gamma production point in the
+     * earth centered coordinates. */
+    sphere2cartesian(data[0], data[1], data[2], &data[6], &data[7], &data[8]);
 
     gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
 
@@ -188,6 +211,61 @@ double get_event_rate()
     return result;
 }
 
+double f3_earth(double theta_earth, void *params)
+{
+    double result, error;
+    gsl_function F;
+    double *data = (double *) params;
+    data[1] = theta_earth;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f4_earth;
+    F.params = params;
+
+    gsl_integration_qags(&F, 0, 2*M_PI, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
+double f2_earth(double r_earth, void *params)
+{
+    double result, error;
+    double data[9];
+    gsl_function F;
+    data[0] = r_earth;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f2;
+    F.params = (void *) data;
+
+    gsl_integration_qags(&F, 0, M_PI, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
+double get_event_rate()
+{
+    double result, error;
+    gsl_function F;
+
+    gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
+
+    F.function = &f2_earth;
+    F.params = NULL;
+
+    gsl_integration_qags(&F, 0, radius_earth, 0, 1e-7, 1000, w, &result, &error);
+
+    gsl_integration_workspace_free(w);
+
+    return result;
+}
+
 int main(int argc, char **argv)
 {
     /* Spherical angles for the SNO detector in the earth frame which has z