diff options
| -rw-r--r-- | src/electron.c | 2 | ||||
| -rw-r--r-- | src/likelihood.c | 47 | ||||
| -rw-r--r-- | src/muon.c | 4 |
3 files changed, 33 insertions, 20 deletions
diff --git a/src/electron.c b/src/electron.c index 73fe24e..60bcec7 100644 --- a/src/electron.c +++ b/src/electron.c @@ -98,7 +98,7 @@ double electron_get_angular_distribution_alpha(double T0) * * where T0 is the initial energy of the electron in MeV and c0, c1, c2, * and c3 are constants which I fit out. */ - return 3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00); + return fmin(100.0,3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00)); } double electron_get_angular_distribution_beta(double T0) 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`. */ @@ -98,7 +98,7 @@ double muon_get_angular_distribution_alpha(double T0) * * where T0 is the initial energy of the muon in MeV and c0, c1, c2, * and c3 are constants which I fit out. */ - return 8.238633e-01 + 3.896665e-03/log(1.581060e-05*T0 + 9.991576e-01); + return fmin(100.0,8.238633e-01 + 3.896665e-03/log(1.581060e-05*T0 + 9.991576e-01)); } double muon_get_angular_distribution_beta(double T0) @@ -144,7 +144,7 @@ void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b) * * where T0 is the initial energy of the muon in MeV and c0, c1, c2, * and c3 are constants which I fit out. */ - *a = 3.463093e-01 + 1.110835e-02/log(5.662948e-06*T0 + 1.009215e+00); + *a = fmin(100.0,3.463093e-01 + 1.110835e-02/log(5.662948e-06*T0 + 1.009215e+00)); *b = 2.297358e-01 + 4.085721e-03/log(8.218774e-06*T0 + 1.007135e+00); return; |