diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 47 |
1 files changed, 30 insertions, 17 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index aab4ede..e3814a1 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -256,6 +256,13 @@ static double get_expected_charge(double x, double beta, double theta0, double * z = 1.0; + R = NORM(pos); + + n = (R <= AV_RADIUS) ? n_d2o : n_h2o; + + *reflected = 0.0; + if (beta < 1/n) return 0.0; + SUB(pmt_dir,pmt_pos,pos); normalize(pmt_dir); @@ -273,13 +280,6 @@ static double get_expected_charge(double x, double beta, double theta0, double * omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); - R = NORM(pos); - - n = (R <= AV_RADIUS) ? n_d2o : n_h2o; - - *reflected = 0.0; - if (beta < 1/n) return 0.0; - theta_pmt = acos(-cos_theta_pmt); f_reflec = get_weighted_pmt_reflectivity(theta_pmt); f = get_weighted_pmt_response(theta_pmt); @@ -392,7 +392,7 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 static double q_indirects[MAX_NPOINTS]; static double ts[MAX_NPOINTS]; static double ts2[MAX_NPOINTS]; - double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray; + double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray, dx, tmp; if (n > MAX_NPOINTS) { fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); @@ -420,19 +420,26 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 ts2[i] = t*t*qs[i]; } - double dx = x[1]-x[0]; + dx = x[1]-x[0]; *mu_direct = trapz(qs,dx,n); *mu_indirect = trapz(q_indirects,dx,n); - *time = trapz(ts,dx,n)/(*mu_direct); - double t2 = trapz(ts2,dx,n)/(*mu_direct); - *sigma = sqrt(t2-(*time)*(*time)); - - if (*mu_direct == 0.0 || *sigma != *sigma) { + if (*mu_direct == 0.0) { *time = 0.0; *sigma = PMT_TTS; } else { - *sigma = sqrt((*sigma)*(*sigma) + PMT_TTS*PMT_TTS); + *time = trapz(ts,dx,n)/(*mu_direct); + + /* Variance in the time = E(t^2) - E(t)^2. */ + tmp = trapz(ts2,dx,n)/(*mu_direct) - (*time)*(*time); + + if (tmp >= 0) { + *sigma = sqrt(tmp + PMT_TTS*PMT_TTS); + } else { + /* This should never happen but does occasionally due to floating + * point rounding error. */ + *sigma = PMT_TTS; + } } } @@ -542,7 +549,10 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa * 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + if (sin_theta_cerenkov != 0) + s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + else + s = R*cos_theta; /* Make sure that the point is somewhere along the track between 0 and * `smax`. */ @@ -738,7 +748,10 @@ static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, * 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + if (sin_theta_cerenkov != 0) + s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + else + s = R*cos_theta; /* Make sure that the point is somewhere along the track between 0 and * `smax`. */ |