aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/event.h2
-rw-r--r--src/find_peaks.c10
-rw-r--r--src/fit.c10
-rw-r--r--src/likelihood.c6
-rw-r--r--src/test-find-peaks.c2
-rw-r--r--src/test.c4
-rw-r--r--src/zdab_utils.c35
7 files changed, 44 insertions, 25 deletions
diff --git a/src/event.h b/src/event.h
index 1c349e3..1f669b4 100644
--- a/src/event.h
+++ b/src/event.h
@@ -40,6 +40,8 @@ typedef struct pmt_hit {
float ehs;
/* ECA calibrated QLX (pedestal subtracted). */
float elx;
+ /* Best charge. */
+ float q;
/* Integrated charge. */
float qhl;
/* Short-time integrated charge. */
diff --git a/src/find_peaks.c b/src/find_peaks.c
index 135e4e3..7d350d5 100644
--- a/src/find_peaks.c
+++ b/src/find_peaks.c
@@ -151,7 +151,7 @@ double get_qhs_avg(event *ev, double *pos, double *dir)
SUB(pmt_dir,pmts[i].pos,pos);
normalize(pmt_dir);
if (DOT(pmt_dir,dir) > 1/n_d2o) {
- qhs_sum += ev->pmt_hits[i].qhs;
+ qhs_sum += ev->pmt_hits[i].q;
n += 1;
}
}
@@ -190,15 +190,15 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
}
for (i = 0; i < MAX_PMTS; i++) {
- if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
- qhs[i] = ev->pmt_hits[i].qhs;
+ if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+ qhs[i] = ev->pmt_hits[i].q;
}
/* Subtract off previous rings. */
for (i = 0; i < len; i++) {
qhs_avg = get_qhs_avg(ev,pos,last+3*i);
for (j = 0; j < MAX_PMTS; j++) {
- if (pmts[j].pmt_type != PMT_NORMAL || !ev->pmt_hits[j].hit) continue;
+ if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
SUB(pmt_dir,pmts[j].pos,pos);
@@ -216,7 +216,7 @@ void get_hough_transform(event *ev, double *pos, double *x, double *y, size_t n,
for (i = 0; i < n*m; i++) result[i] = 0.0;
for (i = 0; i < MAX_PMTS; i++) {
- if (pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
+ if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
SUB(pmt_dir,pmts[i].pos,pos);
diff --git a/src/fit.c b/src/fit.c
index 5e6a572..c0b9858 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -5155,9 +5155,9 @@ double guess_t0(event *ev, double *pos)
/* Add the charge weighted time difference between the PMT hit time and
* the time it would take a photon to hit the PMT from `pos`. */
- t0 += ev->pmt_hits[i].qhs*(ev->pmt_hits[i].t - distance*n/SPEED_OF_LIGHT);
+ t0 += ev->pmt_hits[i].q*(ev->pmt_hits[i].t - distance*n/SPEED_OF_LIGHT);
- qhs_sum += ev->pmt_hits[i].qhs;
+ qhs_sum += ev->pmt_hits[i].q;
}
/* Divide by the total QHS sum. */
@@ -5195,7 +5195,7 @@ void guess_direction(event *ev, double *pos, double *theta, double *phi)
normalize(pmt_dir);
/* Multiply the vector by the QHS charge in the PMT. */
- MUL(pmt_dir,ev->pmt_hits[i].qhs);
+ MUL(pmt_dir,ev->pmt_hits[i].q);
/* Add this to the estimated direction. */
ADD(dir,dir,pmt_dir);
@@ -5229,7 +5229,7 @@ double guess_energy(event *ev, double *pos, double *dir)
SUB(pmt_dir,pmts[i].pos,pos);
normalize(pmt_dir);
if (DOT(pmt_dir,dir) > 1/n_d2o)
- qhs_sum += ev->pmt_hits[i].qhs;
+ qhs_sum += ev->pmt_hits[i].q;
}
return qhs_sum/6.0;
@@ -5571,7 +5571,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id, double maxtime)
qhs_sum = 0.0;
for (i = 0; i < MAX_PMTS; i++) {
if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue;
- qhs_sum += ev->pmt_hits[i].qhs;
+ qhs_sum += ev->pmt_hits[i].q;
}
T0 = qhs_sum/6.0;
diff --git a/src/likelihood.c b/src/likelihood.c
index 2681475..344f320 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -1103,7 +1103,7 @@ double nll_best(event *ev)
* mu[i] = max(p(q|mu))
*
*/
- mu[i] = get_most_likely_mean_pe(ev->pmt_hits[i].qhs);
+ mu[i] = get_most_likely_mean_pe(ev->pmt_hits[i].q);
/* We want to estimate the mean time which is most likely to produce a
* first order statistic of the actual hit time given there were
* approximately mu[i] PE. As far as I know there are no closed form
@@ -1136,7 +1136,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) + 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);
+ logp[j] = get_log_pq(ev->pmt_hits[i].q,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];
@@ -1423,7 +1423,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, 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) + get_log_phit(j) - mu_sum[i] + j*log_mu - lnfact(j);
+ logp[j] = get_log_pq(ev->pmt_hits[i].q,j) + get_log_phit(j) - mu_sum[i] + j*log_mu - lnfact(j);
if (!charge_only) {
if (fast)
diff --git a/src/test-find-peaks.c b/src/test-find-peaks.c
index f9b7fa4..11fa4b4 100644
--- a/src/test-find-peaks.c
+++ b/src/test-find-peaks.c
@@ -106,7 +106,7 @@ void plot_find_peaks(event *ev, double *pos, double *peak_theta, double *peak_ph
theta = acos(pmts[i].pos[2]/r);
phi = atan2(pmts[i].pos[1],pmts[i].pos[0]);
phi = (phi < 0) ? phi + 2*M_PI: phi;
- fprintf(pipe, "%.10g %.10g %.10g %.10g\n", phi, theta, 0.01, ev->pmt_hits[i].qhs);
+ fprintf(pipe, "%.10g %.10g %.10g %.10g\n", phi, theta, 0.01, ev->pmt_hits[i].q);
}
}
fprintf(pipe,"e\n");
diff --git a/src/test.c b/src/test.c
index 0513598..100091c 100644
--- a/src/test.c
+++ b/src/test.c
@@ -1464,7 +1464,7 @@ int test_quad(char *err)
for (j = 0; j < LEN(hits); j++) {
SUB(pmt_dir,pmts[hits[j]].pos,x0);
ev.pmt_hits[hits[j]].hit = 1;
- ev.pmt_hits[hits[j]].qhs = genrand_real2();
+ ev.pmt_hits[hits[j]].q = genrand_real2();
ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT;
}
@@ -1569,7 +1569,7 @@ int test_quad_noise(char *err)
for (j = 0; j < LEN(hits); j++) {
SUB(pmt_dir,pmts[hits[j]].pos,x0);
ev.pmt_hits[hits[j]].hit = 1;
- ev.pmt_hits[hits[j]].qhs = genrand_real2();
+ ev.pmt_hits[hits[j]].q = genrand_real2();
ev.pmt_hits[hits[j]].t = t0 + NORM(pmt_dir)*n_d2o/SPEED_OF_LIGHT + randn()*PMT_TTS;
}
diff --git a/src/zdab_utils.c b/src/zdab_utils.c
index 6d05d96..e91d31b 100644
--- a/src/zdab_utils.c
+++ b/src/zdab_utils.c
@@ -140,14 +140,25 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
if (bpmt.qms)
ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;
- /* Currently, the charge model only deals with QHS, so we flag any
- * hits which have a bad or railed QHS value.
+ /* Determine the best charge to use in the likelihood function (QHS
+ * or QLX).
*
- * FIXME: In the future, it would be nice to use the best charge
- * word (either QHS or QLX) depending on the if QHS is railed or
- * not, but I need to do more work to see how the QLX values are
- * normalized and if the existing charge model is good enough. */
- if (is_mc(&ev_bank))
+ * Technically, the charge model only deals with QHS and so my
+ * original plan was to just flag any PMT hit with a railed or bad
+ * QHS and not include it in the likelihood function, but when
+ * testing I noticed that muons close to the PSUP can accidentally
+ * get reconstructed as electrons inside the AV if I ignore railed
+ * QHS values, so instead we use the QLX values when QHS is railed.
+ *
+ * Although the charge model only deals with QHS, I think it should
+ * be close enough that it's OK to use QLX. Both charges are
+ * normalized in the same way, the only difference would be in
+ * their SPE width.
+ *
+ * If both QHS and QLX are railed or have uncalibrated charges
+ * below 300, we flag the PMT and it's not included in the
+ * likelihood function. */
+ if (is_mc(&ev_bank)) {
/* Skip uncalibrated charge check if this is MC. The reason is
* that for some reason the uncalibrator in SNOMAN needs the
* QSLP tables which aren't available for most runs. According
@@ -163,9 +174,15 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev)
* > doesn't work, don't try to use it".
*
* Therefore, we ignore it for MC. */
- ;
- else if(bpmt.pihs >= 4095 || bpmt.pihs < 300)
+ ev->pmt_hits[id].q = bpmt.phs;
+ } else if (bpmt.pihs < 4095 && bpmt.pihs >= 300) {
+ ev->pmt_hits[id].q = bpmt.phs;
+ } else if (bpmt.pilx < 4095 && bpmt.pilx >= 300) {
+ ev->pmt_hits[id].q = bpmt.plx;
+ } else {
+ ev->pmt_hits[id].q = 0.0;
ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;
+ }
if (pmts[id].pmt_type != pmt_types[i]) {
get_pmt_type_string(pmts[id].pmt_type,pmt_type_string);
generated by cgit (git 2.34.1) at 2025-05-06 14:01:32 +0000