From e89a28cb4e7e6cbbe97709d6703f502e9107977a Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 9 Sep 2019 12:32:55 -0500 Subject: add a first draft of a data cleaning cut to detect breakdown events --- src/dc.c | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) (limited to 'src/dc.c') diff --git a/src/dc.c b/src/dc.c index 5f2d4c6..6195eac 100644 --- a/src/dc.c +++ b/src/dc.c @@ -196,6 +196,8 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev) status |= DC_FTS; if (is_itc(ev)) status |= DC_ITC; + if (is_breakdown(ev)) + status |= DC_BREAKDOWN; return status; } @@ -604,6 +606,47 @@ int is_neck_event(event *ev) return 0; } +/* Returns 1 if the event is classified as a breakdown event. The definition + * for a breakdown event is: + * + * - nhit > 1000 + * + * - a single crate has at least 256 hits and the median TAC for that crate + * is more than 500 from the median TAC for all the other crates with at + * least 20 hits + * + */ +int is_breakdown(event *ev) +{ + int i; + double tac[20][512]; + int nhit[20] = {0}; + double median_tac[20]; + size_t index[20]; + + /* Breakdown event must have an nhit greater than 1000. */ + if (ev->nhit < 1000) return 0; + + for (i = 0; i < MAX_PMTS; i++) { + if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_NORMAL) continue; + + tac[i/512][nhit[i/512]++] = ev->pmt_hits[i].tac; + } + + for (i = 0; i < 19; i++) { + if (nhit[i] <= MIN_NHIT_CRATE) { + median_tac[i] = 0.0; + continue; + } + gsl_sort(&tac[i][0],1,nhit[i]); + median_tac[i] = gsl_stats_median_from_sorted_data(&tac[i][0],1,nhit[i]); + } + + fargsort(median_tac,19,index); + + return (nhit[index[18]] >= MIN_NHIT_BREAKDOWN && median_tac[index[18]] > median_tac[index[17]] + 500); +} + /* 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. -- cgit