aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/likelihood.c')
-rw-r--r--src/likelihood.c80
1 files changed, 38 insertions, 42 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 45274b7..379ac9a 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -210,7 +210,7 @@ void particle_free(particle *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;
+ double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, theta_pmt, charge, constant, prob_abs, prob_sct;
SUB(pmt_dir,pmt_pos,pos);
@@ -234,15 +234,10 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir,
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();
+ 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);
- 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);
+ 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
@@ -256,19 +251,17 @@ static double get_expected_charge_shower(particle *p, double *pos, double *dir,
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);
+ *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
+ *q_indirect_delta_ray = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*(*q_delta_ray) + prob_sct*(*q_delta_ray);
- *q_delta_ray *= f*scatter;
+ *q_delta_ray *= (1.0-prob_abs)*(1.0-prob_sct)*f;
- return f*charge*scatter;
+ return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}
static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o)
{
- double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter;
+ double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, theta_pmt, charge, prob_abs, prob_sct;
z = 1.0;
@@ -300,21 +293,29 @@ static double get_expected_charge(double x, double beta, double theta0, double *
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);
+ /* 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);
- scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o);
+ charge = omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
- *reflected = (f_reflec + 1.0 - scatter)*charge;
+ *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
- return f*charge*scatter;
+ return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}
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)
@@ -544,7 +545,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
*
* `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;
+ 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, 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. */
@@ -643,32 +644,27 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
theta_pmt = acos(-cos_theta_pmt);
f = get_weighted_pmt_response(theta_pmt);
- f_reflected = get_weighted_pmt_reflectivity(theta_pmt);
+ f_reflec = 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();
+ 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_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);
+ charge = 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);
+ charge += 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);
+ 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/n_d2o)*omega/(2*M_PI);
- *mu_reflected = f_reflected*charge;
+ *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
- return f*charge;
+ return (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}
typedef struct betaRootParams {
generated by cgit (git 2.34.1) at 2025-05-06 21:44:43 +0000