diff options
| -rw-r--r-- | src/dc.c | 82 | ||||
| -rw-r--r-- | src/dc.h | 7 |
2 files changed, 42 insertions, 47 deletions
@@ -222,22 +222,25 @@ uint32_t get_dc_word(event *ev, zebraFile *f, zebraBank *bmast, zebraBank *bev) * * - 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 + * - at least 1 OWL hit which has a high charge and is early relative to nearby + * normal PMTs * - * 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_ehs; + size_t i, j; + double t_nearby[MAX_PMTS], q_nearby[MAX_PMTS]; + int n_nearby; + size_t nhit_owl, nhit_owl_early; + double pmt_dir[3]; + double distance; + double t_median, q_median; + + /* Require at least 150 normal hits. */ + if (ev->nhit < MUON_MIN_NHIT) return 0; - nhit = 0; nhit_owl = 0; - max_owl_ehs = 0.0; + nhit_owl_early = 0; for (i = 0; i < MAX_PMTS; i++) { /* Require good calibrations. * @@ -251,47 +254,40 @@ int is_muon(event *ev) * 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; + if (!ev->pmt_hits[i].hit || pmts[i].pmt_type != PMT_OWL || ev->pmt_hits[i].ept <= -100) continue; - switch (pmts[i].pmt_type) { - case PMT_OWL: - if (nhit_owl == 0) { - max_owl_ehs = ev->pmt_hits[i].ehs; - } else if (ev->pmt_hits[i].ehs > max_owl_ehs) { - max_owl_ehs = ev->pmt_hits[i].ehs; - } - nhit_owl += 1; - break; - case PMT_NORMAL: - if (!ev->pmt_hits[i].flags) ept[nhit++] = ev->pmt_hits[i].ept; - break; - } - } + n_nearby = 0; + for (j = 0; j < MAX_PMTS; j++) { + if (!ev->pmt_hits[j].hit || pmts[j].pmt_type != PMT_NORMAL || ev->pmt_hits[j].ept <= -100) continue; - /* 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; + SUB(pmt_dir,pmts[i].pos,pmts[j].pos); + distance = NORM(pmt_dir); - /* Require at least 1 OWL tube to have a QHS >= 20.0. */ - if (max_owl_ehs < MUON_OWL_QHS) return 0; + if (distance < MUON_OWL_NEARBY_DISTANCE) { + t_nearby[n_nearby] = ev->pmt_hits[j].ept; + q_nearby[n_nearby] = ev->pmt_hits[j].ehs; + n_nearby += 1; + } + } - /* Sort the times. */ - gsl_sort(ept,1,nhit); + /* Not sure what to do if there are no nearby PMTs. */ + if (n_nearby == 0) continue; - /* Calculate the time to define early OWL hits. */ - early_time = gsl_stats_quantile_from_sorted_data(ept,1,nhit,MUON_OWL_TIME_PERCENTILE); + gsl_sort(t_nearby,1,n_nearby); + gsl_sort(q_nearby,1,n_nearby); + t_median = gsl_stats_median_from_sorted_data(t_nearby,1,n_nearby); + q_median = gsl_stats_median_from_sorted_data(q_nearby,1,n_nearby); - 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 (ev->pmt_hits[i].ehs/q_median > 1.0 && ev->pmt_hits[i].ept < t_median) + nhit_owl_early += 1; - if (pmts[i].pmt_type == PMT_OWL && ev->pmt_hits[i].ept < early_time) nhit_owl_early += 1; + nhit_owl += 1; } - return nhit_owl_early >= 2; + /* Require at least 5 OWL hits. */ + if (nhit_owl < MUON_MIN_OWL_NHIT) return 0; + + return nhit_owl_early >= 1; } /* Returns 1 if the event fails the junk cut. The definition for the junk cut @@ -71,10 +71,9 @@ /* 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 +/* Maximum distance between OWL and normal PMTs when computing the median time + * and charge. */ +#define MUON_OWL_NEARBY_DISTANCE 300.0 /* QvNHIT ratio threshold. */ #define QRATIO_THRESHOLD 0.25 |