/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago> * * This program is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) * any later version. * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * You should have received a copy of the GNU General Public License along with * this program. If not, see <https://www.gnu.org/licenses/>. */ #include <stdio.h> #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; /* 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; /* 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; double longitude = -81.186755; /* Radius of the earth in mm. */ double radius_earth = 6.371e9; /* Depth of the SNO detector in mm. Don't be fooled by all the digits. I just * converted 6800 feet -> mm. */ double sno_depth = 2072640; /* Fiducial volume in mm. */ double radius_fiducial = 5000; /* Cartesian coordinates of SNO in earth frame. They need to be global since * they are used in some functions. */ double x_sno[3]; double epsabs = 1e-1; double epsrel = 1e-1; double deg2rad(double deg) { return deg*M_PI/180.0; } double rad2deg(double rad) { return rad*180.0/M_PI; } /* Convert spherical coordinates to cartesian coordinates. * * See https://en.wikipedia.org/wiki/Spherical_coordinate_system. */ void sphere2cartesian(double r, double theta, double phi, double *x, double *y, double *z) { *x = r*sin(theta)*cos(phi); *y = r*sin(theta)*sin(phi); *z = r*cos(theta); } /* Convert cartesian coordinates to spherical coordinates. * * See https://en.wikipedia.org/wiki/Spherical_coordinate_system. */ void cartesian2sphere(double x, double y, double z, double *r, double *theta, double *phi) { *r = sqrt(x*x + y*y + z*z); *theta = acos(z/(*r)); *phi = atan2(y,x); } void cross(double *a, double *b, double *c) { c[0] = a[1]*b[2] - a[2]*b[1]; c[1] = a[2]*b[0] - a[0]*b[2]; c[2] = a[0]*b[1] - a[1]*b[0]; } double dot(double *a, double *b) { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; } double norm(double *a) { return sqrt(dot(a,a)); } void normalize(double *a) { double n = norm(a); a[0] /= n; a[1] /= n; a[2] /= n; } /* Rotate a vector x around the vector dir by an angle theta. */ void rotate(double *result, double *x, double *dir, double theta) { double a = dot(dir,x); double b[3]; double sin_theta = sin(theta); double cos_theta = cos(theta); /* Make sure the direction vector is normalized. */ normalize(dir); cross(x,dir,b); result[0] = x[0]*cos_theta + dir[0]*a*(1-cos_theta) + b[0]*sin_theta; result[1] = x[1]*cos_theta + dir[1]*a*(1-cos_theta) + b[1]*sin_theta; result[2] = x[2]*cos_theta + dir[2]*a*(1-cos_theta) + b[2]*sin_theta; } /* Rotate a vector in earth centered coordinates to SNO coordinates (doesn't do * the translation). */ void rotate_earth_to_sno(double *x_earth, double *x_sno) { double dir[3]; double z[3] = {0,0,1}; cross(x_sno, z, dir); /* Normalize. */ normalize(dir); double theta = acos(dot(x_sno,z)/norm(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[5] = phi; double x[3]; double r[3]; double distance; /* 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. */ double f2(double theta, void *params) { double result, error; gsl_function F; double *data = (double *) params; data[4] = theta; gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); F.function = &f3; F.params = params; gsl_integration_qags(&F, 0, 2*M_PI, epsabs, epsrel, 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[3] = r; gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); F.function = &f2; F.params = params; gsl_integration_qags(&F, 0, M_PI, epsabs, epsrel, 1000, w, &result, &error); gsl_integration_workspace_free(w); return result; } double f4_earth(double phi_earth, void *params) { double result, error; 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); F.function = &f1; F.params = params; gsl_integration_qags(&F, 0, radius_fiducial, 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/mass; return dm_cross_section*number_density*flux*result*data[0]*data[0]*sin(data[1]); } 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, epsabs, epsrel, 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 = &f3_earth; F.params = (void *) data; gsl_integration_qags(&F, 0, M_PI, epsabs, epsrel, 1000, w, &result, &error); gsl_integration_workspace_free(w); return result; } /* 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; gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000); F.function = &f2_earth; F.params = NULL; /* 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); return result; } int main(int argc, char **argv) { /* Spherical angles for the SNO detector in the earth frame which has z * along the north and south poles and the x axis passing through Greenwich. * Should double check this. */ double sno_theta = deg2rad(latitude + 90.0); double sno_phi = deg2rad(longitude); sphere2cartesian(radius_earth - sno_depth, sno_theta, sno_phi, x_sno, x_sno+1, x_sno+2); /* 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, 1000e9, 1e-30)); return 0; }