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
|
#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 double xlo = -1.0;
static double xhi = 1.0;
static size_t nx = 200;
static double ylo = 0.0;
static double yhi = 1.0;
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;
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;
}
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);
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 = fabs(sin(acos(cos_theta)));
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);
}
|
