diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/event.h | 2 | ||||
| -rw-r--r-- | src/find_peaks.c | 10 | ||||
| -rw-r--r-- | src/fit.c | 10 | ||||
| -rw-r--r-- | src/likelihood.c | 6 | ||||
| -rw-r--r-- | src/test-find-peaks.c | 2 | ||||
| -rw-r--r-- | src/test.c | 4 | ||||
| -rw-r--r-- | src/zdab_utils.c | 35 |
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); @@ -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"); @@ -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); |