diff options
Diffstat (limited to 'src/sno_charge.c')
| -rw-r--r-- | src/sno_charge.c | 35 |
1 files changed, 35 insertions, 0 deletions
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; |