add some phase space factors for the integral over the earth

1 files changed, 16 insertions, 11 deletions

diff --git a/calculate_limits.c b/calculate_limits.c
index cbcb3a5..6dd37f0 100644
--- a/calculate_limits.c
+++ b/calculate_limits.c
@@ -2,11 +2,14 @@
 #include <gsl/gsl_integration.h>
 #include <math.h> /* For M_PI */
 
+/* Mass of dark matter particle (MeV). */
+double mass = 1000.0;
+
 /* Decay length of mediator V (in mm). */
 double decay_length = 1000e9;
 
-/* Mass of dark matter particle (MeV). */
-double mass = 1000.0;
+/* Cross section for dark matter interaction (in mm^2). */
+double dm_cross_section = 1e-30;
 
 /* Approximate dark matter density in MeV/mm^3. From Tom Caldwell's thesis. */
 double dm_density = 400e3;
@@ -44,7 +47,7 @@ double radius_fiducial = 5000;
 double x_sno[3];
 
 double epsabs = 1e-1;
-double epsrel = 1e-2;
+double epsrel = 1e-1;
 
 double deg2rad(double deg)
 {
@@ -227,7 +230,11 @@ double f4_earth(double phi_earth, void *params)
 
     gsl_integration_workspace_free(w);
 
-    return result;
+    /* For now we assume the event rate is constant throughout the earth, so we
+     * are implicitly assuming that the cross section is pretty small. */
+    double flux = dm_velocity*dm_density/mass;
+
+    return dm_cross_section*number_density*flux*result*data[0]*data[0]*sin(data[1]);
 }
 
 double f3_earth(double theta_earth, void *params)
@@ -281,18 +288,16 @@ double get_event_rate(double dm_mass, double gamma_length, double cs)
     F.function = &f2_earth;
     F.params = NULL;
 
-    /* For now we are storing the gamma decay length as a global variable. */
+    /* For now we just use global variables. */
+    mass = dm_mass;
     decay_length = gamma_length;
+    dm_cross_section = cs;
 
     gsl_integration_qags(&F, 0, radius_earth, epsabs, epsrel, 1000, w, &result, &error);
 
     gsl_integration_workspace_free(w);
 
-    /* For now we assume the event rate is constant throughout the earth, so we
-     * are implicitly assuming that the cross section is pretty small. */
-    double flux = dm_velocity*dm_density/dm_mass;
-
-    return cs*number_density*flux*result;
+    return result;
 }
 
 int main(int argc, char **argv)
@@ -307,7 +312,7 @@ int main(int argc, char **argv)
 
     /* Calculate the event rate for a standard DM candidate with a mass of 1
      * GeV, and a mediator decay length of 1 m. */
-    printf("event rate = %.18e Hz\n", get_event_rate(1000, 1000, 1e-30));
+    printf("event rate = %.18e Hz\n", get_event_rate(1000, 1000e9, 1e-30));
 
     return 0;
 }