diff options
Diffstat (limited to 'calculate_limits.c')
| -rw-r--r-- | calculate_limits.c | 36 |
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; } |