update get_event_rate to return an event rate assuming some simplified model for the dark matter in the galaxy

1 files changed, 32 insertions, 4 deletions

diff --git a/calculate_limits.c b/calculate_limits.c
index eefb6b8..cbcb3a5 100644
--- a/calculate_limits.c
+++ b/calculate_limits.c
@@ -8,6 +8,22 @@ double decay_length = 1000e9;
 /* Mass of dark matter particle (MeV). */
 double mass = 1000.0;
 
+/* Approximate dark matter density in MeV/mm^3. From Tom Caldwell's thesis. */
+double dm_density = 400e3;
+
+/* Approximate dark matter velocity in mm/s. The true distribution is expected
+ * to be a Maxwell Boltzmann distribution which is modulated annually by the
+ * earth's rotation around the sun, but we just assume a single constant
+ * velocity here. From Tom Caldwell's thesis page 26. */
+double dm_velocity = 244e6;
+
+/* Number density of scatterers in the Earth.
+ *
+ * FIXME: Currently just set to the number density of atoms in water. Need to
+ * update this for rock, and in fact this will change near the detector since
+ * there is water outside the AV. */
+double number_density = 30e18; /* In 1/mm^3 */
+
 /* From Google maps. Probably not very accurate, but should be good enough for
  * this calculation. */
 double latitude = 46.471857;
@@ -27,7 +43,7 @@ double radius_fiducial = 5000;
  * they are used in some functions. */
 double x_sno[3];
 
-double epsabs = 1e-2;
+double epsabs = 1e-1;
 double epsrel = 1e-2;
 
 double deg2rad(double deg)
@@ -252,7 +268,10 @@ double f2_earth(double r_earth, void *params)
     return result;
 }
 
-double get_event_rate()
+/* Returns the event rate in SNO for a self-destructing dark matter particle
+ * with a mass of dm_mass, a dark photon decay length of gamma_length, and a
+ * cross section of cs (in mm^2). */
+double get_event_rate(double dm_mass, double gamma_length, double cs)
 {
     double result, error;
     gsl_function F;
@@ -262,11 +281,18 @@ double get_event_rate()
     F.function = &f2_earth;
     F.params = NULL;
 
+    /* For now we are storing the gamma decay length as a global variable. */
+    decay_length = gamma_length;
+
     gsl_integration_qags(&F, 0, radius_earth, epsabs, epsrel, 1000, w, &result, &error);
 
     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/dm_mass;
+
+    return cs*number_density*flux*result;
 }
 
 int main(int argc, char **argv)
@@ -279,7 +305,9 @@ int main(int argc, char **argv)
 
     sphere2cartesian(radius_earth - sno_depth, sno_theta, sno_phi, x_sno, x_sno+1, x_sno+2);
 
-    printf("volume = %.18e\n", get_event_rate());
+    /* 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));
 
     return 0;
 }