#include "likelihood.h" #include /* for size_t */ #include "pmt.h" #include #include "muon.h" #include "misc.h" #include #include "sno.h" #include "vector.h" #include "event.h" #include "optics.h" #include "sno_charge.h" #include "pdg.h" #include "path.h" #include /* for size_t */ #include "scattering.h" #include "solid_angle.h" #include #include #include "pmt_response.h" #include "electron.h" #include "proton.h" #include "id_particles.h" #include #include "find_peaks.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; 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->shower_photons = 0.0; p->delta_ray_a = 0.0; p->delta_ray_b = 0.0; 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); 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 = rad*ELECTRON_PHOTONS_PER_MEV; 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_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], 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 = rad*MUON_PHOTONS_PER_MEV; 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); 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; p->shower_photons = rad*PROTON_PHOTONS_PER_MEV; break; default: fprintf(stderr, "unknown particle id %i\n", id); exit(1); } return p; } 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); } static double get_expected_charge_shower(particle *p, 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 *q_delta_ray, double *q_indirect_delta_ray) { double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant; SUB(pmt_dir,pmt_pos,pos); normalize(pmt_dir); cos_theta_pmt = DOT(pmt_dir,pmt_normal); charge = 0.0; *q_delta_ray = 0.0; *q_indirect_delta_ray = 0.0; *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); 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(); scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*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. */ if (p->shower_photons > 0) charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o); if (p->delta_ray_photons > 0) *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o); scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); *reflected = (f_reflec + 1.0 - scatter)*charge; *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray); *q_delta_ray *= f*scatter; return f*charge*scatter; } 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, scattering_length_h2o, scattering_length_d2o, scatter; z = 1.0; SUB(pmt_dir,pmt_pos,pos); 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); *reflected = 0.0; if (fabs(cos_theta-1.0/(n_d2o*beta))/theta0 > 5) return 0.0; cos_theta_pmt = DOT(pmt_dir,pmt_normal); if (cos_theta_pmt > 0) return 0.0; omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); R = NORM(pos); n = (R <= AV_RADIUS) ? n_d2o : n_h2o; *reflected = 0.0; if (beta < 1/n) return 0.0; 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(); scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); 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)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); *reflected = (f_reflec + 1.0 - scatter)*charge; return f*charge*scatter; } double time_cdf(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; double p, mu_total; p = mu_noise*t/GTVALID; if (t < tmean) ; else if (t < tmean + PSUP_REFLECTION_TIME) p += mu_indirect*(t-tmean)/PSUP_REFLECTION_TIME; else p += mu_indirect; mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { if (mu_direct[i] == 0.0) continue; p += mu_direct[i]*norm_cdf(t,ts[i],sigma); mu_total += mu_direct[i]; } for (i = 0; i < n; i++) { if (mu_shower[i] == 0.0) continue; p += mu_shower[i]*norm_cdf(t,ts_shower[i],ts_sigma[i]); mu_total += mu_shower[i]; } return p/mu_total; } double time_pdf(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`. * * 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 past some time `tmean`. * * 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. */ size_t i; double p, mu_total; p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*sigma))-erf((t-tmean-PSUP_REFLECTION_TIME)/(sqrt(2)*sigma)))/(2*PSUP_REFLECTION_TIME); mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { if (mu_direct[i] == 0.0) continue; p += mu_direct[i]*norm(t,ts[i],sigma); mu_total += mu_direct[i]; } for (i = 0; i < n; i++) { if (mu_shower[i] == 0.0) continue; 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, 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. * * 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. */ if (n == 1) return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); else return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); } 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) { 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, q_delta_ray, q_indirect_delta_ray; 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; 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(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray); qs[i] = qs[i]*pdf[i] + q_delta_ray/p->range; q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/p->range; ts[i] = t*qs[i]; ts2[i] = t*t*qs[i]; } double dx = x[1]-x[0]; *mu_direct = trapz(qs,dx,n); *mu_indirect = trapz(q_indirects,dx,n); *time = trapz(ts,dx,n)/(*mu_direct); double t2 = trapz(ts2,dx,n)/(*mu_direct); *sigma = sqrt(t2-(*time)*(*time)); if (*mu_direct == 0.0 || *sigma != *sigma) { *time = 0.0; *sigma = PMT_TTS; } else { *sigma = sqrt((*sigma)*(*sigma) + PMT_TTS*PMT_TTS); } } static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indirect, double *time) { size_t i; 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; if (p->len > 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; n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); for (i = 0; i < p->len; i++) { x = p->s[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]; t = p->t[i]; if (p->n > 0) { theta0 = p->theta0; } else { theta0 = fmax(MIN_THETA0,p->theta0*sqrt(x)); theta0 = fmin(MAX_THETA0,theta0); } get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o); ts[i] = t*qs[i]; } double dx = p->s[1]-p->s[0]; *mu_direct = trapz(qs,dx,p->len); *mu_indirect = trapz(q_indirects,dx,p->len); *time = trapz(ts,dx,p->len); } 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 n_d2o, double n_h2o, 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, theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, f_reflected, charge, prob, frac; /* 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 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/n_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)*n_d2o*x*beta0*theta0; theta_pmt = acos(-cos_theta_pmt); f = get_weighted_pmt_response(theta_pmt); f_reflected = get_weighted_pmt_reflectivity(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(); get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; /* Assume the particle is travelling at the speed of light. */ *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; double abs_prob = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic); charge = abs_prob*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); /* Add expected number of photons from electromagnetic shower. */ if (p->shower_photons > 0) charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); if (p->delta_ray_photons > 0) charge += abs_prob*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/n_d2o)*omega/(2*M_PI); *mu_reflected = f_reflected*charge; return 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; betaRootParams *pars; pars = (betaRootParams *) params; T = particle_get_energy(x, pars->p); /* Calculate total energy */ E = T + pars->p->mass; mom = sqrt(E*E - pars->p->mass*pars->p->mass); 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-1, 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 n_d2o, double n_h2o, 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. */ 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 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*n_d2o/SPEED_OF_LIGHT + l_h2o*n_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, maxj; static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_shower, ts_sigma, mu_shower; 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 * come from noise hits. */ mu[i] = mu_noise; continue; } /* Technically we should be computing: * * mu[i] = max(p(q|mu)) * * where by max I mean we find the expected number of PE which * maximizes p(q|mu). However, I don't know of any easy analytical * solution to this. So instead we just compute the number of PE which * maximizes p(q|n). */ for (j = 1; j < MAX_PE; j++) { logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j); if (j == 1 || logp[j] > max_logp) { maxj = j; max_logp = logp[j]; } } mu[i] = maxj; ts[i] = ev->pmt_hits[i].t; } mu_shower = 0; ts_shower = 0; 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].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], &mu_shower, 1, &ts[i], &ts_shower, ts[i], PMT_TTS, &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); } double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, 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, nhit; static double logp[MAX_PE], nll[MAX_PMTS]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; double wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov; double logp_path; double mu_reflected; path *path; static double mu_direct[MAX_PMTS][MAX_VERTICES]; static double mu_shower[MAX_PMTS][MAX_VERTICES]; static double ts[MAX_PMTS][MAX_VERTICES]; static double ts_shower[MAX_PMTS][MAX_VERTICES]; static double ts_sigma[MAX_PMTS][MAX_VERTICES]; static double mu[MAX_PMTS]; double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total; double *x, *pdf; double result; double min_ratio; double pos_a, pos_b; double xlo, xhi; size_t npoints, npoints_shower; double log_pt_fast; if (n > MAX_VERTICES) { fprintf(stderr, "maximum number of vertices is %i\n", MAX_VERTICES); exit(1); } /* Initialize static arrays. */ for (i = 0; i < MAX_PMTS; i++) { mu[i] = 0.0; for (j = 0; j < n; j++) { mu_direct[i][j] = 0.0; mu_shower[i][j] = 0.0; ts[i][j] = 0.0; ts_shower[i][j] = 0.0; ts_sigma[i][j] = PMT_TTS; } } for (j = 0; j < n; j++) { 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; /* 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); beta0 = p0/E0; /* 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; x = malloc(sizeof(double)*npoints_shower); pdf = malloc(sizeof(double)*npoints_shower); 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); } /* 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; } } 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. */ wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); cos_theta_cerenkov = 1/(n_d2o*beta0); sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2)); mu_indirect[j] = 0.0; for (i = 0; i < MAX_PMTS; i++) { 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_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov); ts[i][j] += v[j].t0; mu_indirect[j] += mu_reflected; continue; } double q_indirect; integrate_path(path, i, &mu_direct[i][j], &q_indirect, &result); mu_indirect[j] += q_indirect; if (mu_direct[i][j] > 1e-9) { ts[i][j] = v[j].t0 + result/mu_direct[i][j]; } 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] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov); } integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]); mu_indirect[j] += q_indirect; ts_shower[i][j] = v[j].t0 + result; } if (!fast) { path_free(path); free(x); free(pdf); } particle_free(p); } mu_noise = DARK_RATE*GTVALID*1e-9; mu_indirect_total = 0.0; for (j = 0; j < n; j++) { if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION_ELECTRON/10000.0; else mu_indirect_total += mu_indirect[j]*CHARGE_FRACTION_MUON/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; for (j = 0; j < n; j++) { if (fast) mu[i] += mu_direct[i][j] + mu_noise; else mu[i] += mu_direct[i][j] + mu_shower[i][j] + mu_indirect_total + mu_noise; } } 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[i]); if (fast) log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]); if (ev->pmt_hits[i].hit) { for (j = 1; j < MAX_PE; j++) { logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[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_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]); } 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]; if (fast) { nll[nhit++] = -logp[0]; continue; } 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); } } 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; }