aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore9
-rw-r--r--src/electron.c195
-rw-r--r--src/electron.h5
-rw-r--r--src/fit.c85
-rw-r--r--src/likelihood.c318
-rw-r--r--src/likelihood.h10
-rw-r--r--src/misc.c8
-rw-r--r--src/misc.h1
-rw-r--r--src/muon.c34
-rw-r--r--src/muon.h1
-rw-r--r--src/proton.c34
-rw-r--r--src/proton.h1
-rw-r--r--src/scattering.c43
-rw-r--r--src/scattering.h1
-rw-r--r--src/solid_angle.c80
-rw-r--r--src/solid_angle.h3
-rw-r--r--src/test.c219
17 files changed, 928 insertions, 119 deletions
diff --git a/.gitignore b/.gitignore
index 15c3979..7a4d51c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,3 +13,12 @@ test-charge
*.txt
test-path
release.h
+*.cmd
+*.png
+*.aux
+*.log
+*.out
+*.pdf
+*.nav
+*.snm
+*.toc
diff --git a/src/electron.c b/src/electron.c
index e8e24d3..261ba39 100644
--- a/src/electron.c
+++ b/src/electron.c
@@ -15,12 +15,17 @@
#include "scattering.h"
#include "pmt_response.h"
#include "misc.h"
+#include <gsl/gsl_integration.h>
+#include <math.h> /* for fmax() */
static int initialized = 0;
-static double *x, *dEdx, *csda_range;
+static double *x, *dEdx_rad, *dEdx, *csda_range;
static size_t size;
+static gsl_interp_accel *acc_dEdx_rad;
+static gsl_spline *spline_dEdx_rad;
+
static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;
@@ -36,6 +41,165 @@ static gsl_spline *spline_range;
static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33;
static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33;
+void electron_get_position_distribution_parameters(double T0, double *a, double *b)
+{
+ /* Computes the gamma distribution parameters describing the longitudinal
+ * profile of the photons generated from an electromagnetic shower.
+ *
+ * From the PDG "Passage of Particles" section 32.5:
+ *
+ * "The mean longitudinal profile of the energy deposition in an
+ * electromagnetic cascade is reasonably well described by a gamma
+ * distribution."
+ *
+ * Here we use a slightly different form of the gamma distribution:
+ *
+ * f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a)
+ *
+ * I determined the b parameter by simulating high energy electrons using
+ * RAT-PAC and determined that it's roughly equal to the radiation length.
+ * To calculate the a parameter we use the formula from the PDG, i.e.
+ *
+ * tmax = (a-1)/b = ln(E/E_C) - 0.5
+ *
+ * Therefore, we calculate a as:
+ *
+ * a = tmax*b+1.
+ *
+ * `T` should be in units of MeV.
+ *
+ * Example:
+ *
+ * double a, b;
+ * electron_get_position_distribution_parameters(1000.0, &a, &b);
+ * double pdf = gamma_pdf(x, a, b);
+ *
+ * See http://pdg.lbl.gov/2014/reviews/rpp2014-rev-passage-particles-matter.pdf. */
+ double tmax;
+
+ *b = RADIATION_LENGTH;
+ tmax = log(T0/ELECTRON_CRITICAL_ENERGY_D2O) - 0.5;
+ *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1);
+}
+
+double electron_get_angular_distribution_alpha(double T0)
+{
+ /* To describe the angular distribution of photons in an electromagnetic
+ * shower I came up with the heuristic form:
+ *
+ * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
+ *
+ * I simulated high energy electrons using RAT-PAC in heavy water to fit
+ * for the alpha and beta parameters as a function of energy and determined
+ * that a reasonably good fit to the data was:
+ *
+ * alpha = c0 + c1/log(c2*T0 + c3)
+ *
+ * where T0 is the initial energy of the electron in MeV and c0, c1, c2,
+ * and c3 are constants which I fit out. */
+ return 3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00);
+}
+
+double electron_get_angular_distribution_beta(double T0)
+{
+ /* To describe the angular distribution of photons in an electromagnetic
+ * shower I came up with the heuristic form:
+ *
+ * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
+ *
+ * I simulated high energy electrons using RAT-PAC in heavy water to fit
+ * for the alpha and beta parameters as a function of energy and determined
+ * that a reasonably good fit to the data was:
+ *
+ * beta = c0 + c1/log(c2*T0 + c3)
+ *
+ * where T0 is the initial energy of the electron in MeV and c0, c1, c2,
+ * and c3 are constants which I fit out. */
+ return 1.35372e-01 + 2.22344e-01/log(1.96249e-02*T0 + 1.23912e+00);
+}
+
+double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double beta, double mu)
+{
+ return exp(-pow(fabs(cos_theta-mu),alpha)/beta);
+}
+
+static double gsl_electron_get_angular_pdf(double cos_theta, void *params)
+{
+ double alpha = ((double *) params)[0];
+ double beta = ((double *) params)[1];
+ double mu = ((double *) params)[2];
+ return electron_get_angular_pdf_no_norm(cos_theta,alpha,beta,mu);
+}
+
+static double electron_get_angular_pdf_norm(double alpha, double beta, double mu)
+{
+ double params[3];
+ gsl_integration_cquad_workspace *w;
+ gsl_function F;
+ double result, error;
+ int status;
+ size_t nevals;
+
+ w = gsl_integration_cquad_workspace_alloc(100);
+
+ F.function = &gsl_electron_get_angular_pdf;
+ F.params = params;
+
+ params[0] = alpha;
+ params[1] = beta;
+ params[2] = mu;
+
+ status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);
+
+ if (status) {
+ fprintf(stderr, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
+ exit(1);
+ }
+
+ gsl_integration_cquad_workspace_free(w);
+
+ return result;
+}
+
+double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu)
+{
+ /* Returns the probability density that a photon in the wavelength range
+ * 200 nm - 800 nm is emitted at an angle cos_theta.
+ *
+ * The angular distribution is modelled by the pdf:
+ *
+ * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta)
+ *
+ * where alpha and beta are constants which depend on the initial energy of
+ * the particle.
+ *
+ * There is no nice closed form solution for the integral of this function,
+ * so we just compute it on the fly. To make things faster, we keep track
+ * of the last alpha, beta, and mu parameters that were passed in so we
+ * avoid recomputing the normalization factor. */
+ size_t i;
+ static double last_alpha, last_beta, last_mu, norm;
+ static int first = 1;
+ static double x[10000], y[10000];
+
+ if (first || alpha != last_alpha || beta != last_beta || mu != last_mu) {
+ norm = electron_get_angular_pdf_norm(alpha, beta, mu);
+
+ last_alpha = alpha;
+ last_beta = beta;
+ last_mu = mu;
+
+ for (i = 0; i < LEN(x); i++) {
+ x[i] = -1.0 + 2.0*i/(LEN(x)-1);
+ y[i] = electron_get_angular_pdf_no_norm(x[i], alpha, beta, mu)/norm;
+ }
+
+ first = 0;
+ }
+
+ return interp1d(cos_theta,x,y,LEN(x));
+}
+
static int init()
{
int i, j;
@@ -80,6 +244,7 @@ static int init()
}
x = malloc(sizeof(double)*n);
+ dEdx_rad = malloc(sizeof(double)*n);
dEdx = malloc(sizeof(double)*n);
csda_range = malloc(sizeof(double)*n);
size = n;
@@ -112,6 +277,9 @@ static int init()
case 0:
x[n] = value;
break;
+ case 2:
+ dEdx_rad[n] = value;
+ break;
case 3:
dEdx[n] = value;
break;
@@ -128,6 +296,10 @@ static int init()
fclose(f);
+ acc_dEdx_rad = gsl_interp_accel_alloc();
+ spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
+ gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);
+
acc_dEdx = gsl_interp_accel_alloc();
spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
gsl_spline_init(spline_dEdx, x, dEdx, size);
@@ -167,6 +339,27 @@ double electron_get_range(double T, double rho)
return gsl_spline_eval(spline_range, T, acc_range)/rho;
}
+double electron_get_dEdx_rad(double T, double rho)
+{
+ /* Returns the approximate radiative dE/dx for a electron in water with
+ * kinetic energy `T`.
+ *
+ * `T` should be in MeV and `rho` in g/cm^3.
+ *
+ * Return value is in MeV/cm.
+ *
+ * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+ if (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ if (T < spline_dEdx_rad->x[0]) return spline_dEdx_rad->y[0];
+
+ return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
+}
+
double electron_get_dEdx(double T, double rho)
{
/* Returns the approximate dE/dx for a electron in water with kinetic energy
diff --git a/src/electron.h b/src/electron.h
index 1bb589f..cb3f24e 100644
--- a/src/electron.h
+++ b/src/electron.h
@@ -3,7 +3,12 @@
#include <stddef.h> /* for size_t */
+void electron_get_position_distribution_parameters(double T0, double *a, double *b);
+double electron_get_angular_distribution_alpha(double T0);
+double electron_get_angular_distribution_beta(double T0);
+double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu);
double electron_get_range(double T, double rho);
+double electron_get_dEdx_rad(double T, double rho);
double electron_get_dEdx(double T, double rho);
#endif
diff --git a/src/fit.c b/src/fit.c
index eba125e..82586f9 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -22,6 +22,7 @@
#include <signal.h> /* for signal() */
#include "release.h"
#include "id_particles.h"
+#include "misc.h"
char *GitSHA1(void);
char *GitDirty(void);
@@ -40,8 +41,14 @@ typedef struct fitParams {
int n;
int fast;
int print;
+ int id;
} fitParams;
+int particles[] = {
+ IDP_E_MINUS,
+ IDP_MU_MINUS,
+};
+
/* In order to start the fitter close to the minimum, we first do a series of
* "quick" minimizations starting at the following points. We keep track of the
* parameters with the best likelihood value and then start the "real"
@@ -4996,7 +5003,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
z2[0] = x[8];
gettimeofday(&tv_start, NULL);
- fval = nll_muon(fpars->ev, IDP_MU_MINUS, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast);
+ fval = nll_muon(fpars->ev, fpars->id, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast);
gettimeofday(&tv_stop, NULL);
long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -5099,11 +5106,11 @@ void guess_direction(event *ev, double *pos, double *theta, double *phi)
*phi = atan2(dir[1],dir[0]);
}
-int fit_event(event *ev, double *xopt, double *fmin)
+int fit_event(event *ev, double *xopt, double *fmin, int id)
{
size_t i;
fitParams fpars;
- double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin;
+ double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin, mass;
struct timeval tv_start, tv_stop;
opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
@@ -5124,12 +5131,29 @@ int fit_event(event *ev, double *xopt, double *fmin)
n = get_index_snoman_d2o(400.0);
+ switch (id) {
+ case IDP_E_MINUS:
+ case IDP_E_PLUS:
+ mass = ELECTRON_MASS;
+ break;
+ case IDP_MU_MINUS:
+ case IDP_MU_PLUS:
+ mass = MUON_MASS;
+ break;
+ case IDP_PROTON:
+ mass = PROTON_MASS;
+ break;
+ default:
+ fprintf(stderr, "unknown particle id %i\n", id);
+ exit(1);
+ }
+
/* Calculate the Cerenkov threshold for a muon. */
- Tmin = MUON_MASS/sqrt(1.0-1.0/(n*n)) - MUON_MASS;
+ Tmin = mass/sqrt(1.0-1.0/(n*n)) - mass;
/* If our guess is below the Cerenkov threshold, start at the Cerenkov
* threshold. */
- double Tmin2 = sqrt(MUON_MASS*MUON_MASS/(1-pow(0.9,2)))-MUON_MASS;
+ double Tmin2 = sqrt(mass*mass/(1-pow(0.9,2)))-mass;
if (T0 < Tmin2) T0 = Tmin2;
x0[0] = T0;
@@ -5178,6 +5202,7 @@ int fit_event(event *ev, double *xopt, double *fmin)
nlopt_set_initial_step(opt, ss);
fpars.ev = ev;
+ fpars.id = id;
iter = 0;
@@ -5538,32 +5563,38 @@ int main(int argc, char **argv)
rv = get_event(f,&ev,&b);
- gettimeofday(&tv_start, NULL);
- if (fit_event(&ev,xopt,&fmin) == NLOPT_FORCED_STOP) {
- printf("ctrl-c caught. quitting...\n");
- goto end;
- }
- gettimeofday(&tv_stop, NULL);
-
- long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
-
if (fout) {
if (first_ev) fprintf(fout, " ev:\n");
fprintf(fout, " - gtid: %i\n", ev.gtid);
fprintf(fout, " fit:\n");
- fprintf(fout, " -\n");
- fprintf(fout, " energy: %.2f\n", xopt[0]);
- fprintf(fout, " posx: %.2f\n", xopt[1]);
- fprintf(fout, " posy: %.2f\n", xopt[2]);
- fprintf(fout, " posz: %.2f\n", xopt[3]);
- fprintf(fout, " theta: %.4f\n", xopt[4]);
- fprintf(fout, " phi: %.4f\n", xopt[5]);
- fprintf(fout, " t0: %.2f\n", xopt[6]);
- fprintf(fout, " z1: %.2f\n", xopt[7]);
- fprintf(fout, " z2: %.2f\n", xopt[8]);
- fprintf(fout, " fmin: %.2f\n", fmin);
- fprintf(fout, " time: %lld\n", elapsed);
- fflush(fout);
+ }
+
+ for (i = 0; i < LEN(particles); i++) {
+ gettimeofday(&tv_start, NULL);
+ if (fit_event(&ev,xopt,&fmin,particles[i]) == NLOPT_FORCED_STOP) {
+ printf("ctrl-c caught. quitting...\n");
+ goto end;
+ }
+ gettimeofday(&tv_stop, NULL);
+
+ long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
+
+ if (fout) {
+ fprintf(fout, " -\n");
+ fprintf(fout, " id: %i\n", particles[i]);
+ fprintf(fout, " energy: %.2f\n", xopt[0]);
+ fprintf(fout, " posx: %.2f\n", xopt[1]);
+ fprintf(fout, " posy: %.2f\n", xopt[2]);
+ fprintf(fout, " posz: %.2f\n", xopt[3]);
+ fprintf(fout, " theta: %.4f\n", xopt[4]);
+ fprintf(fout, " phi: %.4f\n", xopt[5]);
+ fprintf(fout, " t0: %.2f\n", xopt[6]);
+ fprintf(fout, " z1: %.2f\n", xopt[7]);
+ fprintf(fout, " z2: %.2f\n", xopt[8]);
+ fprintf(fout, " fmin: %.2f\n", fmin);
+ fprintf(fout, " time: %lld\n", elapsed);
+ fflush(fout);
+ }
}
first_ev = 0;
diff --git a/src/likelihood.c b/src/likelihood.c
index 2340fb1..f769216 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -21,6 +21,7 @@
#include "electron.h"
#include "proton.h"
#include "id_particles.h"
+#include <gsl/gsl_cdf.h>
particle *particle_init(int id, double T0, size_t n)
{
@@ -46,6 +47,7 @@ particle *particle_init(int id, double T0, size_t n)
p->x = malloc(sizeof(double)*n);
p->T = malloc(sizeof(double)*n);
p->n = n;
+ p->rad = 0.0;
p->x[0] = 0;
p->T[0] = T0;
@@ -64,6 +66,8 @@ particle *particle_init(int id, double T0, size_t n)
p->x[i] = p->range*i/(n-1);
dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY);
p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ dEdx = electron_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
+ p->rad += dEdx*(p->x[i]-p->x[i-1]);
}
/* Make sure that the energy is zero at the last step. This is so that
* when we try to bisect the point along the track where the speed of
@@ -91,6 +95,8 @@ particle *particle_init(int id, double T0, size_t n)
p->x[i] = p->range*i/(n-1);
dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY);
p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ dEdx = muon_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
+ p->rad += dEdx*(p->x[i]-p->x[i-1]);
}
/* Make sure that the energy is zero at the last step. This is so that
* when we try to bisect the point along the track where the speed of
@@ -114,6 +120,8 @@ particle *particle_init(int id, double T0, size_t n)
p->x[i] = p->range*i/(n-1);
dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY);
p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ dEdx = proton_get_dEdx_rad(p->T[i-1], WATER_DENSITY);
+ p->rad += dEdx*(p->x[i]-p->x[i-1]);
}
/* Make sure that the energy is zero at the last step. This is so that
* when we try to bisect the point along the track where the speed of
@@ -155,6 +163,50 @@ void particle_free(particle *p)
free(p);
}
+static double get_expected_charge_shower(double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double alpha, double beta, double rad)
+{
+ double pmt_dir[3], cos_theta, n, omega, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge;
+
+ SUB(pmt_dir,pmt_pos,pos);
+
+ normalize(pmt_dir);
+
+ cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+
+ *reflected = 0.0;
+ if (cos_theta_pmt > 0) return 0.0;
+
+ /* Calculate the cosine of the angle between the track direction and the
+ * vector to the PMT. */
+ cos_theta = DOT(dir,pmt_dir);
+
+ omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
+
+ R = NORM(pos);
+
+ if (R <= AV_RADIUS) {
+ n = n_d2o;
+ } else {
+ n = n_h2o;
+ }
+
+ theta_pmt = acos(-cos_theta_pmt);
+ f_reflec = get_weighted_pmt_reflectivity(theta_pmt);
+ f = get_weighted_pmt_response(theta_pmt);
+
+ absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
+ absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
+ absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
+
+ l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
+
+ charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI);
+
+ *reflected = f_reflec*charge;
+
+ return f*charge;
+}
+
static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o)
{
double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge;
@@ -174,7 +226,7 @@ static double get_expected_charge(double x, double beta, double theta0, double *
* vector to the PMT. */
cos_theta = DOT(dir,pmt_dir);
- omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
+ omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r);
R = NORM(pos);
@@ -204,7 +256,7 @@ static double get_expected_charge(double x, double beta, double theta0, double *
return f*charge;
}
-double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
+double F(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
{
/* Returns the CDF for the time distribution of photons at time `t`. */
size_t i;
@@ -217,11 +269,15 @@ double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_
p += mu_direct[i]*norm_cdf(t,ts[i],sigma);
mu_total += mu_direct[i];
}
+ for (i = 0; i < n; i++) {
+ p += mu_shower[i]*norm_cdf(t,ts_shower[i],ts_sigma[i]);
+ mu_total += mu_shower[i];
+ }
return p/mu_total;
}
-double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
+double f(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
{
/* Returns the probability that a photon is detected at time `t`.
*
@@ -247,11 +303,15 @@ double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_
p += mu_direct[i]*norm(t,ts[i],sigma);
mu_total += mu_direct[i];
}
+ for (i = 0; i < n; i++) {
+ p += mu_shower[i]*norm(t,ts_shower[i],ts_sigma[i]);
+ mu_total += mu_shower[i];
+ }
return p/mu_total;
}
-double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, size_t n2, double *ts, double tmean, double sigma)
+double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n2, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
{
/* Returns the first order statistic for observing a PMT hit at time `t`
* given `n` hits.
@@ -263,18 +323,73 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
* different transit times across the face of the PMT, it seems better to
* convolve first which is what we do here. In addition, the problem is not
* analytically tractable if you do things the other way around. */
- return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
+ return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+}
+
+static void integrate_path_shower(double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad, double pos_a, double pos_b)
+{
+ size_t i;
+ static double qs[MAX_NPOINTS];
+ static double q_indirects[MAX_NPOINTS];
+ static double ts[MAX_NPOINTS];
+ static double ts2[MAX_NPOINTS];
+ double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, alpha, beta;
+
+ if (n > MAX_NPOINTS) {
+ fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
+ exit(1);
+ }
+
+ alpha = electron_get_angular_distribution_alpha(T0);
+ beta = electron_get_angular_distribution_beta(T0);
+
+ /* FIXME: I just calculate delta assuming 400 nm light. */
+ wavelength0 = 400.0;
+ n_d2o = get_index_snoman_d2o(wavelength0);
+ n_h2o = get_index_snoman_h2o(wavelength0);
+
+ for (i = 0; i < n; i++) {
+ pos[0] = pos0[0] + x[i]*dir0[0];
+ pos[1] = pos0[1] + x[i]*dir0[1];
+ pos[2] = pos0[2] + x[i]*dir0[2];
+
+ get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+ t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+
+ qs[i] = get_expected_charge_shower(pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, alpha, beta, rad)*pdf[i];
+ q_indirects[i] *= pdf[i];
+ ts[i] = t*qs[i];
+ ts2[i] = t*t*qs[i];
+ }
+
+ *mu_direct = trapz(qs,(b-a)/(n-1),n);
+ *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n);
+ *time = trapz(ts,(b-a)/(n-1),n)/(*mu_direct);
+ double t2 = trapz(ts2,(b-a)/(n-1),n)/(*mu_direct);
+
+ *sigma = sqrt(t2-(*time)*(*time));
+
+ if (*mu_direct == 0.0) {
+ *time = 0.0;
+ *sigma = PMT_TTS;
+ } else {
+ *sigma = sqrt((*sigma)*(*sigma) + PMT_TTS*PMT_TTS);
+ }
}
static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time)
{
size_t i;
- double *qs, *q_indirects, *ts;
+ static double qs[MAX_NPOINTS];
+ static double q_indirects[MAX_NPOINTS];
+ static double ts[MAX_NPOINTS];
double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;
- qs = malloc(sizeof(double)*n);
- q_indirects = malloc(sizeof(double)*n);
- ts = malloc(sizeof(double)*n);
+ if (n > MAX_NPOINTS) {
+ fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n");
+ exit(1);
+ }
/* FIXME: I just calculate delta assuming 400 nm light. */
wavelength0 = 400.0;
@@ -297,10 +412,6 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl
*mu_direct = trapz(qs,(b-a)/(n-1),n);
*mu_indirect = trapz(q_indirects,(b-a)/(n-1),n);
*time = trapz(ts,(b-a)/(n-1),n);
-
- free(qs);
- free(q_indirects);
- free(ts);
}
static double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
@@ -424,7 +535,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
sin_theta = sqrt(1-pow(cos_theta,2));
/* Get the solid angle of the PMT at the position `s` along the track. */
- omega = get_solid_angle_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);
+ omega = get_solid_angle_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS);
theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
@@ -526,7 +637,8 @@ static double getKineticEnergy(double x, void *p)
double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, size_t npoints, int fast)
{
size_t i, j, nhit;
- double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
+ static double logp[MAX_PE], nll[MAX_PMTS];
+ double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
particle *p;
double pmt_dir[3], R, cos_theta, sin_theta, wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov, s, l_h2o, l_d2o;
double logp_path;
@@ -535,11 +647,31 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
double mu_reflected;
path *path;
- double mu_direct[MAX_PMTS];
- double ts[MAX_PMTS];
- double mu[MAX_PMTS];
+ static double mu_direct[MAX_PMTS];
+ static double mu_shower[MAX_PMTS] = {0};
+ static double ts[MAX_PMTS];
+ static double ts_shower[MAX_PMTS] = {0};
+ static double ts_sigma[MAX_PMTS] = {0};
+ static double mu[MAX_PMTS];
double mu_noise, mu_indirect;
+ for (i = 0; i < MAX_PMTS; i++) ts_sigma[i] = PMT_TTS;
+
+ double pos_a, pos_b;
+
+ electron_get_position_distribution_parameters(T0, &pos_a, &pos_b);
+
+ /* Upper limit to integrate for the electromagnetic shower. */
+ double xlo = 0.0;
+ double xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b);
+
+ double *x = malloc(sizeof(double)*npoints);
+ double *pdf = malloc(sizeof(double)*npoints);
+ for (i = 0; i < npoints; i++) {
+ x[i] = xlo + i*(xhi-xlo)/(npoints-1);
+ pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
+ }
+
double result;
double min_ratio;
@@ -582,80 +714,87 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
ts[i] += t0;
mu_indirect += mu_reflected;
+ continue;
+ }
+
+ /* First, we try to compute the distance along the track where the
+ * PMT is at the Cerenkov angle. The reason for this is because for
+ * heavy particles like muons which don't scatter much, the
+ * probability distribution for getting a photon hit along the
+ * track looks kind of like a delta function, i.e. the PMT is only
+ * hit over a very narrow window when the angle between the track
+ * direction and the PMT is *very* close to the Cerenkov angle
+ * (it's not a perfect delta function since there is some width due
+ * to dispersion). In this case, it's possible that the numerical
+ * integration completely skips over the delta function and so
+ * predicts an expected charge of 0.
+ *
+ * To fix this, we only integrate 1 meter on both sides of the
+ * point along the track where the PMT is at the Cerenkov angle. */
+
+ /* First, we find the point along the track where the PMT is at the
+ * Cerenkov angle. */
+ SUB(pmt_dir,pmts[i].pos,pos);
+ /* Compute the distance to the PMT. */
+ R = NORM(pmt_dir);
+ normalize(pmt_dir);
+
+ /* Calculate the cosine of the angle between the track direction and the
+ * vector to the PMT at the start of the track. */
+ cos_theta = DOT(dir,pmt_dir);
+ sin_theta = sqrt(1-pow(cos_theta,2));
+
+ /* Now, we compute the distance along the track where the PMT is at the
+ * Cerenkov angle.
+ *
+ * Note: This formula comes from using the "Law of sines" where the three
+ * vertices of the triangle are the starting position of the track, the
+ * point along the track that we want to find, and the PMT position. */
+ s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
+
+ /* Make sure that the point is somewhere along the track between 0 and
+ * `smax`. */
+ if (s < 0) s = 0.0;
+ else if (s > smax) s = smax;
+
+ /* Compute the endpoints of the integration. */
+ a = s - 100.0;
+ b = s + 100.0;
+
+ if (a < 0.0) a = 0.0;
+ if (b > range) b = range;
+
+ double q_indirect;
+ integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result);
+ mu_indirect += q_indirect;
+
+ if (mu_direct[i] > 1e-9) {
+ ts[i] = t0 + result/mu_direct[i];
} else {
- /* First, we try to compute the distance along the track where the
- * PMT is at the Cerenkov angle. The reason for this is because for
- * heavy particles like muons which don't scatter much, the
- * probability distribution for getting a photon hit along the
- * track looks kind of like a delta function, i.e. the PMT is only
- * hit over a very narrow window when the angle between the track
- * direction and the PMT is *very* close to the Cerenkov angle
- * (it's not a perfect delta function since there is some width due
- * to dispersion). In this case, it's possible that the numerical
- * integration completely skips over the delta function and so
- * predicts an expected charge of 0.
- *
- * To fix this, we only integrate 1 meter on both sides of the
- * point along the track where the PMT is at the Cerenkov angle. */
-
- /* First, we find the point along the track where the PMT is at the
- * Cerenkov angle. */
- SUB(pmt_dir,pmts[i].pos,pos);
- /* Compute the distance to the PMT. */
- R = NORM(pmt_dir);
- normalize(pmt_dir);
-
- /* Calculate the cosine of the angle between the track direction and the
- * vector to the PMT at the start of the track. */
- cos_theta = DOT(dir,pmt_dir);
- sin_theta = sqrt(1-pow(cos_theta,2));
-
- /* Now, we compute the distance along the track where the PMT is at the
- * Cerenkov angle.
- *
- * Note: This formula comes from using the "Law of sines" where the three
- * vertices of the triangle are the starting position of the track, the
- * point along the track that we want to find, and the PMT position. */
- s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov;
-
- /* Make sure that the point is somewhere along the track between 0 and
- * `smax`. */
- if (s < 0) s = 0.0;
- else if (s > smax) s = smax;
-
- /* Compute the endpoints of the integration. */
- a = s - 100.0;
- b = s + 100.0;
-
- if (a < 0.0) a = 0.0;
- if (b > range) b = range;
-
- double q_indirect;
- integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result);
- mu_indirect += q_indirect;
+ /* If the expected number of PE is very small, our estimate of
+ * the time can be crazy since the error in the total charge is
+ * in the denominator. Therefore, we estimate the time using
+ * the same technique as in get_total_charge_approx(). See that
+ * function for more details. */
+ n_h2o = get_index_snoman_h2o(wavelength0);
- if (mu_direct[i] > 1e-9) {
- ts[i] = t0 + result/mu_direct[i];
- } else {
- /* If the expected number of PE is very small, our estimate of
- * the time can be crazy since the error in the total charge is
- * in the denominator. Therefore, we estimate the time using
- * the same technique as in get_total_charge_approx(). See that
- * function for more details. */
- n_h2o = get_index_snoman_h2o(wavelength0);
+ /* Compute the position of the particle at a distance `s` along the track. */
+ tmp[0] = pos[0] + s*dir[0];
+ tmp[1] = pos[1] + s*dir[1];
+ tmp[2] = pos[2] + s*dir[2];
- /* Compute the position of the particle at a distance `s` along the track. */
- tmp[0] = pos[0] + s*dir[0];
- tmp[1] = pos[1] + s*dir[1];
- tmp[2] = pos[2] + s*dir[2];
+ SUB(pmt_dir,pmts[i].pos,tmp);
- SUB(pmt_dir,pmts[i].pos,tmp);
+ get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
- get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
+ /* Assume the particle is travelling at the speed of light. */
+ ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+ }
- /* Assume the particle is travelling at the speed of light. */
- ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
- }
+ if (id == IDP_E_MINUS || id == IDP_E_PLUS) {
+ integrate_path_shower(x,pdf,T0,pos,dir,i,0,xhi,npoints,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad,pos_a,pos_b);
+ mu_indirect += q_indirect;
+ ts_shower[i] = t0 + result;
}
}
@@ -677,7 +816,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
if (fast)
mu[i] = mu_direct[i] + mu_noise;
else
- mu[i] = mu_direct[i] + mu_indirect + mu_noise;
+ mu[i] = mu_direct[i] + mu_shower[i] + mu_indirect + mu_noise;
}
if (fast)
@@ -697,7 +836,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
if (ev->pmt_hits[i].hit) {
for (j = 1; j < MAX_PE; j++) {
- logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], ts[i], PMT_TTS);
+ logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
if (logp[j] - max_logp < min_ratio*ln(10)) {
@@ -724,5 +863,8 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
for (i = 0; i < n; i++)
logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1);
+ free(x);
+ free(pdf);
+
return kahan_sum(nll,nhit) - logp_path;
}
diff --git a/src/likelihood.h b/src/likelihood.h
index b326d5d..7622ffe 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -4,6 +4,8 @@
#include "event.h"
#include <stddef.h> /* for size_t */
+#define MAX_NPOINTS 1000
+
/* Maximum number of PE to consider if the PMT was hit. */
#define MAX_PE 10000
/* Maximum number of PE to consider if the PMT wasn't hit.
@@ -35,9 +37,17 @@
#define GTVALID 400.0
#define BETA_MIN 0.8
+/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy
+ * lost to radiation for electrons.
+ *
+ * FIXME: This is just a rough estimate, should use an energy dependent
+ * quantity from simulation. */
+#define PHOTONS_PER_MEV 400.0
+
typedef struct particle {
double mass;
double range;
+ double rad;
double *x;
double *T;
size_t n;
diff --git a/src/misc.c b/src/misc.c
index 605f3ea..d5667e8 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -498,3 +498,11 @@ double std(const double *x, size_t n)
return sqrt(sum/n);
}
+
+double gamma_pdf(double x, double k, double theta)
+{
+ /* Returns the PDF for the gamma distribution.
+ *
+ * See https://en.wikipedia.org/wiki/Gamma_distribution. */
+ return pow(x,k-1)*exp(-x/theta)/(gsl_sf_gamma(k)*pow(theta,k));
+}
diff --git a/src/misc.h b/src/misc.h
index 3a55c2f..94ccee7 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -26,5 +26,6 @@ double norm(double x, double mu, double sigma);
double norm_cdf(double x, double mu, double sigma);
double mean(const double *x, size_t n);
double std(const double *x, size_t n);
+double gamma_pdf(double x, double k, double theta);
#endif
diff --git a/src/muon.c b/src/muon.c
index 02cf55a..76e5223 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -18,9 +18,12 @@
static int initialized = 0;
-static double *x, *dEdx, *csda_range;
+static double *x, *dEdx_rad, *dEdx, *csda_range;
static size_t size;
+static gsl_interp_accel *acc_dEdx_rad;
+static gsl_spline *spline_dEdx_rad;
+
static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;
@@ -83,6 +86,7 @@ static int init()
}
x = malloc(sizeof(double)*n);
+ dEdx_rad = malloc(sizeof(double)*n);
dEdx = malloc(sizeof(double)*n);
csda_range = malloc(sizeof(double)*n);
size = n;
@@ -117,6 +121,9 @@ static int init()
case 0:
x[n] = value;
break;
+ case 6:
+ dEdx_rad[n] = value;
+ break;
case 7:
dEdx[n] = value;
break;
@@ -133,6 +140,10 @@ static int init()
fclose(f);
+ acc_dEdx_rad = gsl_interp_accel_alloc();
+ spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
+ gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);
+
acc_dEdx = gsl_interp_accel_alloc();
spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
gsl_spline_init(spline_dEdx, x, dEdx, size);
@@ -172,6 +183,27 @@ double muon_get_range(double T, double rho)
return gsl_spline_eval(spline_range, T, acc_range)/rho;
}
+double muon_get_dEdx_rad(double T, double rho)
+{
+ /* Returns the approximate radiative dE/dx for a muon in water with kinetic
+ * energy `T`.
+ *
+ * `T` should be in MeV and `rho` in g/cm^3.
+ *
+ * Return value is in MeV/cm.
+ *
+ * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+ if (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ if (T < spline_dEdx_rad->x[0]) return spline_dEdx_rad->y[0];
+
+ return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
+}
+
double muon_get_dEdx(double T, double rho)
{
/* Returns the approximate dE/dx for a muon in water with kinetic energy
diff --git a/src/muon.h b/src/muon.h
index 7717501..586d67d 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -6,6 +6,7 @@
#define EULER_CONSTANT 0.57721
double muon_get_range(double T, double rho);
+double muon_get_dEdx_rad(double T, double rho);
double muon_get_dEdx(double T, double rho);
#endif
diff --git a/src/proton.c b/src/proton.c
index 922c552..f5a6b61 100644
--- a/src/proton.c
+++ b/src/proton.c
@@ -18,9 +18,12 @@
static int initialized = 0;
-static double *x, *dEdx, *csda_range;
+static double *x, *dEdx_rad, *dEdx, *csda_range;
static size_t size;
+static gsl_interp_accel *acc_dEdx_rad;
+static gsl_spline *spline_dEdx_rad;
+
static gsl_interp_accel *acc_dEdx;
static gsl_spline *spline_dEdx;
@@ -71,6 +74,7 @@ static int init()
}
x = malloc(sizeof(double)*n);
+ dEdx_rad = malloc(sizeof(double)*n);
dEdx = malloc(sizeof(double)*n);
csda_range = malloc(sizeof(double)*n);
size = n;
@@ -103,6 +107,9 @@ static int init()
case 0:
x[n] = value;
break;
+ case 2:
+ dEdx_rad[n] = value;
+ break;
case 3:
dEdx[n] = value;
break;
@@ -119,6 +126,10 @@ static int init()
fclose(f);
+ acc_dEdx_rad = gsl_interp_accel_alloc();
+ spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size);
+ gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size);
+
acc_dEdx = gsl_interp_accel_alloc();
spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size);
gsl_spline_init(spline_dEdx, x, dEdx, size);
@@ -158,6 +169,27 @@ double proton_get_range(double T, double rho)
return gsl_spline_eval(spline_range, T, acc_range)/rho;
}
+double proton_get_dEdx_rad(double T, double rho)
+{
+ /* Returns the approximate radiative dE/dx for a proton in water with
+ * kinetic energy `T`.
+ *
+ * `T` should be in MeV and `rho` in g/cm^3.
+ *
+ * Return value is in MeV/cm.
+ *
+ * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+ if (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ if (T < spline_dEdx_rad->x[0]) return spline_dEdx_rad->y[0];
+
+ return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho;
+}
+
double proton_get_dEdx(double T, double rho)
{
/* Returns the approximate dE/dx for a proton in water with kinetic energy
diff --git a/src/proton.h b/src/proton.h
index 8e19e9c..008f835 100644
--- a/src/proton.h
+++ b/src/proton.h
@@ -4,6 +4,7 @@
#include <stddef.h> /* for size_t */
double proton_get_range(double T, double rho);
+double proton_get_dEdx_rad(double T, double rho);
double proton_get_dEdx(double T, double rho);
#endif
diff --git a/src/scattering.c b/src/scattering.c
index c4bdd57..68256dc 100644
--- a/src/scattering.c
+++ b/src/scattering.c
@@ -14,6 +14,8 @@
#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;
@@ -36,6 +38,9 @@ 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;
@@ -73,6 +78,30 @@ static double prob_scatter2(double wavelength, void *params)
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;
@@ -149,6 +178,20 @@ void init_interpolation(void)
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);
}
diff --git a/src/scattering.h b/src/scattering.h
index 65765b0..69d4c96 100644
--- a/src/scattering.h
+++ b/src/scattering.h
@@ -9,6 +9,7 @@
void init_interpolation(void);
double get_probability(double beta, double cos_theta, double theta0);
double get_probability2(double beta);
+double get_weighted_quantum_efficiency(void);
void free_interpolation(void);
#endif
diff --git a/src/solid_angle.c b/src/solid_angle.c
index 7201533..354547e 100644
--- a/src/solid_angle.c
+++ b/src/solid_angle.c
@@ -4,6 +4,7 @@
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include "vector.h"
+#include "misc.h"
static double hd[11] = {0.1,0.125,0.150,0.175,0.2,0.3,0.4,0.5,0.6,0.8,1.0};
static double Rd[13] = {0.1,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0};
@@ -40,7 +41,8 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r)
* regimes, the approximation is not very good and so it is modified by a
* correction factor which comes from a lookup table.
*
- * This formula comes from "The Solid Angle Subtended at a Point by a * Circular Disk." Gardener et al. 1969. */
+ * This formula comes from "The Solid Angle Subtended at a Point by a
+ * Circular Disk." Gardener et al. 1969. */
double dir[3];
double h, r0, D, a, u1, u2, F;
static gsl_spline2d *spline;
@@ -78,6 +80,82 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r)
}
}
+double get_solid_angle2(double L2, double r2)
+{
+ double k = 2*sqrt(r2/(pow(L2,2)+pow(1+r2,2)));
+ double a2 = 4*r2/pow(1+r2,2);
+ double L_Rmax = L2/sqrt(pow(L2,2)+pow(1+r2,2));
+ double r0_r = (r2-1)/(r2+1);
+ double K, P;
+
+ if (fabs(r2-1) < 1e-5) {
+ /* If r0 is very close to r, the GSL routines below will occasionally
+ * throw a domain error. */
+ K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+ return M_PI - 2*L_Rmax*K;
+ } else if (r2 <= 1) {
+ K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+ P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE);
+ return 2*M_PI - 2*L_Rmax*K + (2*L_Rmax)*r0_r*P;
+ } else {
+ K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+ P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE);
+ return -2*L_Rmax*K + (2*L_Rmax)*r0_r*P;
+ }
+}
+
+double get_solid_angle_lookup(double L2, double r2)
+{
+ size_t i, j;
+ static int initialized = 0;
+ static double xs[1000];
+ static double ys[1000];
+ static double zs[1000*1000];
+ static double xlo = 0.0;
+ static double xhi = 1000;
+ static double ylo = 0.0;
+ static double yhi = 1000;
+
+ if (!initialized) {
+ for (i = 0; i < LEN(xs); i++) {
+ xs[i] = xlo + (xhi-xlo)*i/(LEN(xs)-1);
+ }
+ for (i = 0; i < LEN(ys); i++) {
+ ys[i] = ylo + (yhi-ylo)*i/(LEN(ys)-1);
+ }
+ for (i = 0; i < LEN(xs); i++) {
+ for (j = 0; j < LEN(ys); j++) {
+ zs[i*LEN(ys) + j] = get_solid_angle2(xs[i],ys[j]);
+ }
+ }
+ initialized = 1;
+ }
+
+ if (L2 < xlo || L2 > xhi || r2 < ylo || r2 > yhi) {
+ /* If the arguments are out of bounds of the lookup table, just call
+ * get_solid_angle2(). */
+ return get_solid_angle2(L2,r2);
+ }
+
+ return interp2d(L2,r2,xs,ys,zs,LEN(xs),LEN(ys));
+}
+
+double get_solid_angle_fast(double *pos, double *pmt, double *n, double r)
+{
+ double dir[3];
+ double L, r0, R;
+
+ dir[0] = pos[0] - pmt[0];
+ dir[1] = pos[1] - pmt[1];
+ dir[2] = pos[2] - pmt[2];
+
+ L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+ R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+ r0 = sqrt(R*R - L*L);
+
+ return get_solid_angle_lookup(L/r,r0/r);
+}
+
double get_solid_angle(double *pos, double *pmt, double *n, double r)
{
/* Returns the solid angle subtended by a circular disk of radius r at a
diff --git a/src/solid_angle.h b/src/solid_angle.h
index aacb692..d0c287c 100644
--- a/src/solid_angle.h
+++ b/src/solid_angle.h
@@ -4,6 +4,9 @@
#include <math.h>
double get_solid_angle_approx(double *pos, double *pmt, double *n, double r);
+double get_solid_angle2(double L2, double r2);
+double get_solid_angle_lookup(double L2, double r2);
+double get_solid_angle_fast(double *pos, double *pmt, double *n, double r);
double get_solid_angle(double *pos, double *pmt, double *n, double r);
#endif
diff --git a/src/test.c b/src/test.c
index f8cedb9..7179f6b 100644
--- a/src/test.c
+++ b/src/test.c
@@ -19,6 +19,8 @@
#include "likelihood.h"
#include "id_particles.h"
#include <gsl/gsl_spline2d.h>
+#include <gsl/gsl_errno.h> /* for gsl_strerror() */
+#include <gsl/gsl_randist.h>
typedef int testFunction(char *err);
@@ -970,6 +972,218 @@ err:
return 1;
}
+static double gsl_electron_get_angular_pdf(double cos_theta, void *params)
+{
+ double alpha = ((double *) params)[0];
+ double beta = ((double *) params)[1];
+ double mu = ((double *) params)[2];
+ return electron_get_angular_pdf(cos_theta,alpha,beta,mu);
+}
+
+int test_electron_get_angular_pdf(char *err)
+{
+ /* Tests that the electron_get_angular_pdf() function integrates to 1. */
+ size_t i;
+ double params[3];
+ gsl_integration_cquad_workspace *w;
+ gsl_function F;
+ double result, error;
+ int status;
+ size_t nevals;
+ double T0;
+
+ w = gsl_integration_cquad_workspace_alloc(100);
+
+ F.function = &gsl_electron_get_angular_pdf;
+ F.params = params;
+
+ init_genrand(0);
+
+ for (i = 0; i < 100; i++) {
+ T0 = genrand_real2()*1000;
+
+ params[0] = electron_get_angular_distribution_alpha(T0);
+ params[1] = electron_get_angular_distribution_beta(T0);
+ params[2] = genrand_real2();
+
+ status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals);
+
+ if (status) {
+ sprintf(err, "error integrating electron angular distribution: %s\n", gsl_strerror(status));
+ goto err;
+ }
+
+ if (!isclose(result, 1.0, 0, 1e-3)) {
+ sprintf(err, "integral of electron_get_angular_pdf() returned %.5g, but expected %.5g", result, 1.0);
+ goto err;
+ }
+ }
+
+ gsl_integration_cquad_workspace_free(w);
+
+ return 0;
+
+err:
+ gsl_integration_cquad_workspace_free(w);
+
+ return 1;
+}
+
+int test_gamma_pdf(char *err)
+{
+ /* Tests the gamma_pdf() function. */
+ size_t i;
+ double x, k, theta, expected;
+ double result;
+
+ init_genrand(0);
+
+ for (i = 0; i < 100; i++) {
+ x = genrand_real2();
+ k = genrand_real2();
+ theta = genrand_real2();
+
+ result = gamma_pdf(x,k,theta);
+
+ expected = gsl_ran_gamma_pdf(x,k,theta);
+
+ if (!isclose(result, expected, 0, 1e-9)) {
+ sprintf(err, "gamma_pdf() returned %.5g, but expected %.5g", result, expected);
+ goto err;
+ }
+ }
+
+ return 0;
+
+err:
+ return 1;
+}
+
+int test_solid_angle2(char *err)
+{
+ /* Tests the get_solid_angle()2 function. */
+ size_t i;
+ double pmt[3] = {0,0,0};
+ double pos[3] = {0,0,1};
+ double n[3] = {0,0,1};
+ double r = 1.0;
+ double solid_angle;
+ double dir[3];
+ double L, r0, R;
+ double expected;
+
+ init_genrand(0);
+
+ for (i = 0; i < 100; i++) {
+ rand_sphere(pmt);
+ rand_sphere(n);
+ r = genrand_real2();
+
+ dir[0] = pos[0] - pmt[0];
+ dir[1] = pos[1] - pmt[1];
+ dir[2] = pos[2] - pmt[2];
+
+ L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+ R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+ r0 = sqrt(R*R - L*L);
+
+ solid_angle = get_solid_angle2(L/r,r0/r);
+
+ expected = get_solid_angle(pos,pmt,n,r);
+
+ if (!isclose(solid_angle, expected, 1e-4, 0)) {
+ sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
+int test_solid_angle_lookup(char *err)
+{
+ /* Tests the get_solid_angle_lookup() function. */
+ size_t i;
+ double pmt[3] = {0,0,0};
+ double pos[3] = {0,0,1};
+ double n[3] = {0,0,1};
+ double r = 1.0;
+ double solid_angle;
+ double dir[3];
+ double L, r0, R;
+ double expected;
+
+ init_genrand(0);
+
+ for (i = 0; i < 100; i++) {
+ rand_sphere(pmt);
+ MUL(pmt,800);
+ rand_sphere(n);
+ r = genrand_real2();
+
+ dir[0] = pos[0] - pmt[0];
+ dir[1] = pos[1] - pmt[1];
+ dir[2] = pos[2] - pmt[2];
+
+ L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+ R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+ r0 = sqrt(R*R - L*L);
+
+ solid_angle = get_solid_angle_lookup(L/r,r0/r);
+
+ expected = get_solid_angle(pos,pmt,n,r);
+
+ if (!isclose(solid_angle, expected, 1e-2, 0)) {
+ sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
+int test_solid_angle_fast(char *err)
+{
+ /* Tests the get_solid_angle_lookup() function. */
+ size_t i;
+ double pmt[3] = {0,0,0};
+ double pos[3] = {0,0,1};
+ double n[3] = {0,0,1};
+ double r = 1.0;
+ double solid_angle;
+ double dir[3];
+ double L, r0, R;
+ double expected;
+
+ init_genrand(0);
+
+ for (i = 0; i < 100; i++) {
+ rand_sphere(pmt);
+ MUL(pmt,800);
+ rand_sphere(n);
+ r = genrand_real2();
+
+ dir[0] = pos[0] - pmt[0];
+ dir[1] = pos[1] - pmt[1];
+ dir[2] = pos[2] - pmt[2];
+
+ L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]);
+ R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
+ r0 = sqrt(R*R - L*L);
+
+ solid_angle = get_solid_angle_fast(pos,pmt,n,r);
+
+ expected = get_solid_angle(pos,pmt,n,r);
+
+ if (!isclose(solid_angle, expected, 1e-2, 0)) {
+ sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, expected);
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
struct tests {
testFunction *test;
char *name;
@@ -999,6 +1213,11 @@ struct tests {
{test_log_norm, "test_log_norm"},
{test_trapz, "test_trapz"},
{test_interp2d, "test_interp2d"},
+ {test_electron_get_angular_pdf, "test_electron_get_angular_pdf"},
+ {test_gamma_pdf, "test_gamma_pdf"},
+ {test_solid_angle2, "test_solid_angle2"},
+ {test_solid_angle_lookup, "test_solid_angle_lookup"},
+ {test_solid_angle_fast, "test_solid_angle_fast"},
};
int main(int argc, char **argv)
generated by cgit (git 2.34.1) at 2025-05-07 00:09:13 +0000