aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/likelihood.c')
-rw-r--r--src/likelihood.c34
1 files changed, 17 insertions, 17 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 3f27915..25563af 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -202,7 +202,7 @@ particle *particle_init(int id, double T0, size_t n)
}
if (p->shower_photons) {
- norm = electron_get_angular_pdf_norm(p->a, p->b, 1/avg_index_d2o);
+ 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);
@@ -215,7 +215,7 @@ particle *particle_init(int id, double T0, size_t n)
*
* 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/avg_index_d2o)/norm + EPSILON;
+ 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
@@ -236,7 +236,7 @@ particle *particle_init(int id, double T0, size_t n)
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/avg_index_d2o);
+ 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);
@@ -249,7 +249,7 @@ particle *particle_init(int id, double T0, size_t n)
*
* 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/avg_index_d2o)/norm + EPSILON;
+ 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
@@ -521,7 +521,7 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
*reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
- *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+ *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;
}
@@ -553,7 +553,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
R = NORM(pos);
- n = (R <= AV_RADIUS) ? avg_index_d2o : avg_index_h2o;
+ n = (R <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
*q = 0.0;
*t = 0.0;
@@ -583,7 +583,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
- *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+ *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:
*
@@ -872,7 +872,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
*t = 0.0;
*mu_reflected = 0.0;
- if (beta < 1/avg_index_d2o) return 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. */
@@ -903,7 +903,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
- frac = sqrt(2)*avg_index_d2o*x*beta0*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);
@@ -914,15 +914,15 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
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*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/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 = avg_index_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+avg_index_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+avg_index_d2o*z*beta0)/frac))/(b+avg_index_d2o*beta0)/(4*M_PI);
+ 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/avg_index_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/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/avg_index_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/get_avg_index_d2o())*omega/(2*M_PI);
*mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
@@ -1049,7 +1049,7 @@ static double guess_time(double *pos, double *dir, double *pmt_pos, double smax,
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*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/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)
@@ -1273,8 +1273,8 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in
/* 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);
+ if (beta0 > 1/get_avg_index_d2o())
+ get_smax(p, 1/get_avg_index_d2o(), range, &cerenkov_range);
else
cerenkov_range = range;
@@ -1298,7 +1298,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in
*
* 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/(avg_index_d2o*beta0);
+ cos_theta_cerenkov = 1/(get_avg_index_d2o()*beta0);
sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2));
mu_indirect[j] = 0.0;
generated by cgit (git 2.34.1) at 2025-05-06 18:51:17 +0000