diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 14 |
1 files changed, 11 insertions, 3 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index fc1cf54..92c8306 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -17,6 +17,7 @@ #include "solid_angle.h" #include <gsl/gsl_roots.h> #include <gsl/gsl_errno.h> +#include "pmt_response.h" typedef struct intParams { path *p; @@ -152,7 +153,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl * * `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, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0; + double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length, wavelength0; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -173,7 +174,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl theta = acos(cos_theta); /* Compute the Cerenkov angle at the start of the track. */ - n = get_index_snoman_d2o(400.0); + wavelength0 = 400.0; + n = get_index_snoman_d2o(wavelength0); theta_cerenkov = acos(1/(n*beta0)); /* Now, we compute the distance along the track where the PMT is at the @@ -236,6 +238,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl SUB(pmt_dir,pmts[i].pos,tmp); + cos_theta_pmt = DOT(pmt_dir,pmts[i].normal)/NORM(pmt_dir); + /* Calculate the sine of the angle between the track direction and the PMT * at the position `s` along the track. */ sin_theta = fabs(sin(acos(DOT(dir,pmt_dir)/NORM(pmt_dir)))); @@ -247,7 +251,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, int i, doubl double frac = sqrt(2)*n*x*beta0*theta0; - return n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI); + f = get_weighted_pmt_response(acos(-cos_theta_pmt)); + + absorption_length = get_absorption_length_snoman_d2o(wavelength0); + + return f*exp(-a/absorption_length)*n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI); } double getKineticEnergy(double x, double T0) |