diff options
| -rw-r--r-- | src/likelihood.c | 12 |
1 files changed, 7 insertions, 5 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index ef6622d..f78eee8 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, theta, x, z, s, a, b, beta, E, p, T, omega, 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, theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, 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; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -170,7 +170,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); - theta_cerenkov = acos(1/(n_d2o*beta0)); + cos_theta_cerenkov = 1/(n_d2o*beta0); + theta_cerenkov = acos(cos_theta_cerenkov); /* Now, we compute the distance along the track where the PMT is at the * Cerenkov angle. @@ -178,7 +179,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy * 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-theta)/sin(theta_cerenkov); + s = R*sin(theta_cerenkov-theta)/sqrt(1-pow(cos_theta_cerenkov,2)); /* Make sure that the point is somewhere along the track between 0 and * `smax`. */ @@ -203,7 +204,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy z = R*cos_theta; /* `x` is the perpendicular distance from the PMT position to the track. */ - x = R*fabs(sin(theta)); + x = R*sqrt(1-pow(cos_theta,2)); /* `b` is the second coefficient in the Taylor expansion. */ b = (s-z)/a; @@ -233,7 +234,8 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy /* 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)))); + cos_theta = DOT(dir,pmt_dir)/NORM(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_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); |