1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
|
#include <stdio.h>
#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;
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 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[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
* 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);
printf("volume = %.2g\n", get_event_rate());
return 0;
}
|
