diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 71 |
1 files changed, 66 insertions, 5 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index db59040..e045512 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -333,6 +333,7 @@ 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 mu_direct[MAX_PMTS]; double ts[MAX_PMTS]; @@ -376,15 +377,75 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i]); ts[i] += t0; } else { + /* First, we try to compute the distance along the track where the + * PMT is at the Cerenkov angle. The reason for this is because for + * heavy particles like muons which don't scatter much, the + * probability distribution for getting a photon hit along the + * track looks kind of like a delta function, i.e. the PMT is only + * hit over a very narrow window when the angle between the track + * direction and the PMT is *very* close to the Cerenkov angle + * (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. */ + + /* First, we find the point along the track where the PMT is at the + * Cerenkov angle. */ + SUB(pmt_dir,pmts[i].pos,pos); + /* Compute the distance to the PMT. */ + R = NORM(pmt_dir); + normalize(pmt_dir); + + /* Calculate the cosine of the angle between the track direction and the + * vector to the PMT at the start of the track. */ + cos_theta = DOT(dir,pmt_dir); + /* Compute the angle between the track direction and the PMT. */ + theta = acos(cos_theta); + + /* Compute the Cerenkov angle at the start of the track. */ + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); + theta_cerenkov = acos(1/(n_d2o*beta0)); + + /* Now, we compute the distance along the track where the PMT is at the + * Cerenkov angle. + * + * 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); + F.function = &gsl_muon_charge; - gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals); - mu_direct[i] = result; + + 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; + } + + F.function = &gsl_muon_time; ts[i] = t0; if (mu_direct[i] > 1e-9) { - F.function = &gsl_muon_time; - gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals); - ts[i] += result/mu_direct[i]; + 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]; + } } } |