diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 12 |
1 files changed, 6 insertions, 6 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 66e72e4..012243c 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -124,10 +124,10 @@ double getKineticEnergy(double x, double T0) double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel) { - size_t i, j; + size_t i, j, nhit; intParams params; double total_charge; - double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0; + double logp[MAX_PE], nll[MAX_PMTS], range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0; double tmean = 0.0; int npmt = 0; @@ -236,7 +236,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl mu[i] = mu_direct[i] + mu_indirect + mu_noise; } - nll = 0; + nhit = 0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue; @@ -245,15 +245,15 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl for (j = 1; j < MAX_PE; j++) { logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], tmean, 1.5); } - nll -= logsumexp(logp, sizeof(logp)/sizeof(double)); + nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double)); } else { logp[0] = -mu[i]; for (j = 1; j < MAX_PE; j++) { logp[j] = log(get_pmiss(j)) - mu[i] + j*log(mu[i]) - gsl_sf_lnfact(j); } - nll -= logsumexp(logp, sizeof(logp)/sizeof(double)); + nll[nhit++] = -logsumexp(logp, sizeof(logp)/sizeof(double)); } } - return nll; + return kahan_sum(nll,nhit); } |