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
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
|
#include "scattering.h"
#include "quantum_efficiency.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include "optics.h"
#include "sno.h"
#include <stddef.h> /* for size_t */
#include <stdlib.h> /* for exit() */
#include <stdio.h> /* for fprintf() */
#include <gsl/gsl_errno.h> /* for gsl_strerror() */
#include "pdg.h"
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
static int initialized = 0;
static double xlo = -1.0;
static double xhi = 1.0;
static size_t nx = 1000;
static double ylo = 0.0;
static double yhi = MAX_THETA0;
static size_t ny = 1000;
static double *x;
static double *y;
static double *z;
static gsl_spline2d *spline;
static gsl_interp_accel *xacc;
static gsl_interp_accel *yacc;
static double x2lo = 0.0;
static double x2hi = 1.0;
static size_t nx2 = 1000;
static double *x2;
static double *y2;
/* Quantum efficiency weighted by the Cerenkov spectrum. */
static double weighted_qe;
static gsl_spline *spline2;
static gsl_interp_accel *xacc2;
static double prob_scatter(double wavelength, void *params)
{
/* Calculate the number of photons emitted per unit solid angle per cm at
* an angle `theta` for a particle travelling with velocity `beta` with an
* angular distribution of width `theta0`. */
double qe, delta, index;
double beta_cos_theta = ((double *) params)[0];
double beta_sin_theta_theta0 = ((double *) params)[1];
qe = get_quantum_efficiency(wavelength);
index = get_index_snoman_d2o(wavelength);
delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0;
return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI);
}
static double prob_scatter2(double wavelength, void *params)
{
/* Calculate the number of photons emitted per per cm for a particle
* travelling with velocity `beta`. */
double qe, index;
double beta = ((double *) params)[0];
qe = get_quantum_efficiency(wavelength);
index = get_index_snoman_d2o(wavelength);
return 2*M_PI*FINE_STRUCTURE_CONSTANT*(1-(1/(beta*beta*index*index)))*qe/pow(wavelength,2)*1e7;
}
double get_weighted_quantum_efficiency(void)
{
/* Returns the quantum efficiency weighted by the Cerenkov wavelength
* distribution, i.e. the probability that a photon randomly sampled from
* the Cerenkov wavelenght distribution between 200 and 800 nm is detected. */
if (!initialized) {
fprintf(stderr, "quantum efficiency hasn't been initialized!\n");
exit(1);
}
return weighted_qe;
}
static double gsl_quantum_efficiency(double wavelength, void *params)
{
/* Returns the quantum efficiency times the Cerenkov wavelength
* distribution. */
double qe;
qe = get_quantum_efficiency(wavelength);
return qe/pow(wavelength,2);
}
void init_interpolation(void)
{
size_t i, j;
double params[2];
double result, error;
size_t nevals;
int status;
gsl_integration_cquad_workspace *w;
gsl_function F;
x = malloc(nx*sizeof(double));
y = malloc(ny*sizeof(double));
z = malloc(nx*ny*sizeof(double));
spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, nx, ny);
xacc = gsl_interp_accel_alloc();
yacc = gsl_interp_accel_alloc();
for (i = 0; i < nx; i++) {
x[i] = xlo + (xhi-xlo)*i/(nx-1);
}
for (i = 0; i < ny; i++) {
y[i] = ylo + (yhi-ylo)*i/(ny-1);
}
w = gsl_integration_cquad_workspace_alloc(100);
F.function = &prob_scatter;
F.params = params;
for (i = 0; i < nx; i++) {
for (j = 0; j < ny; j++) {
params[0] = x[i];
params[1] = y[j];
status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
if (status) {
fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
exit(1);
}
gsl_spline2d_set(spline, z, i, j, result);
}
}
gsl_spline2d_init(spline, x, y, z, nx, ny);
x2 = malloc(nx2*sizeof(double));
y2 = malloc(nx2*sizeof(double));
spline2 = gsl_spline_alloc(gsl_interp_linear, nx2);
xacc2 = gsl_interp_accel_alloc();
for (i = 0; i < nx2; i++) {
x2[i] = x2lo + (x2hi-x2lo)*i/(nx2-1);
}
F.function = &prob_scatter2;
F.params = params;
for (i = 0; i < nx2; i++) {
params[0] = x2[i];
status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
if (status) {
fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
exit(1);
}
y2[i] = result;
}
gsl_spline_init(spline2, x2, y2, nx2);
F.function = &gsl_quantum_efficiency;
F.params = params;
status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
if (status) {
fprintf(stderr, "error integrating quantum efficiency distribution: %s\n", gsl_strerror(status));
exit(1);
}
weighted_qe = result*800/3.0;
initialized = 1;
gsl_integration_cquad_workspace_free(w);
}
double get_probability(double beta, double cos_theta, double theta0)
{
double sin_theta;
/* Technically this isn't defined up to a sign, but it doesn't matter since
* we are going to square it everywhere. */
sin_theta = sqrt(1-pow(cos_theta,2));
return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta);
}
double get_probability2(double beta)
{
return gsl_spline_eval(spline2, beta, xacc2);
}
void free_interpolation(void)
{
free(x);
free(y);
free(z);
gsl_spline2d_free(spline);
gsl_interp_accel_free(xacc);
gsl_interp_accel_free(yacc);
}
|
