
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

/* 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 "scattering.h"
#include "quantum_efficiency.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_integration.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 "misc.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 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 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 FINE_STRUCTURE_CONSTANT*qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI)/beta_sin_theta_theta0;
}

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));

    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);
            }

            z[i*ny + j] = result;
        }
    }

    x2 = malloc(nx2*sizeof(double));
    y2 = malloc(nx2*sizeof(double));

    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;
    }

    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 sin_theta, double theta0)
{
    /* Make sure theta0 is less than MAX_THETA0, otherwise it's possible that
     * interp2d() will quit. */
    if (theta0 > MAX_THETA0) theta0 = MAX_THETA0;

    return beta*interp2d(beta*cos_theta, beta*sin_theta*theta0, x, y, z, nx, ny);
}

double get_probability2(double beta)
{
    return interp1d(beta, x2, y2, nx2);
}

void free_interpolation(void)
{
    free(x);
    free(y);
    free(z);

    free(x2);
    free(y2);
}