diff options
| -rw-r--r-- | src/likelihood.c | 4 | ||||
| -rw-r--r-- | src/sno_charge.c | 35 | ||||
| -rw-r--r-- | src/sno_charge.h | 1 |
3 files changed, 38 insertions, 2 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 442a8b3..f0430cf 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -941,7 +941,7 @@ double nll_best(event *ev) 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], 1, &ts[i], ts[i], &ts_sigma); + logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) + get_log_phit(j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], 1, &ts[i], ts[i], &ts_sigma); if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; @@ -1232,7 +1232,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, const in 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_sum[i] + j*log_mu - lnfact(j); + logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) + get_log_phit(j) - mu_sum[i] + j*log_mu - lnfact(j); if (!charge_only) { if (fast) diff --git a/src/sno_charge.c b/src/sno_charge.c index f78da38..4b51954 100644 --- a/src/sno_charge.c +++ b/src/sno_charge.c @@ -54,6 +54,7 @@ static double qmean, qstd; * init_charge() by integrating the charge PDF up to the mean threshold for all * channels. */ static double log_pmiss[MAX_PE]; +static double log_phit[MAX_PE]; /* Keep track of whether init_charge() has been called so we can exit if a * function is called before being initialized. */ @@ -164,6 +165,28 @@ double get_pq(double q, int n) return interp1d(q,x,pq[n-1],nq); } +double get_log_phit(int n) +{ + /* Returns the log of the probability that the channel discriminator fired + * given `n` PE were generated. */ + if (!initialized) { + fprintf(stderr, "charge interpolation hasn't been initialized!\n"); + exit(1); + } + + if (n == 0) { + /* If no PE were generated, there is a 100% chance that the + * discriminator didn't fire. */ + return -INFINITY; + } else if (n > MAX_PE) { + /* For n > MAX_PE we assume there is a 100% chance that the + * discriminator fires. */ + return 0.0; + } + + return log_phit[n-1]; +} + double get_log_pmiss(int n) { /* Returns the log of the probability that the channel discriminator didn't @@ -290,6 +313,18 @@ void init_charge(void) params[0] = i; gsl_integration_cquad(&F, 0, MEAN_THRESH/MEAN_HIPT, 0, 1e-9, w, &result, &error, &nevals); log_pmiss[i-1] = log(result); + log_phit[i-1] = log(1.0-result); + } + + /* Set the PDF for the charge to 0 before the discriminator level. */ + for (i = 1; i <= MAX_PE; i++) { + for (j = 0; j < nq; j++) { + if (x[j] < MEAN_THRESH/MEAN_HIPT) pq[i-1][j] = 0.0; + } + double norm = trapz(&pq[i-1][0],x[1]-x[0],nq); + for (j = 0; j < nq; j++) { + pq[i-1][j] /= norm; + } } F.function = &gsl_smear; diff --git a/src/sno_charge.h b/src/sno_charge.h index 7f7caa2..141ad28 100644 --- a/src/sno_charge.h +++ b/src/sno_charge.h @@ -20,6 +20,7 @@ void init_charge(void); double get_log_pq(double q, int n); double get_pq(double q, int n); +double get_log_phit(int n); double get_log_pmiss(int n); double spe_pol2exp(double q); |