sddm - Code for the self-destructing dark matter search in SNO

diff options

Diffstat (limited to 'src')

0 files changed, 0 insertions, 0 deletions

/* 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 "likelihood.h" #include <stdlib.h> /* for size_t */ #include "pmt.h" #include <gsl/gsl_integration.h> #include "muon.h" #include "misc.h" #include <gsl/gsl_sf_gamma.h> #include "sno.h" #include "vector.h" #include "event.h" #include "optics.h" #include "sno_charge.h" #include "pdg.h" #include "path.h" #include <stddef.h> /* for size_t */ #include "scattering.h" #include "solid_angle.h" #include <gsl/gsl_roots.h> #include <gsl/gsl_errno.h> #include "pmt_response.h" #include "electron.h" #include "proton.h" #include "id_particles.h" #include <gsl/gsl_cdf.h> #include "find_peaks.h" #include <gsl/gsl_interp.h> #include <gsl/gsl_spline.h> particle *particle_init(int id, double T0, size_t n) { /* Returns a struct with arrays of the particle position and kinetic * energy. This struct can then be passed to particle_get_energy() to * interpolate the particle's kinetic energy at any point along the track. * For example: * * particle *p = particle_init(IDP_MU_MINUS, 1000.0, 100); * double T = particle_get_energy(x, p); * * To compute the kinetic energy as a function of distance we need to solve * the following integral equation: * * T(x) = T0 - \int_0^x dT(T(x))/dx * * which we solve by dividing the track up into `n` segments and then * numerically computing the energy at each step. */ size_t i; double dEdx, rad, norm, pdf, pdf_last; particle *p = malloc(sizeof(particle)); p->id = id; p->x = malloc(sizeof(double)*n); p->T = malloc(sizeof(double)*n); p->n = n; p->a = 0.0; p->b = 0.0; p->cos_theta = calloc(n, sizeof(double)); p->cdf_shower = calloc(n, sizeof(double)); p->spline_shower = gsl_spline_alloc(gsl_interp_linear, n); p->acc_shower = gsl_interp_accel_alloc(); p->shower_photons = 0.0; p->delta_ray_a = 0.0; p->delta_ray_b = 0.0; p->cdf_delta = calloc(n, sizeof(double)); p->x_shower = calloc(n, sizeof(double)); p->gamma_pdf = calloc(n, sizeof(double)); p->spline_delta = gsl_spline_alloc(gsl_interp_linear, n); p->acc_delta = gsl_interp_accel_alloc(); p->delta_ray_photons = 0.0; p->x[0] = 0; p->T[0] = T0; switch (id) { case IDP_E_MINUS: case IDP_E_PLUS: p->mass = ELECTRON_MASS; /* We use the density of light water here since we don't have the * tables for heavy water for electrons. */ p->range = electron_get_range(T0, WATER_DENSITY); p->a = electron_get_angular_distribution_alpha(T0); p->b = electron_get_angular_distribution_beta(T0); electron_get_position_distribution_parameters(T0, &p->pos_a, &p->pos_b); p->delta_ray_photons = 0.0; /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ rad = 0.0; for (i = 1; i < n; i++) { 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); 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 * the particle is equal to BETA_MIN, we guarantee that there is a * point along the track where the speed drops below BETA_MIN. * * A possible future improvement would be to dynamically compute the * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; p->shower_photons = electron_get_shower_photons(T0, rad); break; case IDP_MU_MINUS: case IDP_MU_PLUS: p->mass = MUON_MASS; /* We use the density of heavy water here since we only load the heavy * water table for muons. * * FIXME: It would be nice in the future to load both tables and then * load the correct table based on the media. */ p->range = muon_get_range(T0, HEAVY_WATER_DENSITY); p->a = muon_get_angular_distribution_alpha(T0); p->b = muon_get_angular_distribution_beta(T0); muon_get_position_distribution_parameters(T0, &p->pos_a, &p->pos_b); muon_get_delta_ray_distribution_parameters(T0, &p->delta_ray_a, &p->delta_ray_b); p->delta_ray_photons = muon_get_delta_ray_photons(T0); /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ rad = 0.0; for (i = 1; i < n; i++) { 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], HEAVY_WATER_DENSITY); 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 * the particle is equal to BETA_MIN, we guarantee that there is a * point along the track where the speed drops below BETA_MIN. * * A possible future improvement would be to dynamically compute the * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; p->shower_photons = muon_get_shower_photons(T0, rad); break; case IDP_PROTON: p->mass = PROTON_MASS; /* We use the density of light water here since we don't have the * tables for heavy water for protons. */ p->range = proton_get_range(T0, WATER_DENSITY); /* FIXME: add delta ray photons for muons. */ p->delta_ray_photons = 0.0; /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ rad = 0.0; for (i = 1; i < n; i++) { 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); 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 * the particle is equal to BETA_MIN, we guarantee that there is a * point along the track where the speed drops below BETA_MIN. * * A possible future improvement would be to dynamically compute the * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; /* FIXME: Add shower photons for protons. */ p->shower_photons = 0.0; break; default: fprintf(stderr, "unknown particle id %i\n", id); exit(1); } if (p->shower_photons) { norm = electron_get_angular_pdf_norm(p->a, p->b, 1/get_avg_index_d2o()); for (i = 0; i < n; i++) { p->cos_theta[i] = -1.0 + 2.0*i/(n-1); /* Note: We assume here that the peak of the angular distribution * is at the Cerenkov angle for a particle with beta = 1. This * seems to be an OK assumption for high energy showers, but is not * exactly correct for showers from lower energy electrons or for * delta rays from lower energy muons. It seems good enough, but in * the future it would be nice to parameterize this. * * Note: We add EPSILON to the pdf since the cdf values are * required to be strictly increasing in order to use gsl splines. */ pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/get_avg_index_d2o())/norm + EPSILON; if (i > 0) p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; else p->cdf_shower[i] = 0.0; pdf_last = pdf; } gsl_spline_init(p->spline_shower, p->cdf_shower, p->cos_theta, n); p->xlo_shower = gsl_cdf_gamma_Pinv(0.01,p->pos_a,p->pos_b); p->xhi_shower = gsl_cdf_gamma_Pinv(0.99,p->pos_a,p->pos_b); for (i = 0; i < n; i++) { p->x_shower[i] = p->xlo_shower + (p->xhi_shower-p->xlo_shower)*i/(n-1); p->gamma_pdf[i] = gamma_pdf(p->x_shower[i],p->pos_a,p->pos_b); } } if (p->delta_ray_photons) { /* Calculate the normalization constant for the angular distribution. */ norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/get_avg_index_d2o()); for (i = 0; i < n; i++) { p->cos_theta[i] = -1.0 + 2.0*i/(n-1); /* Note: We assume here that the peak of the angular distribution * is at the Cerenkov angle for a particle with beta = 1. This * seems to be an OK assumption for high energy showers, but is not * exactly correct for showers from lower energy electrons or for * delta rays from lower energy muons. It seems good enough, but in * the future it would be nice to parameterize this. * * Note: We add EPSILON to the pdf since the cdf values are * required to be strictly increasing in order to use gsl splines. */ pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/get_avg_index_d2o())/norm + EPSILON; if (i > 0) p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; else p->cdf_delta[i] = 0.0; pdf_last = pdf; } gsl_spline_init(p->spline_delta, p->cdf_delta, p->cos_theta, n); } return p; } /* Returns the cosine of the angle between a particle travelling in a straight * line and a PMT after a distance `x` in cm. * * `R` is the distance between the particle and the PMT at the start of the * path, and `cos_gamma` is the cosine of the angle between the particle * direction and the PMT at the start of the path. */ double x_to_cos_theta(double x, double R, double cos_gamma) { return (R*cos_gamma-x)/sqrt(R*R+x*x-2*x*R*cos_gamma); } /* Returns the x position along a straight particle path given the cosine of * the angle between the particle direction and a PMT. * * `R` is the distance between the particle and the PMT at the start of the * path, `cos_gamma` is the cosine of the angle between the particle * direction and the PMT at the start of the path, and `sin_gamma` is the sine * of the angle (this is passed simply to avoid an extra call to sin()). */ double cos_theta_to_x(double cos_theta, double R, double cos_gamma, double sin_gamma) { double sin_theta; /* If cos(theta) == 1.0 we must be at the beginning of the track. */ if (cos_theta == 1.0) return R*cos_gamma; /* We don't care about the sign here because theta is in the range [0,pi] * and so sin(theta) is always positive. */ sin_theta = sqrt(1-cos_theta*cos_theta); /* This result comes from using the law of sines. */ return R*(cos_gamma - cos_theta*sin_gamma/sin_theta); } /* Returns the x points and weights for numerically integrating the delta ray * photons. * * The integral we want to solve here is: * * \int x n(x) f(x) g(x) * * where n(x) is the number of delta ray photons emitted per cm, f(x) is the * angular distribution of the delta ray photons, and g(x) contains all the * other factors like solid angle, absorption, quantum efficiency, etc. * * We approximate this integral as: * * \sum_i w_i g(x_i) * * where the weights and points are returned by this function. * * We currently assume the number of delta ray photons is constant along the * particle track, i.e. * * n(x) = 1/range * * This is an OK approximation, but in the future it might be nice to * parameterize the number of photons produced as a function of x. * * `distance_to_psup` is the distance to the PSUP. If the distance to the PSUP * is shorter than the range then we only integrate to the PSUP. * * Other arguments: * * R - distance to PMT from start of track * cos_gamma - cosine of angle between track and PMT at start of track * sin_gamma - sine of angle between track and PMT at start of track * * Returns 1 if the integration can be skipped, and 0 otherwise. */ int get_delta_ray_weights(particle *p, double range, double distance_to_psup, double *x, double *w, size_t n, double R, double cos_gamma, double sin_gamma) { size_t i; double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last; /* First we need to figure out the bounds of integration in x. */ x1 = 0.0; if (range < distance_to_psup) x2 = range; else x2 = distance_to_psup; /* Map bounds of integration in x -> cos(theta) and then find the values of * the CDF at the endpoints. */ cdf1 = interp1d(x_to_cos_theta(x1,R,cos_gamma),p->cos_theta,p->cdf_delta,p->n); cdf2 = interp1d(x_to_cos_theta(x2,R,cos_gamma),p->cos_theta,p->cdf_delta,p->n); delta = (cdf2-cdf1)/n; /* Now, we loop over different values of the CDF spaced evenly by delta and * compute cos(theta) for each of these values, map that value back to x, * and then compute the weight associated with that point. */ for (i = 0; i < n + 1; i++) { if (i < n) cdf = cdf1 + i*delta; else cdf = cdf2; cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); /* Occasionally due to floating point rounding error, it's possible * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, * we just make sure that xcdf is bounded by x1 and x2. */ if (xcdf < x1) xcdf = x1; if (xcdf > x2) xcdf = x2; if (i > 0) { /* For each interval we approximate the integral as: * * \int f(x) g(x) ~ <f(x)> \int g(x) * * and so we are free to choose best how to approximate the second * integral. Since g(x) is assumed to be slowly varying as a * function of x, it shouldn't matter too much, however based on * some testing I found that it was best to use the midpoint, i.e. * * \int g(x) ~ g((a+b)/2) * * as opposed to evaluating g(x) at the beginning, end, or * averaging the function at the two points. */ x[i-1] = (xcdf + xcdf_last)/2; w[i-1] = p->delta_ray_photons*delta*(1/range)*(xcdf - xcdf_last)/(cos_theta - cos_theta_last); } cos_theta_last = cos_theta; xcdf_last = xcdf; } return 0; } /* Returns the x points and weights for numerically integrating the shower * photons. * * The integral we want to solve here is: * * \int x n(x) f(x) g(x) * * where n(x) is the number of shower photons emitted per cm, f(x) is the * angular distribution of the shower photons, and g(x) contains all the other * factors like solid angle, absorption, quantum efficiency, etc. * * We approximate this integral as: * * \sum_i w_i g(x_i) * * where the weights and points are returned by this function. * * `distance_to_psup` is the distance to the PSUP. If the distance to the PSUP is * shorter than the shower range then we only integrate to the PSUP. * * Other arguments: * * R - distance to PMT from start of track * cos_gamma - cosine of angle between track and PMT at start of track * sin_gamma - sine of angle between track and PMT at start of track * * Returns 1 if the integration can be skipped, and 0 otherwise. */ int get_shower_weights(particle *p, double distance_to_psup, double *x, double *w, size_t n, double R, double cos_gamma, double sin_gamma) { size_t i; double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last; /* First we need to figure out the bounds of integration in x. */ /* If the start of the shower is past the PSUP, no light will reach the * PMTs, so we just return 1. */ if (p->xlo_shower > distance_to_psup) return 1; x1 = p->xlo_shower; if (p->xhi_shower > distance_to_psup) x2 = distance_to_psup; else x2 = p->xhi_shower; /* Map bounds of integration in x -> cos(theta) and then find the values of * the CDF at the endpoints. */ cdf1 = interp1d(x_to_cos_theta(x1,R,cos_gamma),p->cos_theta,p->cdf_shower,p->n); cdf2 = interp1d(x_to_cos_theta(x2,R,cos_gamma),p->cos_theta,p->cdf_shower,p->n); delta = (cdf2-cdf1)/n; /* Now, we loop over different values of the CDF spaced evenly by delta and * compute cos(theta) for each of these values, map that value back to x, * and then compute the weight associated with that point. */ for (i = 0; i < n + 1; i++) { if (i < n) cdf = cdf1 + i*delta; else cdf = cdf2; cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); /* Occasionally due to floating point rounding error, it's possible * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, * we just make sure that xcdf is bounded by x1 and x2. */ if (xcdf < x1) xcdf = x1; if (xcdf > x2) xcdf = x2; if (i > 0) { /* For each interval we approximate the integral as: * * \int f(x) g(x) ~ <f(x)> \int g(x) * * and so we are free to choose best how to approximate the second * integral. Since g(x) is assumed to be slowly varying as a * function of x, it shouldn't matter too much, however based on * some testing I found that it was best to use the midpoint, i.e. * * \int g(x) ~ g((a+b)/2) * * as opposed to evaluating g(x) at the beginning, end, or * averaging the function at the two points. */ x[i-1] = (xcdf + xcdf_last)/2; w[i-1] = p->shower_photons*delta*interp1d(x[i-1],p->x_shower,p->gamma_pdf,p->n)*(xcdf-xcdf_last)/(cos_theta - cos_theta_last); } cos_theta_last = cos_theta; xcdf_last = xcdf; } return 0; } double particle_get_energy(double x, particle *p) { /* Returns the approximate kinetic energy of a particle in water after * travelling `x` cm with an initial kinetic energy `T`. * * Return value is in MeV. */ double T; T = interp1d(x,p->x,p->T,p->n); if (T < 0) return 0; return T; } void particle_free(particle *p) { free(p->x); free(p->T); free(p->cos_theta); free(p->cdf_shower); free(p->x_shower); free(p->gamma_pdf); free(p->cdf_delta); gsl_spline_free(p->spline_shower); gsl_spline_free(p->spline_delta); gsl_interp_accel_free(p->acc_shower); gsl_interp_accel_free(p->acc_delta); free(p); } static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t) { double pmt_dir[3], omega, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o; SUB(pmt_dir,pmts[pmt].pos,pos); normalize(pmt_dir); cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal); omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS); f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); f = get_weighted_pmt_response(-cos_theta_pmt); get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS); charge = get_weighted_quantum_efficiency()*omega/(2*M_PI); /* Note: We multiply the amount of Rayleigh scattered light here by 2 since * we are only calculating the fraction of Rayleigh scattered light that * would hit the PMT and the coverage in SNO+ is about 50%. */ *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge; *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT; *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } /* Returns the angular width of the PMT bucket from the center of the PMT to * the edge of the bucket in the plane formed by the particle direction and the * direction to the PMT. * * `R` is the distance to the PMT, `r` is the PMT radius, and `sin_theta_pmt` * is the sine of the angle between the vector from the particle position to * the PMT and the PMT normal vector. * * This function is called get_theta0_min because we use this angular width to * set a minimum value for the RMS scattering angle of the particle as a kind * of hack to deal with the fact that we assume that the angular distribution * is constant across the face of the PMT. By introducing a minimum value for * the scattering RMS we broaden the angular distribution such that it * effectively averages across the face of a PMT. */ double get_theta0_min(double R, double r, double sin_theta_pmt) { return fast_acos((R-r*sin_theta_pmt)/sqrt(r*r + R*R - 2*r*R*sin_theta_pmt)); } static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t) { double pmt_dir[3], cos_theta, sin_theta, n, omega, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov, distance_to_pmt; /* Previously the index of refraction was calculated based on the position, * i.e.: * * n = (NORM(pos) <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o(); * * but this was causing a discontinuity in the likelihood function when the * particle track crossed the AV since we're numerically integrating over a * finite set of points. */ n = get_avg_index_d2o(); cos_theta_cerenkov = 1/(beta*n); *q = 0.0; *t = 0.0; *reflected = 0.0; if (cos_theta_cerenkov > 1) return; SUB(pmt_dir,pmts[pmt].pos,pos); distance_to_pmt = NORM(pmt_dir); normalize(pmt_dir); /* Calculate the cosine of the angle between the track direction and the * vector to the PMT. */ cos_theta = DOT(dir,pmt_dir); sin_theta = sqrt(1-cos_theta*cos_theta); cos_theta_pmt = -DOT(pmt_dir,pmts[pmt].normal); sin_theta_pmt = sqrt(1-cos_theta_pmt*cos_theta_pmt); theta0 = fmax(theta0,get_theta0_min(distance_to_pmt,PMT_RADIUS,sin_theta_pmt)); #ifdef FAST_GET_EXPECTED_CHARGE /* This next line is used to skip out of calculating the expected charge if * abs((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0)) > 5. The idea * here is that later in the likelihood calculation the PDF for Cerenkov * light has a term like * * exp(((cos(theta)-cos_theta_cherenkov)/(sin(theta)*theta0))**2) * * which will be less than 10^-25 if the following inequality holds and * therefore the charge should be negligible. * * However! I noticed that this line was causing discontinuities in the * likelihood when fitting low energy muons. I realized there are two * potential issues with this. One is that the PDF is multiplied by two * other terms: 1/(sin_theta*theta0) and the solid angle of the PMT. * Therefore if these are large it may have an impact. Second, the PDF * actually integrates over the Cerenkov spectrum and correctly takes into * account the change in the index as a function of wavelength. * * For now we don't use this to prevent the discontinuities. However, it * would be nice to use it in the future since it provides a massive ~2.5x * speedup. One idea for trying to eventually include it again is to update * the PDF for Cerenkov light to not include the wavelength dependence * (which would make it slightly less accurate, but it may be worth the * 2.5x speedup). */ if (fabs(cos_theta-cos_theta_cerenkov) > 5*sin_theta*theta0) return; #endif omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS); f_reflec = get_weighted_pmt_reflectivity(cos_theta_pmt); f = get_weighted_pmt_response(cos_theta_pmt); get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT; /* Probability that a photon is absorbed. We calculate this by computing: * * 1.0 - P(not absorbed in D2O)*P(not absorbed in H2O)*P(not absorbed in acrylic) * * since if we worked with the absorption probabilities directly it would * be more complicated, i.e. * * P(absorbed in D2O) + P(absorbed in acrylic|not absorbed in D2O)*P(not absorbed in D2O) + ... * */ prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); /* Similiar calculation for the probability that a photon is scattered. * * Technically we should compute this conditionally on the probability that * a photon is not absorbed, but since the probability of scattering is * pretty low, this is expected to be a very small effect. */ prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS); charge = omega*(1-cos_theta_cerenkov*cos_theta_cerenkov)*get_probability(beta, cos_theta, sin_theta, theta0); /* Note: We multiply the amount of Rayleigh scattered light here by 2 since * we are only calculating the fraction of Rayleigh scattered light that * would hit the PMT and the coverage in SNO+ is about 50%. */ *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge; *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } double time_cdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma) { /* Returns the CDF for the time distribution of photons at time `t`. */ size_t i; double p, mu_total; p = mu_noise*t/GTVALID; if (t < tmean) ; else if (t < tmean + 2*PSUP_REFLECTION_TIME) p += mu_indirect*(t-tmean)/(2*PSUP_REFLECTION_TIME); else p += mu_indirect; mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { if (mu[i] == 0.0) continue; p += mu[i]*norm_cdf(t,ts[i],ts_sigma[i]); mu_total += mu[i]; } return p/mu_total; } /* Returns the probability that a photon is detected at time `t`. * * The probability distribution is the sum of three different components: dark * noise, indirect light, and direct light. The dark noise is assumed to be * constant in time. The direct light is assumed to be a delta function around * the times `ts`, where each element of `ts` comes from a different particle. * This assumption is probably valid for particles like muons which don't * scatter much, and the hope is that it is *good enough* for electrons too. * The probability distribution for indirect light is assumed to be a step * function lasting 2*PSUP_REFLECTION_TIME past some time `tmean`. * * Note: Initially I had the indirect light be a step function which only * lasted PSUP_REFLECTION_TIME (which is set to the time it takes light to * cross from one side of the PSUP to the other) because I thought that this * represented the maximum amount of time it took for light to traverse the * detector (and the PMT hit time plot for a laserball run has a bump after the * main peak which lasts for about 80 ns). However, I later realized that for * events near the edge of the detector, like say a flasher or muon, the * reflected light from the other side of the detector will take twice as long. * So, now we use twice the PSUP reflection time. * * In the future it might be cool to experiment with new ways of defining the * time PDF. For example, we could get a much more accurate approximation for * the reflected light time distribution by using something like a kernel * density estimator and for each PMT we could loop over the other PMTs and add * a point to our PDF at the light travel time between the PMTs weighted by the * second PMT's charge. I think currently this would be way too slow, and it's * not obvious how much it would help. * * The probability returned is calculated by taking the sum of these three * components and convolving it with a gaussian with standard deviation `sigma` * which should typically be the PMT transit time spread. */ double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma) { size_t i; double p, mu_total; p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*PMT_TTS))-erf((t-tmean-2*PSUP_REFLECTION_TIME)/(sqrt(2)*PMT_TTS)))/(4*PSUP_REFLECTION_TIME); mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { if (mu[i] == 0.0) continue; p += mu[i]*norm(t,ts[i],ts_sigma[i]); mu_total += mu[i]; } return p/mu_total; } /* Returns the first order statistic for observing a PMT hit at time `t` given * `n` hits. * * The first order statistic is computed from the probability distribution * above. It's not obvious whether one should take the first order statistic * before or after convolving with the PMT transit time spread. Since at least * some of the transit time spread in SNO comes from the 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. * * Note: Technically we should include a final convolution representing both * the ECA+PCA uncertainty and the digitization noise. I was initially worried * about this causing an issue because as you take higher and higher order * statistics the width of the time PDF becomes smaller and smaller and these * other uncertainties may become relevant. However, the maximum number of PE * we should be calculating the time PDF for is when QLX rails which is at * approximately: * * (4096-600)*12/30 ~= 1000 * * and the width at that point is ~0.6 ns. The digitization noise is * approximately 4096/400.0 ~= 0.1 ns and according to Javi: * * > The ECA+PCA uncertainties are negligible with respect to the TTS. * > Basically, the width of the prompt peak for N16 calibration events (mainly * > SPEs) is 1.6ns, compatible with the PMT TTS. * > * > That said, long ago I estimated the precision of the ECA calibration to be * > about 0.4ns. The PCA delays are pretty stable and precisely measured, so * > probably negligible compared to the ECA uncertainties. * * So therefore, even at the highest possible charge we should still be * dominated by the width of the PDF. * * It is probably too expensive to try and actually calculate another * convolution at this point since we'd have to do it numerically. I think a * better option would just be to cut off the first order statistic at some * point, i.e. to calculate: * * logp[j] += log_pt(ev->pmt_hits[i].t, j > MAX_PE_TIME_PDF ? MAX_PE_TIME_PDF : j, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]); * */ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu, size_t n2, double *ts, double tmean, double *ts_sigma) { if (n == 1) return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)); else return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)); } static void integrate_path_shower(particle *p, double *x, double *w, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma) { size_t i; double pos[3], t, q, r, tmp, q_sum, r_sum, t_sum, t2_sum; q_sum = 0.0; r_sum = 0.0; t_sum = 0.0; t2_sum = 0.0; 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_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &t); t += x[i]/SPEED_OF_LIGHT; q_sum += q*w[i]; r_sum += r*w[i]; t_sum += t*q*w[i]; t2_sum += t*t*q*w[i]; } *mu_direct = q_sum; *mu_indirect = r_sum; if (*mu_direct == 0.0) { *time = 0.0; *sigma = PMT_TTS; } else { *time = t_sum/(*mu_direct); /* Variance in the time = E(t^2) - E(t)^2. */ tmp = t2_sum/(*mu_direct) - (*time)*(*time); if (tmp >= 0) { *sigma = sqrt(tmp + PMT_TTS*PMT_TTS); } else { /* This should never happen but does occasionally due to floating * point rounding error. */ *sigma = PMT_TTS; } } } static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indirect, double *time) { size_t i; double dir[3], pos[3], t, theta0, beta, q, r, q_sum, r_sum, t_sum, dx; q_sum = 0.0; r_sum = 0.0; t_sum = 0.0; for (i = 0; i < p->len; i++) { pos[0] = p->x[i]; pos[1] = p->y[i]; pos[2] = p->z[i]; dir[0] = p->dx[i]; dir[1] = p->dy[i]; dir[2] = p->dz[i]; beta = p->beta[i]; if (p->n > 0) { theta0 = p->theta0; } else { theta0 = fmax(MIN_THETA0,p->theta0*sqrt(p->s[i])); theta0 = fmin(MAX_THETA0,theta0); } get_expected_charge(beta, theta0, pos, dir, pmt, &q, &r, &t); t += p->t[i]; if (i == 0 || i == p->len - 1) { q_sum += q; r_sum += r; t_sum += t*q; } else { q_sum += 2*q; r_sum += 2*r; t_sum += 2*t*q; } } dx = p->s[1] - p->s[0]; *mu_direct = q_sum*dx*0.5; *mu_indirect = r_sum*dx*0.5; *time = t_sum*dx*0.5; } static double get_total_charge_approx(double beta0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double cos_theta_cerenkov, double sin_theta_cerenkov) { /* Returns the approximate expected number of photons seen by PMT `i` using * an analytic formula. * * To come up with an analytic formula for the expected number of photons, * it was necessary to make the following approximations: * * - the index of refraction is constant * - the particle track is a straight line * - the integral along the particle track is dominated by the gaussian * term describing the angular distribution of the light * * With these approximations and a few other ones (like using a Taylor * expansion for the distance to the PMT), it is possible to pull * everything out of the track integral and assume it's equal to it's value * along the track where the exponent of the gaussian dominates. * * The point along the track where the exponent dominates is calculated by * finding the point along the track where the angle between the track * direction and the PMT is equal to the Cerenkov angle. If this point is * before the start of the track, we use the start of the track and if it's * past the end of `smax` we use `smax`. * * Since the integral over the track also contains a term like * (1-1/(beta**2*n**2)) which is not constant near the end of the track, it * is necessary to define `smax` as the point along the track where the * particle velocity drops below some threshold. * * `smax` is currently calculated as the point where the particle velocity * drops to 0.8 times the speed of light. */ double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, f, cos_theta_pmt, l_h2o, l_d2o, f_reflec, charge, prob, frac, prob_abs, prob_sct; /* 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. */ if (sin_theta_cerenkov != 0) s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; else s = R*cos_theta; /* 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 vector from the point `s` along the track to the PMT. */ tmp[0] = pmts[i].pos[0] - (pos[0] + s*dir[0]); tmp[1] = pmts[i].pos[1] - (pos[1] + s*dir[1]); tmp[2] = pmts[i].pos[2] - (pos[2] + s*dir[2]); /* To do the integral analytically, we expand the distance to the PMT along * the track in a Taylor series around `s0`, i.e. * * r(s) = a + b*(s-s0) * * Here, we calculate `a` which is the distance to the PMT at the point * `s`. */ a = NORM(tmp); /* `z` is the distance to the PMT projected onto the track direction. */ z = R*cos_theta; /* `x` is the perpendicular distance from the PMT position to the track. */ x = R*sqrt(1-pow(cos_theta,2)); /* `b` is the second coefficient in the Taylor expansion. */ b = (s-z)/a; /* Compute the kinetic energy at the point `s` along the track. */ T = particle_get_energy(s,p); /* Calculate the particle velocity at the point `s`. */ E = T + p->mass; mom = sqrt(E*E - p->mass*p->mass); beta = mom/E; *t = 0.0; *mu_reflected = 0.0; if (beta < 1/get_avg_index_d2o()) return 0.0; /* `prob` is the number of photons emitted per cm by the particle at a * distance `s` along the track. */ prob = get_probability2(beta); /* 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); normalize(pmt_dir); cos_theta_pmt = DOT(pmt_dir,pmts[i].normal); *t = 0.0; *mu_reflected = 0.0; if (cos_theta_pmt > 0) return 0.0; /* Calculate the sine of the angle between the track direction and the PMT * at the position `s` along the track. */ cos_theta = DOT(dir,pmt_dir); 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_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); theta0 = fmax(theta0*sqrt(s),MIN_THETA0); frac = sqrt(2)*get_avg_index_d2o()*x*beta0*theta0; f = get_weighted_pmt_response(-cos_theta_pmt); f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS); prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS); /* Assume the particle is travelling at the speed of light. */ *t = s/SPEED_OF_LIGHT + l_d2o*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT; charge = get_avg_index_d2o()*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+get_avg_index_d2o()*(smax-z)*beta0)/frac) + erf((-a+b*s+get_avg_index_d2o()*z*beta0)/frac))/(b+get_avg_index_d2o()*beta0)/(4*M_PI); /* Add expected number of photons from electromagnetic shower. */ if (p->shower_photons > 0) charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/get_avg_index_d2o())*omega/(2*M_PI); if (p->delta_ray_photons > 0) charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/get_avg_index_d2o())*omega/(2*M_PI); /* Note: We multiply the amount of Rayleigh scattered light here by 2 since * we are only calculating the fraction of Rayleigh scattered light that * would hit the PMT and the coverage in SNO+ is about 50%. */ *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + 2*prob_sct*charge; return (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } typedef struct betaRootParams { particle *p; double beta_min; } betaRootParams; static double beta_root(double x, void *params) { /* Function used to find at what point along a track a particle crosses * some threshold in beta. */ double T, E, mom, beta, tmp; betaRootParams *pars; pars = (betaRootParams *) params; T = particle_get_energy(x, pars->p); /* Calculate total energy */ E = T + pars->p->mass; tmp = E*E - pars->p->mass*pars->p->mass; if (tmp <= 0) { beta = 0.0; } else { mom = sqrt(tmp); beta = mom/E; } return beta - pars->beta_min; } static int get_smax(particle *p, double beta_min, double range, double *smax) { /* Find the point along the track at which the particle's velocity drops to * `beta_min`. */ int status; betaRootParams pars; gsl_root_fsolver *s; gsl_function F; int iter = 0, max_iter = 100; double r, x_lo, x_hi; s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent); pars.p = p; pars.beta_min = beta_min; F.function = &beta_root; F.params = &pars; gsl_root_fsolver_set(s, &F, 0.0, range); do { iter++; status = gsl_root_fsolver_iterate(s); r = gsl_root_fsolver_root(s); x_lo = gsl_root_fsolver_x_lower(s); x_hi = gsl_root_fsolver_x_upper(s); /* Find the root to within 1 mm. */ status = gsl_root_test_interval(x_lo, x_hi, 1e-10, 0); if (status == GSL_SUCCESS) break; } while (status == GSL_CONTINUE && iter < max_iter); gsl_root_fsolver_free(s); *smax = r; return status; } static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, double cos_theta_cerenkov, double sin_theta_cerenkov) { /* Returns an approximate time at which a PMT is most likely to get hit * from Cerenkov light. * * To do this, we try to compute the distance along the track where the PMT * is at the Cerenkov angle. */ double pmt_dir[3], tmp[3]; double R, cos_theta, sin_theta, s, l_d2o, l_h2o; /* First, we find the point along the track where the PMT is at the * Cerenkov angle. */ SUB(pmt_dir,pmt_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. */ if (sin_theta_cerenkov != 0) s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; else s = R*cos_theta; /* 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 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,pmt_pos,tmp); get_path_length(tmp,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); /* Assume the particle is travelling at the speed of light. */ return s/SPEED_OF_LIGHT + l_d2o*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT; } static double getKineticEnergy(double x, void *p) { return particle_get_energy(x, (particle *) p); } double nll_best(event *ev) { /* Returns the negative log likelihood of the "best hypothesis" for the event `ev`. * * By "best hypothesis" I mean we restrict the model to assign a single * mean number of PE for each PMT in the event assuming that the number of * PE hitting each PMT is poisson distributed. In addition, the model picks * a single time for the light to arrive with the PDF being given by a sum * of dark noise and a Gaussian around this time with a width equal to the * single PE transit time spread. * * This calculation is intended to be used as a goodness of fit test for * the fitter by computing a likelihood ratio between the best fit * hypothesis and this ideal case. See Chapter 9 in Jayne's "Probability * Theory: The Logic of Science" for more details. */ size_t i, j, nhit; static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_sigma; double log_mu, max_logp, mu_noise, mu_indirect_total, min_ratio; mu_noise = DARK_RATE*GTVALID*1e-9; /* Compute the "best" number of expected PE for each PMT. */ for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; if (!ev->pmt_hits[i].hit) { /* If the PMT wasn't hit we assume the expected number of PE just * comes from noise hits. */ mu[i] = mu_noise; continue; } /* Set the expected number of PE to whatever maximizes P(q|mu), i.e. * * mu[i] = max(p(q|mu)) * */ mu[i] = get_most_likely_mean_pe(ev->pmt_hits[i].q); /* We want to estimate the mean time which is most likely to produce a * first order statistic of the actual hit time given there were * approximately mu[i] PE. As far as I know there are no closed form * solutions for the first order statistic of a Gaussian, but there is * an approximate form in the paper "Expected Normal Order Statistics * (Exact and Approximate)" by J. P. Royston which is what I use here. * * I found this via the stack exchange question here: * https://stats.stackexchange.com/questions/9001/approximate-order-statistics-for-normal-random-variables. */ double alpha = 0.375; ts[i] = ev->pmt_hits[i].t - gsl_cdf_ugaussian_Pinv((1-alpha)/(mu[i]-2*alpha+1))*PMT_TTS; } ts_sigma = PMT_TTS; mu_indirect_total = 0.0; min_ratio = MIN_RATIO; /* Now, we actually compute the negative log likelihood for the best * hypothesis. * * Currently this calculation is identical to the one in nll() so it should * probably be made into a function. */ nhit = 0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; log_mu = log(mu[i]); if (ev->pmt_hits[i].hit) { for (j = 1; j < MAX_PE; j++) { logp[j] = get_log_pq(ev->pmt_hits[i].q,j) + get_log_phit(j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], 1, &ts[i], ts[i], &ts_sigma); if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; if (logp[j] - max_logp < min_ratio*ln(10)) { j++; break; } } nll[nhit++] = -logsumexp(logp+1, j-1); } else { logp[0] = -mu[i]; for (j = 1; j < MAX_PE_NO_HIT; j++) { logp[j] = get_log_pmiss(j) - mu[i] + j*log_mu - lnfact(j); } nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT); } } return kahan_sum(nll,nhit); } /* Returns the negative log likelihood for event `ev` given a particle with id * `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and * initial time `t0`. * * `dx` is the distance the track is sampled at to compute the likelihood. * * `ns` is the number of points to use when calculating the number of photons * from shower and delta ray particles. * * If `fast` is nonzero, returns an approximate faster version of the * likelihood. * * `z1` and `z2` should be arrays of length `n` and represent coefficients in * the Karhunen Loeve expansion of the track path. These are only used when * `fast` is zero. */ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, int charge_only, int hit_only) { size_t i, j, k, nhit; static double logp[MAX_PE], nll[MAX_PMTS], tmp[2]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; double cos_theta_cerenkov, sin_theta_cerenkov; double logp_path; double mu_reflected; path *path; /* Array representing the expected number of photons at each PMT for each * vertex. The last axis represents the contributions from direct, shower, * and delta ray photons respectively, i.e. * * mu[0][2][0] - expected number of direct photons from vertex 2 at PMT 0 * mu[0][2][1] - expected number of shower photons from vertex 2 at PMT 0 * mu[0][2][2] - expected number of delta ray photons from vertex 2 at PMT 0 * ... * etc. * */ static double mu[MAX_PMTS][MAX_VERTICES][3]; /* Array representing the average time of arrival at each PMT for each * vertex. See above for the how the array is indexed. */ static double ts[MAX_PMTS][MAX_VERTICES][3]; /* Array representing the standard deviation of the time at each PMT for * each vertex. See above for the how the array is indexed. */ static double ts_sigma[MAX_PMTS][MAX_VERTICES][3]; static double mu_sum[MAX_PMTS]; double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total; /* Points and weights to integrate the shower and delta ray photon * distributions. */ static double x[MAX_NPOINTS], w[MAX_NPOINTS]; double result; double min_ratio; size_t npoints; double log_pt_fast; double distance_to_psup; double q_indirect; double R, cos_gamma, sin_gamma; double pmt_dir[3]; double cerenkov_range; if (n > MAX_VERTICES) { fprintf(stderr, "maximum number of vertices is %i\n", MAX_VERTICES); exit(1); } if (ns > MAX_NPOINTS) { fprintf(stderr, "maximum number of points is %i\n", MAX_NPOINTS); exit(1); } for (i = 0; i < LEN(mu_indirect); i++) { mu_indirect[i] = 0.0; } /* Initialize static arrays. */ for (i = 0; i < MAX_PMTS; i++) { mu_sum[i] = 0.0; for (j = 0; j < n; j++) { for (k = 0; k < 3; k++) { mu[i][j][k] = 0.0; ts[i][j][k] = 0.0; ts_sigma[i][j][k] = PMT_TTS; } } } for (j = 0; j < n; j++) { /* Find the distance from the particle's starting position to the PSUP. * * We calculate this here to make sure that we don't integrate outside * of the PSUP because otherwise there is no code to model the fact * that the PSUP is optically separated from the water outside. */ if (NORM(v[j].pos) > PSUP_RADIUS || !intersect_sphere(v[j].pos,v[j].dir,PSUP_RADIUS,&distance_to_psup)) { /* If the particle is outside the PSUP or it doesn't intersect the * PSUP we just assume it produces no light. */ fprintf(stderr, "particle doesn't intersect PSUP!\n"); continue; } if (fast) p = particle_init(v[j].id, v[j].T0, 1000); else p = particle_init(v[j].id, v[j].T0, 10000); range = p->range; /* Calculate total energy */ E0 = v[j].T0 + p->mass; p0 = sqrt(E0*E0 - p->mass*p->mass); beta0 = p0/E0; /* FIXME: is this formula valid for muons? */ theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2); /* Compute the distance at which the particle falls below the cerenkov * threshold. */ if (beta0 > 1/get_avg_index_d2o()) get_smax(p, 1/get_avg_index_d2o(), range, &cerenkov_range); else cerenkov_range = range; if (distance_to_psup < cerenkov_range) cerenkov_range = distance_to_psup; /* Number of points to sample along the particle track. */ npoints = (size_t) (cerenkov_range/dx + 0.5); if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS; if (!fast) path = path_init(v[j].pos, v[j].dir, v[j].T0, cerenkov_range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass); if (beta0 > BETA_MIN) get_smax(p, BETA_MIN, range, &smax); else smax = 0.0; /* Compute the Cerenkov angle at the start of the track. * * These are computed at the beginning of this function and then passed to * the different functions to avoid recomputing them on the fly. */ cos_theta_cerenkov = 1/(get_avg_index_d2o()*beta0); sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2)); mu_indirect[j] = 0.0; for (i = 0; i < MAX_PMTS; i++) { /* Note: We loop over all normal PMTs here, even ones that are * flagged since they will contribute to the reflected and * scattered light. */ if (pmts[i].pmt_type != PMT_NORMAL) continue; /* Skip PMTs which weren't hit when doing the "fast" likelihood * calculation. */ if (hit_only && !ev->pmt_hits[i].hit) continue; if (fast) { mu[i][j][0] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j][0], &mu_reflected, cos_theta_cerenkov, sin_theta_cerenkov); ts[i][j][0] += v[j].t0; mu_indirect[j] += mu_reflected; continue; } integrate_path(path, i, &mu[i][j][0], &q_indirect, &result); mu_indirect[j] += q_indirect; if (mu[i][j][0] > 1e-9) { ts[i][j][0] = v[j].t0 + result/mu[i][j][0]; } 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. * * Note: I don't really know how much this affects the likelihood. * I should test it to see if we can get away with something even * simpler (like just computing the distance from the start of the * track to the PMT). */ ts[i][j][0] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,cos_theta_cerenkov,sin_theta_cerenkov); } /* Calculate the distance and angle to the PMT from the starting * position. These are used to convert between the particle's x * position along the straight track and the cosine of the angle * with the PMT. */ SUB(pmt_dir,pmts[i].pos,v[j].pos); R = NORM(pmt_dir); cos_gamma = DOT(v[j].dir,pmt_dir)/R; sin_gamma = sqrt(1-cos_gamma*cos_gamma); if (p->shower_photons > 0) { get_shower_weights(p,distance_to_psup,x,w,ns,R,cos_gamma,sin_gamma); integrate_path_shower(p,x,w,v[j].T0,v[j].pos,v[j].dir,i,ns,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]); mu_indirect[j] += q_indirect; ts[i][j][1] = v[j].t0 + result; } if (p->delta_ray_photons > 0 && !get_delta_ray_weights(p,range,distance_to_psup,x,w,ns,R,cos_gamma,sin_gamma)) { integrate_path_shower(p,x,w,v[j].T0,v[j].pos,v[j].dir,i,ns,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]); mu_indirect[j] += q_indirect; ts[i][j][2] = v[j].t0 + result; } } if (!fast) { path_free(path); } particle_free(p); } mu_noise = DARK_RATE*GTVALID*1e-9; mu_indirect_total = 0.0; for (j = 0; j < n; j++) { mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION/10000.0; } /* Compute the expected number of photons reaching each PMT by adding up * the contributions from the noise hits and the direct, indirect, and * shower light. */ for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; /* Skip PMTs which weren't hit when doing the "fast" likelihood * calculation. */ if (hit_only && !ev->pmt_hits[i].hit) continue; mu_sum[i] = mu_indirect_total + mu_noise; for (j = 0; j < n; j++) { for (k = 0; k < 3; k++) { mu_sum[i] += mu[i][j][k]; } } } if (fast) min_ratio = MIN_RATIO_FAST; else min_ratio = MIN_RATIO; /* Now, we actually compute the negative log likelihood. */ nhit = 0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; /* Skip PMTs which weren't hit when doing the "fast" likelihood * calculation. */ if (hit_only && !ev->pmt_hits[i].hit) continue; log_mu = log(mu_sum[i]); if (fast) log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]); if (ev->pmt_hits[i].hit) { for (j = 1; j < MAX_PE; j++) { logp[j] = get_log_pq(ev->pmt_hits[i].q,j) + get_log_phit(j) - mu_sum[i] + j*log_mu - lnfact(j); if (!charge_only) { if (fast) logp[j] += log_pt_fast; else logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]); } if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; if (logp[j] - max_logp < min_ratio*ln(10)) { j++; break; } } /* Now, we include the probability that the channel is miscalibrated as: * * P(q,t) = P(q,t|calibrated)P(calibrated) + P(q,t|miscalibrated)P(miscalibrated) * * For P(q,t|miscalibrated) we just assume a uniform distribution * over all 4096 possible values. */ tmp[0] = logsumexp(logp+1, j-1) + log(1-P_MISCALIBRATION); tmp[1] = log(1/4096.0) + log(1/4096.0) + log(P_MISCALIBRATION); nll[nhit++] = -logsumexp(tmp, 2); } else { logp[0] = -mu_sum[i]; if (fast) { nll[nhit++] = -logp[0]; continue; } for (j = 1; j < MAX_PE_NO_HIT; j++) { logp[j] = get_log_pmiss(j) - mu_sum[i] + j*log_mu - lnfact(j); } /* Now, we include the probability that the channel is miscalibrated as: * * P(not hit) = P(not hit|calibrated)P(calibrated) + P(not hit|miscalibrated)P(miscalibrated) * * where by miscalibrated here we really mean a PMT that wasn't * read out or included in the event for whatever reason. */ tmp[0] = logsumexp(logp, MAX_PE_NO_HIT) + log(1-P_MISCALIBRATION); tmp[1] = log(P_MISCALIBRATION); nll[nhit++] = -logsumexp(tmp, 2); } } logp_path = 0.0; for (j = 0; j < n; j++) for (i = 0; i < v[j].n; i++) logp_path += log_norm(v[j].z1[i],0,1) + log_norm(v[j].z2[i],0,1); /* Add up the negative log likelihood terms from each PMT. I use a kahan * sum here instead of just summing all the values to help prevent round * off error. I actually don't think this makes any difference here since * the numbers are all of the same order of magnitude, so I should actually * test to see if it makes any difference. */ return kahan_sum(nll,nhit) - logp_path; }


context:
space:
mode: