diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 622 |
1 files changed, 392 insertions, 230 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 6f79325..34f379a 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -39,6 +39,8 @@ #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) { @@ -58,7 +60,7 @@ particle *particle_init(int id, double T0, size_t n) * 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; + double dEdx, rad, norm, pdf, pdf_last; particle *p = malloc(sizeof(particle)); p->id = id; @@ -67,9 +69,18 @@ particle *particle_init(int id, double T0, size_t 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; @@ -86,6 +97,8 @@ particle *particle_init(int id, double T0, size_t n) 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 @@ -109,6 +122,34 @@ particle *particle_init(int id, double T0, size_t n) p->shower_photons = electron_get_shower_photons(T0, rad); + norm = electron_get_angular_pdf_norm(p->a, p->b, 1/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. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm; + 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); + } + break; case IDP_MU_MINUS: case IDP_MU_PLUS: @@ -123,6 +164,8 @@ particle *particle_init(int id, double T0, size_t n) 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); @@ -145,7 +188,29 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; - p->shower_photons = rad*MUON_PHOTONS_PER_MEV; + /* FIXME: add shower photons for muons. */ + p->shower_photons = 0.0; + + /* Calculate the normalization constant for the angular distribution. */ + norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/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. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm; + 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); break; case IDP_PROTON: @@ -154,6 +219,7 @@ particle *particle_init(int id, double T0, size_t n) * 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 @@ -175,6 +241,7 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + /* FIXME: Is this correct? */ p->shower_photons = rad*PROTON_PHOTONS_PER_MEV; break; @@ -186,77 +253,243 @@ particle *particle_init(int id, double T0, size_t n) return p; } -double particle_get_energy(double x, particle *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) { - /* 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; + return (R*cos_gamma-x)/sqrt(R*R+x*x-2*x*R*cos_gamma); +} - T = interp1d(x,p->x,p->T,p->n); +/* 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 (T < 0) return 0; + /* If cos(theta) == 1.0 we must be at the beginning of the track. */ + if (cos_theta == 1.0) return R*cos_gamma; - return T; -} + /* 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); -void particle_free(particle *p) -{ - free(p->x); - free(p->T); - free(p); + /* This result comes from using the law of sines. */ + return R*(cos_gamma - cos_theta*sin_gamma/sin_theta); } -static void get_expected_charge_delta_ray(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t) +/* 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) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o; + size_t i; + double x1, x2, cdf1, cdf2, cdf, delta, cos_theta, cos_theta_last, xcdf, xcdf_last; - SUB(pmt_dir,pmts[pmt].pos,pos); + /* First we need to figure out the bounds of integration in x. */ + x1 = 0.0; - normalize(pmt_dir); + 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++) { + cdf = cdf1 + i*delta; + cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta); + xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + 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; + } - cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal); + return 0; +} - *t = 0.0; - *q = 0.0; - *reflected = 0.0; - if (cos_theta_pmt > 0) return; +/* 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; - /* Calculate the cosine of the angle between the track direction and the - * vector to the PMT. */ - cos_theta = DOT(dir,pmt_dir); + /* First we need to figure out the bounds of integration in x. */ - omega = get_solid_angle_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS); + /* 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; - f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); - f = get_weighted_pmt_response(-cos_theta_pmt); + x1 = p->xlo_shower; - get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); + 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++) { + cdf = cdf1 + i*delta; + cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower); + xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + 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; + } - 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); + return 0; +} - constant = get_weighted_quantum_efficiency()*omega/(2*M_PI); +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; - /* 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. */ - charge = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o); + T = interp1d(x,p->x,p->T,p->n); - *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; + if (T < 0) return 0; - *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT; + return T; +} - *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge; +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], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o; + 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); @@ -264,15 +497,6 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal); - *t = 0.0; - *q = 0.0; - *reflected = 0.0; - if (cos_theta_pmt > 0) return; - - /* 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,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS); f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt); @@ -283,15 +507,7 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in 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); - constant = get_weighted_quantum_efficiency()*omega/(2*M_PI); - - /* 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. */ - charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o); + charge = get_weighted_quantum_efficiency()*omega/(2*M_PI); *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; @@ -300,9 +516,14 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge; } +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, n, omega, z, R, f, f_reflec, cos_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o; + double pmt_dir[3], cos_theta, sin_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o; z = 1.0; @@ -322,10 +543,14 @@ static void get_expected_charge(double beta, double theta0, double *pos, double /* Calculate the cosine of the angle between the track direction and the * vector to the PMT. */ cos_theta = DOT(dir,pmt_dir); - - if (fabs(cos_theta-1.0/(n*beta))/theta0 > 5) return; + 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(R,PMT_RADIUS,sin_theta_pmt)); + + if (fabs(cos_theta-1.0/(n*beta))/(sin_theta*theta0) > 5) return; if (cos_theta_pmt > 0) return; @@ -398,7 +623,7 @@ double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_ * 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 + * enough* for electrons too. The probability distribution for indirect * light is assumed to be a step function past some time `tmean`. * * The probability returned is calculated by taking the sum of these three @@ -437,66 +662,10 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m 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_delta_ray(particle *p, double *x, double *pdf, 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, dx, 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_delta_ray(p, pos, dir0, pmt, &q, &r, &t); - - q = q*pdf[i]; - r = r*pdf[i]; - t += x[i]/SPEED_OF_LIGHT; - - if (i == 0 || i == (n - 1)) { - q_sum += q; - r_sum += r; - t_sum += t*q; - t2_sum += t*t*q; - } else { - q_sum += 2*q; - r_sum += 2*r; - t_sum += 2*t*q; - t2_sum += 2*t*t*q; - } - } - - dx = x[1] - x[0]; - *mu_direct = q_sum*dx*0.5; - *mu_indirect = r_sum*dx*0.5; - - if (*mu_direct == 0.0) { - *time = 0.0; - *sigma = PMT_TTS; - } else { - *time = t_sum*dx*0.5/(*mu_direct); - - /* Variance in the time = E(t^2) - E(t)^2. */ - tmp = t2_sum*dx*0.5/(*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_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *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, dx, tmp, q_sum, r_sum, t_sum, t2_sum; + double pos[3], t, q, r, tmp, q_sum, r_sum, t_sum, t2_sum; q_sum = 0.0; r_sum = 0.0; @@ -509,35 +678,25 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &t); - q = q*pdf[i]; - r = r*pdf[i]; t += x[i]/SPEED_OF_LIGHT; - if (i == 0 || i == (n - 1)) { - q_sum += q; - r_sum += r; - t_sum += t*q; - t2_sum += t*t*q; - } else { - q_sum += 2*q; - r_sum += 2*r; - t_sum += 2*t*q; - t2_sum += 2*t*t*q; - } + q_sum += q*w[i]; + r_sum += r*w[i]; + t_sum += t*q*w[i]; + t2_sum += t*t*q*w[i]; } - dx = x[1]-x[0]; - *mu_direct = q_sum*dx*0.5; - *mu_indirect = r_sum*dx*0.5; + *mu_direct = q_sum; + *mu_indirect = r_sum; if (*mu_direct == 0.0) { *time = 0.0; *sigma = PMT_TTS; } else { - *time = t_sum*dx*0.5/(*mu_direct); + *time = t_sum/(*mu_direct); /* Variance in the time = E(t^2) - E(t)^2. */ - tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time); + tmp = t2_sum/(*mu_direct) - (*time)*(*time); if (tmp >= 0) { *sigma = sqrt(tmp + PMT_TTS*PMT_TTS); @@ -757,7 +916,7 @@ 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; + double T, E, mom, beta, tmp; betaRootParams *pars; pars = (betaRootParams *) params; @@ -766,8 +925,15 @@ static double beta_root(double x, void *params) /* Calculate total energy */ E = T + pars->p->mass; - mom = sqrt(E*E - pars->p->mass*pars->p->mass); - beta = mom/E; + + 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; } @@ -801,7 +967,7 @@ static int get_smax(particle *p, double beta_min, double range, double *smax) 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-1, 0); + status = gsl_root_test_interval(x_lo, x_hi, 1e-10, 0); if (status == GSL_SUCCESS) break; } while (status == GSL_CONTINUE && iter < max_iter); @@ -964,19 +1130,23 @@ double nll_best(event *ev) return kahan_sum(nll,nhit); } -double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const int fast, int charge_only, int hit_only) +/* 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) { - /* 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. - * - * 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. */ size_t i, j, k, nhit; static double logp[MAX_PE], nll[MAX_PMTS]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; @@ -1007,26 +1177,38 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in static double mu_sum[MAX_PMTS]; double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total; - static double x[MAX_NPOINTS], pdf[MAX_NPOINTS]; - static double x_delta[MAX_NPOINTS], pdf_delta[MAX_NPOINTS]; + /* 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; - double pos_a, pos_b; + size_t npoints; - double xlo, xhi; + double log_pt_fast; - size_t npoints, npoints_shower; + double distance_to_psup; - double log_pt_fast; + 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); + } + /* Initialize static arrays. */ for (i = 0; i < MAX_PMTS; i++) { mu_sum[i] = 0.0; @@ -1040,6 +1222,18 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in } 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 (!intersect_sphere(v[j].pos,v[j].dir,PSUP_RADIUS,&distance_to_psup)) { + /* If the particle is outside 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 @@ -1047,12 +1241,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in range = p->range; - /* Number of points to sample along the particle track. */ - npoints = (size_t) (range/dx + 0.5); - - if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS; - if (npoints > MAX_NPOINTS) npoints = MAX_NPOINTS; - /* Calculate total energy */ E0 = v[j].T0 + p->mass; p0 = sqrt(E0*E0 - p->mass*p->mass); @@ -1061,63 +1249,23 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in /* FIXME: is this formula valid for muons? */ theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2); - if (!fast) { - path = path_init(v[j].pos, v[j].dir, v[j].T0, range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass); - - switch (v[j].id) { - case IDP_E_MINUS: - case IDP_E_PLUS: - electron_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); - break; - case IDP_MU_MINUS: - case IDP_MU_PLUS: - muon_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); - break; - default: - fprintf(stderr, "unknown particle id %i\n", v[j].id); - exit(1); - } - - /* Lower and upper limits to integrate for the electromagnetic - * shower. */ - xlo = 0.0; - xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); - - /* Number of points to sample along the longitudinal direction for - * the electromagnetic shower. */ - npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5); - - if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS; - if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS; + /* Compute the distance at which the particle falls below the cerenkov + * threshold. */ + if (beta0 > 1/avg_index_d2o) + get_smax(p, 1/avg_index_d2o, range, &cerenkov_range); + else + cerenkov_range = range; - for (i = 0; i < npoints_shower; i++) { - x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); - pdf[i] = gamma_pdf(x[i],pos_a,pos_b); - } + if (distance_to_psup < cerenkov_range) + cerenkov_range = distance_to_psup; - /* Make sure the PDF is normalized. */ - double norm = trapz(pdf,x[1]-x[0],npoints_shower); - for (i = 0; i < npoints_shower; i++) { - pdf[i] /= norm; - } + /* Number of points to sample along the particle track. */ + npoints = (size_t) (cerenkov_range/dx + 0.5); - for (i = 0; i < npoints_shower; i++) { - x_delta[i] = i*p->range/(npoints_shower-1); - /* For now, we assume the delta rays are produced uniformly - * along the path. This isn't quite correct, but the - * distribution is close to flat. - * - * FIXME: Figure out a good parameterization for the - * distribution of delta ray photons along the path. */ - pdf_delta[i] = 1.0; - } + if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS; - /* Make sure the PDF is normalized. */ - norm = trapz(pdf_delta,x_delta[1]-x_delta[0],npoints_shower); - for (i = 0; i < npoints_shower; i++) { - pdf_delta[i] /= norm; - } - } + 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); @@ -1133,6 +1281,9 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in 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 @@ -1146,7 +1297,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in continue; } - double q_indirect; integrate_path(path, i, &mu[i][j][0], &q_indirect, &result); mu_indirect[j] += q_indirect; @@ -1166,14 +1316,26 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in 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) { - integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]); + 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) { - integrate_path_delta_ray(p,x_delta,pdf_delta,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]); + 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; } |