diff options
| -rw-r--r-- | Makefile | 2 | ||||
| -rw-r--r-- | calculate_limits.c | 78 |
2 files changed, 79 insertions, 1 deletions
@@ -1,3 +1,3 @@ -LDLIBS=-lm +LDLIBS=-lm -lgsl -lgslcblas calculate_limits: calculate_limits.c diff --git a/calculate_limits.c b/calculate_limits.c index 47ed1c9..c4919d0 100644 --- a/calculate_limits.c +++ b/calculate_limits.c @@ -2,6 +2,12 @@ #include <gsl/gsl_integration.h> #include <math.h> /* For M_PI */ +/* Decay length of mediator V (in mm). */ +double decay_length = 1000.0; + +/* Mass of dark matter particle (MeV). */ +double mass = 1000.0; + /* From Google maps. Probably not very accurate, but should be good enough for * this calculation. */ double latitude = 46.471857; @@ -112,6 +118,76 @@ void rotate_earth_to_sno(double *x_earth, double *x_sno) rotate(x_sno, x_earth, dir, theta); } +/* Integral over phi. */ +double f3(double phi, void *params) +{ + double result, error; + gsl_function F; + double *data = (double *) params; + data[2] = phi; + + return 1; +} + +/* Integral over theta. */ +double f2(double theta, void *params) +{ + double result, error; + gsl_function F; + double *data = (double *) params; + data[1] = theta; + + gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); + + F.function = &f3; + 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; +} + +/* Integral over r. */ +double f1(double r, void *params) +{ + double result, error; + gsl_function F; + double *data = (double *) params; + data[0] = r; + + gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); + + F.function = &f2; + F.params = params; + + gsl_integration_qags(&F, 0, M_PI, 0, 1e-7, 1000, w, &result, &error); + + gsl_integration_workspace_free(w); + + return result; +} + +/* Integral over r. */ +double get_event_rate() +{ + double result, error; + double params[3]; + gsl_function F; + + gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); + + F.function = &f2; + F.params = params; + + gsl_integration_qags(&F, 0, radius_fiducial, 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 @@ -122,5 +198,7 @@ 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 = %.2g\n", get_event_rate()); + return 0; } |