diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/dc.c | 371 | ||||
| -rw-r--r-- | src/dc.h | 12 | ||||
| -rw-r--r-- | src/fit.c | 7 | ||||
| -rw-r--r-- | src/zdab_utils.c | 8 |
4 files changed, 314 insertions, 84 deletions
@@ -24,6 +24,98 @@ #include <gsl/gsl_sort.h> #include "zdab_utils.h" #include "zebra.h" +#include <gsl/gsl_statistics_double.h> + +/* Returns 1 if the event is tagged as an incoming muon. The goal of this cut + * is to tag external muons *without* tagging muons which start within the PSUP + * and then exit. In order to achieve that we use the ECA time of the OWL + * tubes. + * + * Ideally we would have the full calibrated time of the OWL tubes, but there + * was never any PCA for the OWL tubes. I also discussed this with Josh and he + * said that the cables were never changed for the OWL tubes so we should at + * least be able to rely on them being consistent. + * + * Since the average PCA correction is of order +/- 5 ns and the round trip + * time for light from the edge of the PSUP to the cavity and back should be + * longer than 15 ns (4 meters), the ECA time should be enough to distinguish + * incoming vs outgoing muons. + * + * This function tags an event as an incoming muon if the following criteria + * are met: + * + * - calibrated nhit >= 150 + * - at least 5 calibrated OWL hits + * - at least 1 OWL hit with a QHS >= 20.0 + * - at least 2 OWL hits with ECA times earlier than the 10th percentile of the + * PSUP PMT times + * + * Idea: Should we also require at least one of the early tubes to have a high + * charge? */ +int is_muon(event *ev) +{ + size_t i; + double ept[MAX_PMTS]; + size_t nhit, nhit_owl, nhit_owl_early; + double early_time, max_owl_qhs; + + nhit = 0; + nhit_owl = 0; + for (i = 0; i < MAX_PMTS; i++) { + /* Require good calibrations. + * + * Note: Typically I would do this with the following check: + * + * ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL) + * + * However, Javi mentioned in an email that for very early hits which + * are in the TAC curl region, the calibration processor will flag + * them, set the time to -9999.0, and set the KPF_BAD_CAL bit. However, + * in this case for some reason the ECA time still looks reasonable. + * Therefore, since we only use the ECA time here we instead just make + * sure that the ECA time is something reasonable. */ + if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue; + + switch (pmts[i].pmt_type) { + case PMT_OWL: + if (nhit_owl == 0) { + max_owl_qhs = ev->pmt_hits[i].qhs; + } else if (ev->pmt_hits[i].qhs > max_owl_qhs) { + max_owl_qhs = ev->pmt_hits[i].qhs; + } + nhit_owl += 1; + break; + case PMT_NORMAL: + if (!ev->pmt_hits[i].flags) ept[nhit++] = ev->pmt_hits[i].ept; + break; + } + } + + /* Require at least 150 normal hits. */ + if (nhit < MUON_MIN_NHIT) return 0; + + /* Require at least 5 OWL hits. */ + if (nhit_owl < MUON_MIN_OWL_NHIT) return 0; + + /* Require at least 1 OWL tube to have a QHS >= 20.0. */ + if (max_owl_qhs < MUON_OWL_QHS) return 0; + + /* Sort the times. */ + gsl_sort(ept,1,nhit); + + /* Calculate the time to define early OWL hits. */ + early_time = gsl_stats_quantile_from_sorted_data(ept,1,nhit,MUON_OWL_TIME_PERCENTILE); + + nhit_owl_early = 0; + for (i = 0; i < MAX_PMTS; i++) { + /* Require good calibrations. */ + if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue; + + if (pmts[i].pmt_type == PMT_OWL && ev->pmt_hits[i].ept < early_time) nhit_owl_early += 1; + } + + return nhit_owl_early >= 2; +} /* Returns 1 if the event fails the junk cut. The definition for the junk cut * comes from the SNOMAN companion: @@ -31,11 +123,14 @@ * This will remove orphans, ECA data and events containing the same PMT * more than once. * + * I also added a check that will flag the event if it contains more than + * MAX_PMTS PMT hits (which shouldn't be possible). + * * Note: This function may return -1 if there is an error traversing the PMT * banks. */ int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev) { - size_t i; + int i; int ids[MAX_PMTS]; int nhit; static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK}; @@ -68,7 +163,14 @@ int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev) while (1) { unpack_pmt(b.data, &bpmt); + ids[nhit++] = bpmt.pin; + + if (nhit >= MAX_PMTS) { + fprintf(stderr, "event with more than %i PMT hits\n", MAX_PMTS); + return 1; + } + if (!b.next) break; rv = zebra_get_bank(f,&b,b.next); @@ -240,14 +342,17 @@ int is_neck_event(event *ev) double avg_ept_psup; int nhit; - /* First, we compute the average ECA time for all normal PMTs. */ + /* First, we compute the average ECA time for all the PSUP PMTs with z <= 0. */ nhit = 0; avg_ept_psup = 0.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 || ev->pmt_hits[i].ept <= -100) continue; - if (ev->pmt_hits[i].hit) { + /* Require good calibrations. */ + if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue; + + if (ev->pmt_hits[i].hit && pmts[i].pos[2] < 0) { avg_ept_psup += ev->pmt_hits[i].ept; nhit += 1; } @@ -263,7 +368,10 @@ int is_neck_event(event *ev) n = 0; high_charge = 0; for (i = 0; i < MAX_PMTS; i++) { - if (pmts[i].pmt_type == PMT_NECK && ev->pmt_hits[i].hit) { + /* Require good calibrations. */ + if (!ev->pmt_hits[i].hit || ev->pmt_hits[i].ept <= -100) continue; + + if (pmts[i].pmt_type == PMT_NECK) { n += 1; /* FIXME: The SNOMAN companion says "Early if defined by the neck @@ -280,6 +388,81 @@ int is_neck_event(event *ev) return 0; } +/* Returns 1 if the average time for the hits in the same card as + * `flasher_pmt_id` is at least 10 ns earlier than all the hits within 4 + * meters of the potential flasher. + * + * This cut is designed to help flag flashers in the is_flasher() function for + * events in which the flashing channel may be missing from the event, so there + * is no high charge channel to flag it. */ +int is_slot_early(event *ev, int flasher_pmt_id) +{ + int i; + double t_slot, t_nearby; + int n_slot, n_nearby; + double flasher_pos[3]; + double pmt_dir[3]; + double distance; + + COPY(flasher_pos,pmts[flasher_pmt_id].pos); + + n_slot = 0; + n_nearby = 0; + t_slot = 0.0; + t_nearby = 0.0; + for (i = 0; i < MAX_PMTS; i++) { + if (!ev->pmt_hits[i].hit) continue; + + /* Require good calibrations. + * + * Note: For some reason PMT hits in the flasher slot sometimes are + * marked as having bad calibrations and get a time of -9999. I still + * don't understand why this is, so for now we don't check the + * calibration bits and instead just check to see that the ECA time is + * reasonable. + * + * Here is what Javi said via email: + * + * If SNO calibrations are the same than SNO+ (and I'm pretty sure they + * are, at last for this matter) there is a hit-level time validation, + * i.e., if the TAC value is determined to be in the curl region of the + * time slope (high TAC = very early hits), they get flagged and a + * calibration is not provided (which is compatible with a flashing + * channel). If this is the case, you can fall back to raw TACs and QHS + * when calibration is not available. */ + if (ev->pmt_hits[i].ept < -100) continue; + + if (i/32 == flasher_pmt_id/32) { + /* This hit is in the same slot as the potential flasher. */ + t_slot += ev->pmt_hits[i].ept; + n_slot += 1; + continue; + } + + /* Calculate the distance from the current channel to the high charge + * channel. */ + SUB(pmt_dir,pmts[i].pos,flasher_pos); + distance = NORM(pmt_dir); + + if (distance < 400.0) { + t_nearby += ev->pmt_hits[i].ept; + n_nearby += 1; + continue; + } + } + + /* If there are no PMTs within 4 meters of the slot, it's almost certainly + * a flasher. */ + if (n_slot > 4 && n_nearby == 0) return 1; + + if (n_slot == 0 || n_nearby == 0) return 0; + + t_slot /= n_slot; + t_nearby /= n_nearby; + + return t_slot < t_nearby + 10.0; +} + /* Returns 1 if the event is a flasher and 0 if it isn't. The definition of * a flasher event comes from * http://detector.sno.laurentian.ca/~detector/private/snoop/doc/eventIDs.html, @@ -317,18 +500,19 @@ int is_neck_event(event *ev) int is_flasher(event *ev) { int hits_pc[1280] = {0}; - int qhs, qhl; + int qhs, qhl, qlx; int qhl_pc[8] = {0}; int qhs_pc[8] = {0}; + int qlx_pc[8] = {0}; double t_pc[8] = {0}; + int channel_pc[8] = {0}; size_t crate, card, channel, id; int i, j; - size_t index[1280], index_qhs[8], index_qhl[8]; + size_t index[1280], index_qhs[8], index_qhl[8], index_qlx[8]; double t; int flasher_pc; double flasher_pos[3]; double pmt_dir[3]; - int flasher = 0; double distance; int nhit, nhit_late, nhit_far; @@ -358,12 +542,12 @@ int is_flasher(event *ev) for (j = 1279; j >= 0; j--) { id = index[j]; - /* If this paddle card doesn't have as many hits as the most hit paddle - * card and we haven't found a flasher yet, then it's not a flasher. */ - if (hits_pc[id] < hits_pc[index[1279]]) return 0; + /* Make sure this paddle card has at least 4 hits. */ + if (hits_pc[id] < 4) return 0; + nhit = 0; for (i = 0; i < MAX_PMTS; i++) { - if (!ev->pmt_hits[i].hit) continue; + if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue; crate = i/512; card = (i % 512)/32; @@ -376,101 +560,122 @@ int is_flasher(event *ev) /* Get the uncalibrated QHS and QHL charges. */ qhs = ev->pmt_hits[i].qihs; qhl = ev->pmt_hits[i].qihl; + qlx = ev->pmt_hits[i].qilx; if (qhs < 300) { /* QHS values below 300 are set to 4095. */ - qhs_pc[channel % 8] = 4095; + qhs_pc[nhit] = 4095; } else { - qhs_pc[channel % 8] = qhs; + qhs_pc[nhit] = qhs; } if (qhl < 300) { /* QHL values below 300 are set to 4095. */ - qhl_pc[channel % 8] = 4095; + qhl_pc[nhit] = 4095; } else { - qhl_pc[channel % 8] = qhl; + qhl_pc[nhit] = qhl; } - t_pc[channel % 8] = ev->pmt_hits[i].t; - } + if (qlx < 300) { + /* QLX values below 300 are set to 4095. */ + qlx_pc[nhit] = 4095; + } else { + qlx_pc[nhit] = qlx; + } - /* Sort the QHS and QHL values by size. */ - argsort(qhs_pc, LEN(qhs_pc), index_qhs); - argsort(qhl_pc, LEN(qhl_pc), index_qhl); - - /* Check if this paddle card has QHS or QHL values > 4000, and the QHS - * and QHL values for the other channels in the card are less than 700. */ - if (qhs_pc[index_qhs[7]] > 4000 || qhl_pc[index_qhl[7]] > 4000) { - if (qhs_pc[index_qhs[6]] < 800 && qhl_pc[index_qhl[6]] < 800) { - /* Flasher! */ - flasher = 1; - break; - } + t_pc[nhit] = ev->pmt_hits[i].t; + + channel_pc[nhit] = channel; + + nhit += 1; } - } - /* If we didn't find a flasher, return 0. */ - if (!flasher) return 0; - - flasher_pc = index[j]; - crate = index[j]/64; - card = (index[j] % 64)/4; - if (qhs_pc[index_qhs[7]] > 4000) { - channel = (index[j] % 4)*8 + index_qhs[7]; - t = t_pc[index_qhs[7]]; - } else { - channel = (index[j] % 4)*8 + index_qhl[7]; - t = t_pc[index_qhl[7]]; - } - COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos); + if (nhit < 2) { + fprintf(stderr, "only 2 hits on %zu/%zu/PC%zu but expected at least 4!\n", + id/64, (id % 64)/4, id % 4); + continue; + } - /* Check that 70% of the regular PMTs in the event which are not in the - * flasher slot fired more than 50 ns after the high charge channel and are - * more than 12 meters away from the high charge channel. */ - nhit = 0; - /* Number of PMTs which fired 50 ns after the high charge channel. */ - nhit_late = 0; - /* Number of PMTs which are 12 meters or more away from the high charge - * channel. */ - nhit_far = 0; - for (i = 0; i < MAX_PMTS; i++) { - /* Skip PMTs which weren't hit. */ - if (!ev->pmt_hits[i].hit) continue; + /* Sort the QHS and QHL values by size. */ + argsort(qhs_pc, nhit, index_qhs); + argsort(qhl_pc, nhit, index_qhl); + argsort(qlx_pc, nhit, index_qlx); - /* Only count regular PMTs without any flags. */ - if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + flasher_pc = index[j]; + crate = index[j]/64; + card = (index[j] % 64)/4; + channel = channel_pc[index_qhs[nhit-1]]; - crate = i/512; - card = (i % 512)/32; - channel = i % 32; + /* Check if this paddle card has QHS, QHL, or QLX values > 4000, and the QHS, + * QHL, and QLX values for the other channels in the card are less than 700. + * + * Alternatively, if there isn't an obvious charge outlier, it's still + * possible that the flasher PMT hit is missing from the event or that + * the QHS/QHL/QLX values were just outside of the threshold. In this + * case, we look to see if the TAC for the proposed slot is early with + * respect to neighboring PMTs. */ + if ((qhs_pc[index_qhs[nhit-1]] < 4000 || qhs_pc[index_qhs[nhit-2]] > 800) && + (qhl_pc[index_qhl[nhit-1]] < 4000 || qhl_pc[index_qhl[nhit-2]] > 800) && + (qlx_pc[index_qlx[nhit-1]] < 4000 || qlx_pc[index_qlx[nhit-2]] > 800) && + !is_slot_early(ev, crate*512 + card*32 + channel)) continue; + + if (qhs_pc[index_qhs[nhit-1]] > 4000) { + channel = channel_pc[index_qhs[nhit-1]]; + t = t_pc[index_qhs[nhit-1]]; + } else if (qhl_pc[index_qhl[nhit-1]] > 4000) { + channel = channel_pc[index_qhl[nhit-1]]; + t = t_pc[index_qhl[nhit-1]]; + } else { + channel = channel_pc[index_qlx[nhit-1]]; + t = t_pc[index_qlx[nhit-1]]; + } + COPY(flasher_pos,pmts[crate*512 + card*32 + channel].pos); + + /* Check that 70% of the regular PMTs in the event which are not in the + * flasher slot fired more than 50 ns after the high charge channel and are + * more than 12 meters away from the high charge channel. */ + nhit = 0; + /* Number of PMTs which fired 50 ns after the high charge channel. */ + nhit_late = 0; + /* Number of PMTs which are 12 meters or more away from the high charge + * channel. */ + nhit_far = 0; + for (i = 0; i < MAX_PMTS; i++) { + /* Skip PMTs which weren't hit. */ + if (!ev->pmt_hits[i].hit) continue; - id = crate*64 + card*4 + channel/8; + /* Only count regular PMTs without any flags. */ + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; - /* Skip PMTs in the same paddle card as the high charge channel. */ - if (id == flasher_pc) continue; + crate = i/512; + card = (i % 512)/32; + channel = i % 32; - nhit += 1; + id = crate*64 + card*4 + channel/8; - /* Calculate the distance from the current channel to the high charge - * channel. */ - SUB(pmt_dir,pmts[i].pos,flasher_pos); - distance = NORM(pmt_dir); + /* Skip PMTs in the same card as the high charge channel. */ + if (id/4 == flasher_pc/4) continue; - /* If this channel is 12 meters or more away, increment nhit_far. */ - if (distance >= 1200.0) nhit_far += 1; + nhit += 1; - /* If this channel fired more than 50 ns after the high charge channel, - * increment nhit_late. */ - if (ev->pmt_hits[i].t > t + 50.0) nhit_late += 1; - } + /* Calculate the distance from the current channel to the high charge + * channel. */ + SUB(pmt_dir,pmts[i].pos,flasher_pos); + distance = NORM(pmt_dir); + + /* If this channel is 12 meters or more away, increment nhit_far. */ + if (distance >= 1200.0) nhit_far += 1; - /* If less than 70% of the regular PMTs fired within 50 ns of the high - * charge channel, then it's not a flasher. */ - if (nhit_late < nhit*0.7) return 0; + /* If this channel fired more than 50 ns after the high charge channel, + * increment nhit_late. */ + if (ev->pmt_hits[i].t > t + 50.0) nhit_late += 1; + } - /* If less than 70% of the regular PMTs were closer than 12 meters away - * from the high charge channel, then it's not a flasher. */ - if (nhit_far < nhit*0.7) return 0; + /* If at least 70% of the regular PMTs fired within 50 ns and were at + * least 12 meters away from the high charge channel, then it's a + * flasher! */ + if (nhit_late >= 7*nhit/10 && nhit_far >= 7*nhit/10) return 1; + } - return 1; + return 0; } @@ -20,9 +20,21 @@ #include "event.h" #include "zebra.h" +/* Minimum number of normal PMTs which must be hit to be tagged as an incoming + * muon. */ +#define MUON_MIN_NHIT 150 +/* Minimum number of OWL PMTs which must be hit to be tagged as an incoming + * muon. */ +#define MUON_MIN_OWL_NHIT 5 +/* Minimum value for the maximum OWL QHS value. */ +#define MUON_OWL_QHS 20.0 +/* Percentile of PSUP PMT times to define early OWL hits. */ +#define MUON_OWL_TIME_PERCENTILE 0.1 + /* QvNHIT ratio threshold. */ #define QRATIO_THRESHOLD 0.25 +int is_muon(event *ev); int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev); int crate_isotropy(event *ev); int qvnhit(event *ev); @@ -6172,7 +6172,7 @@ skip_mc: last_run = ev.run; } - rv = get_event(f,&ev,&b); + if (get_event(f,&ev,&b)) goto skip_event; nhit = get_nhit(&ev); @@ -6181,6 +6181,11 @@ skip_mc: fprintf(fout, " gtid: %i\n", ev.gtid); fprintf(fout, " nhit: %zu\n", nhit); fprintf(fout, " is_flasher: %i\n", is_flasher(&ev)); + fprintf(fout, " is_muon: %i\n", is_muon(&ev)); + fprintf(fout, " is_junk: %i\n", junk_cut(f, &bmast, &b)); + fprintf(fout, " is_crate_isotropy: %i\n", crate_isotropy(&ev)); + fprintf(fout, " is_qvnhit: %i\n", qvnhit(&ev)); + fprintf(fout, " is_neck_event: %i\n", is_neck_event(&ev)); } if (nhit < min_nhit) goto skip_event; diff --git a/src/zdab_utils.c b/src/zdab_utils.c index ba76b63..6dffe0f 100644 --- a/src/zdab_utils.c +++ b/src/zdab_utils.c @@ -78,6 +78,11 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev) channel = bpmt.pin % 32; id = crate*512 + card*32 + channel; + if (id >= MAX_PMTS) { + fprintf(stderr, "PMT hit from %i/%i/%i\n", crate, card, channel); + return -1; + } + if (ev->pmt_hits[id].hit) { fprintf(stderr, "%i/%i/%i is in the PMT bank twice!\n", crate, card, channel); } @@ -105,6 +110,9 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev) fprintf(stderr, "%i/%i/%i has PMT type %s but expected %s based on bank\n", crate, card, channel, pmt_type_string, pmt_names[i]); } + /* For now we assume SNOMAN has it correct. */ + pmts[id].pmt_type = pmt_types[i]; + if (!b.next) break; rv = zebra_get_bank(f,&b,b.next); |