diff options
| -rw-r--r-- | src/likelihood.c | 69 |
1 files changed, 42 insertions, 27 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 5f11f95..43fe336 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -333,8 +333,10 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp; double tmean = 0.0; muon_energy *m; - double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, theta_cerenkov, s; + double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o; double logp_path; + double a, b; + double tmp[3]; double mu_direct[MAX_PMTS]; double ts[MAX_PMTS]; @@ -388,13 +390,10 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl * (it's not a perfect delta function since there is some width due * to dispersion). In this case, it's possible that the numerical * integration completely skips over the delta function and so - * predicts an expected charge of 0. To fix this, we compute the - * integral in two steps, one up to the point along the track where - * the PMT is at the Cerenkov angle and another from that point to - * the end of the track. Since the integration routine always - * samples points near the beginning and end of the integral, this - * allows the routine to correctly compute that the integral is non - * zero. */ + * predicts an expected charge of 0. + * + * To fix this, we only integrate 1 meter on both sides of the + * point along the track where the PMT is at the Cerenkov angle. */ /* First, we find the point along the track where the PMT is at the * Cerenkov angle. */ @@ -422,31 +421,47 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl * point along the track that we want to find, and the PMT position. */ s = R*sin(theta_cerenkov-theta)/sin(theta_cerenkov); + /* Make sure that the point is somewhere along the track between 0 and + * `smax`. */ + if (s < 0) s = 0.0; + else if (s > smax) s = smax; + + /* Compute the endpoints of the integration. */ + a = s - 100.0; + b = s + 100.0; + + if (a < 0.0) a = 0.0; + if (b > range) b = range; + F.function = &gsl_muon_charge; - if (s > 0 && s < range) { - gsl_integration_cquad(&F, 0, s, 0, epsrel, w, &result, &error, &nevals); - mu_direct[i] = result; - gsl_integration_cquad(&F, s, range, 0, epsrel, w, &result, &error, &nevals); - mu_direct[i] += result; - } else { - gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals); - mu_direct[i] = result; - } + gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); + mu_direct[i] = result; F.function = &gsl_muon_time; - ts[i] = t0; if (mu_direct[i] > 1e-9) { - if (s > 0 && s < range) { - gsl_integration_cquad(&F, 0, s, 0, epsrel, w, &result, &error, &nevals); - ts[i] += result/mu_direct[i]; - gsl_integration_cquad(&F, s, range, 0, epsrel, w, &result, &error, &nevals); - ts[i] += result/mu_direct[i]; - } else { - gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals); - ts[i] += result/mu_direct[i]; - } + gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); + ts[i] = t0 + result/mu_direct[i]; + } else { + /* If the expected number of PE is very small, our estimate of + * the time can be crazy since the error in the total charge is + * in the denominator. Therefore, we estimate the time using + * the same technique as in get_total_charge_approx(). See that + * function for more details. */ + n_h2o = get_index_snoman_h2o(wavelength0); + + /* Compute the position of the particle at a distance `s` along the track. */ + tmp[0] = pos[0] + s*dir[0]; + tmp[1] = pos[1] + s*dir[1]; + tmp[2] = pos[2] + s*dir[2]; + + SUB(pmt_dir,pmts[i].pos,tmp); + + get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); + + /* Assume the particle is travelling at the speed of light. */ + ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; } } |