diff options
| -rw-r--r-- | src/fit.c | 36 | ||||
| -rw-r--r-- | src/likelihood.c | 111 | ||||
| -rw-r--r-- | src/likelihood.h | 3 | ||||
| -rwxr-xr-x | utils/plot | 13 |
4 files changed, 153 insertions, 10 deletions
@@ -43,6 +43,7 @@ typedef struct fitParams { int fast; int print; int id; + int charge_only; } fitParams; int particles[] = { @@ -4975,7 +4976,7 @@ static struct startingParameters { { 800.0, 800.0, 800.0}, }; -static double nopt_nll(unsigned int n, const double *x, double *grad, void *params) +static double nlopt_nll(unsigned int n, const double *x, double *grad, void *params) { fitParams *fpars = (fitParams *) params; double theta, phi; @@ -5007,7 +5008,7 @@ static double nopt_nll(unsigned int n, const double *x, double *grad, void *para v.n = 1; gettimeofday(&tv_start, NULL); - fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast); + fval = nll(fpars->ev, &v, 1, fpars->dx, fpars->dx_shower, fpars->fast, fpars->charge_only); gettimeofday(&tv_stop, NULL); long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; @@ -5118,7 +5119,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id) struct timeval tv_start, tv_stop; opt = nlopt_create(NLOPT_LN_BOBYQA, 9); - nlopt_set_min_objective(opt,nopt_nll,&fpars); + nlopt_set_min_objective(opt,nlopt_nll,&fpars); /* Make a guess as to the energy. Right now we just use a simple * approximation that the muon produces approximately 6 photons/MeV. @@ -5219,6 +5220,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id) fpars.dx_shower = 10.0; fpars.fast = 1; fpars.print = 0; + fpars.charge_only = 0; nlopt_set_ftol_abs(opt, 1.0); nlopt_set_maxeval(opt, 20); @@ -5313,6 +5315,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id) fpars.dx_shower = 10.0; fpars.fast = 0; fpars.print = 1; + fpars.charge_only = 0; nlopt_set_ftol_abs(opt, 1e-5); nlopt_set_maxeval(opt, 1000); @@ -5408,6 +5411,25 @@ end: return rv; } +size_t get_nhit(event *ev) +{ + /* Returns the number of PMT hits in event `ev`. + * + * Note: Only hits on normal PMTs which aren't flagged are counted. */ + size_t i, nhit; + + nhit = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + + if (!ev->pmt_hits[i].hit) continue; + + nhit++; + } + + return nhit; +} + int main(int argc, char **argv) { int i; @@ -5425,6 +5447,8 @@ int main(int argc, char **argv) FILE *fout = NULL; int skip_second_event = 0; struct timeval tv_start, tv_stop; + long long elapsed; + double fmin_best; for (i = 1; i < argc; i++) { if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) { @@ -5581,6 +5605,7 @@ int main(int argc, char **argv) if (fout) { if (first_ev) fprintf(fout, " ev:\n"); fprintf(fout, " - gtid: %i\n", ev.gtid); + fprintf(fout, " nhit: %zu\n", get_nhit(&ev)); fprintf(fout, " fit:\n"); } @@ -5592,7 +5617,9 @@ int main(int argc, char **argv) } gettimeofday(&tv_stop, NULL); - long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; + elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; + + fmin_best = nll_best(&ev); if (fout) { fprintf(fout, " %i:\n", particles[i]); @@ -5607,6 +5634,7 @@ int main(int argc, char **argv) fprintf(fout, " z2: %.2f\n", xopt[8]); fprintf(fout, " fmin: %.2f\n", fmin); fprintf(fout, " time: %lld\n", elapsed); + fprintf(fout, " psi: %.2f\n", fmin-fmin_best); fflush(fout); } } diff --git a/src/likelihood.c b/src/likelihood.c index 0e78800..64f28d4 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -720,7 +720,104 @@ static double getKineticEnergy(double x, void *p) return particle_get_energy(x, (particle *) p); } -double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast) +double nll_best(event *ev) +{ + /* Returns the negative log likelihood of the "best hypothesis" for the event `ev`. + * + * By "best hypothesis" I mean we restrict the model to assign a single + * mean number of PE for each PMT in the event assuming that the number of + * PE hitting each PMT is poisson distributed. In addition, the model picks + * a single time for the light to arrive with the PDF being given by a sum + * of dark noise and a Gaussian around this time with a width equal to the + * single PE transit time spread. + * + * This calculation is intended to be used as a goodness of fit test for + * the fitter by computing a likelihood ratio between the best fit + * hypothesis and this ideal case. See Chapter 9 in Jayne's "Probability + * Theory: The Logic of Science" for more details. */ + size_t i, j, nhit, maxj; + static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_shower, ts_sigma, mu_shower; + double log_mu, max_logp, mu_noise, mu_indirect_total, min_ratio; + + mu_noise = DARK_RATE*GTVALID*1e-9; + + /* Compute the "best" number of expected PE for each PMT. */ + for (i = 0; i < MAX_PMTS; i++) { + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + + if (!ev->pmt_hits[i].hit) { + /* If the PMT wasn't hit we assume the expected number of PE just + * come from noise hits. */ + mu[i] = mu_noise; + continue; + } + + /* Technically we should be computing: + * + * mu[i] = max(p(q|mu)) + * + * where by max I mean we find the expected number of PE which + * maximizes p(q|mu). However, I don't know of any easy analytical + * solution to this. So instead we just compute the number of PE which + * maximizes p(q|n). */ + for (j = 1; j < MAX_PE; j++) { + logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j); + + if (j == 1 || logp[j] > max_logp) { + maxj = j; + max_logp = logp[j]; + } + } + + mu[i] = maxj; + ts[i] = ev->pmt_hits[i].t; + } + + mu_shower = 0; + ts_shower = 0; + ts_sigma = PMT_TTS; + + mu_indirect_total = 0.0; + + min_ratio = MIN_RATIO; + + /* Now, we actually compute the negative log likelihood for the best + * hypothesis. + * + * Currently this calculation is identical to the one in nll() so it should + * probably be made into a function. */ + nhit = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + + log_mu = log(mu[i]); + + if (ev->pmt_hits[i].hit) { + for (j = 1; j < MAX_PE; j++) { + logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], &mu_shower, 1, &ts[i], &ts_shower, ts[i], PMT_TTS, &ts_sigma); + + if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; + + if (logp[j] - max_logp < min_ratio*ln(10)) { + j++; + break; + } + } + + nll[nhit++] = -logsumexp(logp+1, j-1); + } else { + logp[0] = -mu[i]; + for (j = 1; j < MAX_PE_NO_HIT; j++) { + logp[j] = get_log_pmiss(j) - mu[i] + j*log_mu - lnfact(j); + } + nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT); + } + } + + return kahan_sum(nll,nhit); +} + +double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only) { /* Returns the negative log likelihood for event `ev` given a particle with * id `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and @@ -959,10 +1056,14 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast if (ev->pmt_hits[i].hit) { for (j = 1; j < MAX_PE; j++) { - if (fast) - logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast; - else - logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]); + logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j); + + if (!charge_only) { + if (fast) + logp[j] += log_pt_fast; + else + logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]); + } if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; diff --git a/src/likelihood.h b/src/likelihood.h index 69d6896..1ef2736 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -83,6 +83,7 @@ double particle_get_energy(double x, particle *p); void particle_free(particle *p); double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma); double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma); -double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast); +double nll_best(event *ev); +double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast, int charge_only); #endif @@ -72,10 +72,14 @@ if __name__ == '__main__': likelihood_ratio = [] t_electron = [] t_muon = [] + psi = [] for event in data['data']: # get the particle ID id = event['mctk'][-1]['id'] + if 'ev' not in event: + continue + if id not in event['ev'][0]['fit']: continue @@ -115,6 +119,11 @@ if __name__ == '__main__': t_electron.append(event['ev'][0]['fit'][IDP_E_MINUS]['time']) if IDP_MU_MINUS in event['ev'][0]['fit']: t_muon.append(event['ev'][0]['fit'][IDP_MU_MINUS]['time']) + nhit = event['ev'][0]['nhit'] + psi.append(event['ev'][0]['fit'][id]['psi']/nhit) + + if len(t_electron) < 2 or len(t_muon) < 2: + continue mean, mean_error, std, std_error = get_stats(dT) print("dT = %.2g +/- %.2g" % (mean, mean_error)) @@ -155,6 +164,9 @@ if __name__ == '__main__': plt.figure(8) plot_hist(np.array(t_muon)/1e3/60.0, label=filename) plt.xlabel(r"Muon Fit time (minutes)") + plt.figure(9) + plot_hist(psi, label=filename) + plt.xlabel(r"$\Psi$/Nhit") plot_legend(1) plot_legend(2) @@ -164,4 +176,5 @@ if __name__ == '__main__': plot_legend(6) plot_legend(7) plot_legend(8) + plot_legend(9) plt.show() |