diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 7 |
1 files changed, 5 insertions, 2 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 99258a9..000389b 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -146,7 +146,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy * * `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, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, wavelength0; + double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -247,13 +247,16 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0); absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0); + 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; - return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o)*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); + return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*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); } typedef struct betaRootParams { |