From b4be76539cbc0254df33e8080f965df3a4ce92a4 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 25 Sep 2018 10:21:21 -0500 Subject: update integration bounds in likelihood calculation This commit updates the bounds of the track integration in the likelihood function to integrate up to 1 meter around the point at which the PMT is at the Cerenkov angle from the track. This fixes an issue I was seeing where a *very* small change in the fit paramters would cause the likelihood to jump by a large amount. I eventually tracked it down to the same issue I was seeing before which I solved by splitting up the integration into two intervals. However that fix did not seem to completely fix the issue. Based on initial tests with 500 MeV muons, this fix seems to do a much better job. --- src/likelihood.c | 69 ++++++++++++++++++++++++++++++++++---------------------- 1 file changed, 42 insertions(+), 27 deletions(-) (limited to 'src/likelihood.c') 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; } } -- cgit