diff options
| author | tlatorre <tlatorre@uchicago.edu> | 2018-09-26 08:59:57 -0500 |
|---|---|---|
| committer | tlatorre <tlatorre@uchicago.edu> | 2018-09-26 08:59:57 -0500 |
| commit | e5c270c8c442136c7fd44751e1ccac887f49bd06 (patch) | |
| tree | 854259150133a810c68808e435eae45f30eb3378 /src/likelihood.c | |
| parent | 8bdcdbc04913df2d0d4ce4d1d8d85c8da2185e59 (diff) | |
| download | sddm-e5c270c8c442136c7fd44751e1ccac887f49bd06.tar.gz sddm-e5c270c8c442136c7fd44751e1ccac887f49bd06.tar.bz2 sddm-e5c270c8c442136c7fd44751e1ccac887f49bd06.zip | |
speed up fast likelihood calculation
This commit speeds up the fast likelihood calculation by avoiding calls to
trigonometric functions where possible. Specifically we calculate
sin(a) = sqrt(1-pow(cos(a),2));
instead of
sin(a) = sin(acos(cos(a)));
Diffstat (limited to 'src/likelihood.c')
| -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); |