diff options
Diffstat (limited to 'src/dc.c')
| -rw-r--r-- | src/dc.c | 259 |
1 files changed, 259 insertions, 0 deletions
@@ -14,12 +14,271 @@ * this program. If not, see <https://www.gnu.org/licenses/>. */ +#include "dc.h" #include "event.h" #include "sort.h" #include <stddef.h> /* for size_t */ #include "misc.h" #include "pmt.h" #include "vector.h" +#include <gsl/gsl_sort.h> +#include "zdab_utils.h" +#include "zebra.h" + +/* Returns 1 if the event fails the junk cut. The definition for the junk cut + * comes from the SNOMAN companion: + * + * This will remove orphans, ECA data and events containing the same PMT + * more than once. + * + * 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 ids[MAX_PMTS]; + int nhit; + static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK}; + static char *pmt_names[] = {"PMT","OWL","LG","FECD","BUTT","NECK"}; + size_t index[MAX_PMTS]; + PMTBank bpmt; + zebraBank b; + int rv; + + /* Fail orphans. */ + if (bev->status & KEV_ORPHAN) return 1; + + /* Fail any events with the ECA bit set. */ + if (bmast->status & KMAST_ECA) return 1; + + /* Next we check to see if any PMTs got hit more than once. To do this we + * construct an array of the PMT PIN numbers, sort it, and then check for + * duplicates. */ + + nhit = 0; + for (i = 0; i < LEN(pmt_links); i++) { + if (bev->links[pmt_links[i]-1] == 0) continue; + + rv = zebra_get_bank(f,&b,bev->links[pmt_links[i]-1]); + + if (rv) { + fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err); + return -1; + } + + while (1) { + unpack_pmt(b.data, &bpmt); + ids[nhit++] = bpmt.pin; + if (!b.next) break; + + rv = zebra_get_bank(f,&b,b.next); + + if (rv) { + fprintf(stderr, "error getting %s bank: %s\n", pmt_names[i], zebra_err); + return -1; + } + } + } + + argsort(ids,nhit,index); + + for (i = 0; i < nhit - 1; i++) { + /* Check if we have duplicate PMT hits. */ + if (ids[i] == ids[i+1]) return 1; + } + + return 0; +} + +/* Returns 1 if the event fails the crate isotropy test. The definition of + * this cut comes from the SNOMAN companion: + * + * This tags events where more than 70% of the hits occur on one crate and + * more than 80% of those hits occur on two adjacent cards, including a + * wrap around to count cards 0 and 15 together. + * + */ +int crate_isotropy(event *ev) +{ + size_t i; + int crate_nhit[20]; + int card_nhit[16]; + int nhit; + int crate, card; + int max_crate; + + for (i = 0; i < LEN(crate_nhit); i++) crate_nhit[i] = 0; + + /* Loop over all the hits and add one to the crate_nhit array for each hit + * in a crate. */ + + nhit = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (!ev->pmt_hits[i].hit) continue; + + crate = i/512; + + crate_nhit[crate] += 1; + + nhit += 1; + } + + /* Check if any crate has more than 70% of the hits. */ + max_crate = -1; + for (i = 0; i < LEN(crate_nhit); i++) { + if (crate_nhit[i] > 7*nhit/10) max_crate = i; + } + + /* If no crates had more than 70% of the hits, return 0. */ + if (max_crate == -1) return 0; + + for (i = 0; i < LEN(card_nhit); i++) card_nhit[i] = 0; + + /* Now we count up the number of hits in adjacent cards. We store the + * number of hits in the card_nhit array as follows: Hits from card 15 or + * card 0 will be added to card_nhit[0], hits from cards 0 or 1 will go + * into card_nhit[1], etc. */ + + nhit = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (!ev->pmt_hits[i].hit) continue; + + crate = i/512; + + if (crate != max_crate) continue; + + card = (i % 512)/32; + + card_nhit[card] += 1; + card_nhit[(card + 1) % LEN(card_nhit)] += 1; + + nhit += 1; + } + + /* Check to see if 80% of the hits in this crate came from two adjacent + * cards. */ + for (i = 0; i < LEN(card_nhit); i++) { + if (card_nhit[i] > 8*nhit/10) return 1; + } + + return 0; +} + +/* Returns 1 if the event is tagged by the QvNHIT cut. The definition of + * this cut comes from the SNOMAN companion: + * + * Using pedestal subtracted charge and a hardwired gain of 32.3 first + * throw away the 10% of the tubes that have the largest charge and then + * divide this by 0.9* NHIT. Tag the event if this ratio is below 0.25. + * + * I copied the logic from the SNOMAN file flt_q_nhit_cut.for. */ +int qvnhit(event *ev) +{ + size_t i; + double ehl[MAX_PMTS]; + int nhit, max; + double sum, qhl_ratio; + + nhit = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + + /* Require good calibrations. */ + if (ev->pmt_hits[i].pf & (KPF_NO_CAL | KPF_BAD_CAL)) continue; + + /* Remove unphysical charge hits. */ + if (ev->pmt_hits[i].qhl <= -100) continue; + + /* FIXME: SNOMAN code checks for pmt terminator. I don't know where + * that info is stored. */ + + /* Note: We use QHL instead of QHS here because that's how SNOMAN did + * it. In the SNOMAN code there is a comment which says: + * + * we use QHL since it is a better measure of pickup (long integrate) + * + */ + if (ev->pmt_hits[i].hit) { + ehl[nhit++] = ev->pmt_hits[i].ehl/32.3; + } + } + + gsl_sort(ehl,1,nhit); + + max = nhit*9/10; + + /* Check that we have enough calibrated tubes. */ + if (max < 5) return 0; + + sum = 0.0; + for (i = 0; i < max; i++) { + sum += ehl[i]; + } + + qhl_ratio = sum/max; + + return qhl_ratio < QRATIO_THRESHOLD; +} + +/* Returns 1 if the event is classified as a neck event. The definition of + * a neck event comes from the SNOMAN companion: + * + * This cuts events containing neck tubes. It requires that either both + * tubes in the neck fire, or that one of those tubes fires and it has a + * high charge and is early. High charge is defined by a neck tube having a + * pedestal subtracted charge greater than 70 or less than -110. Early if + * defined by the neck tube having an ECA time 70ns or more before the + * average ECA time of the PSUP PMTS with z les than 0. After the cable + * changes to the neck tubes this time difference changes to 15ns. + * + */ +int is_neck_event(event *ev) +{ + size_t i; + int n; + int high_charge; + double avg_ept_psup; + int nhit; + + /* First, we compute the average ECA time for all normal PMTs. */ + + 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) { + avg_ept_psup += ev->pmt_hits[i].ept; + nhit += 1; + } + } + + if (nhit < 1) return 0; + + avg_ept_psup /= nhit; + + /* Now, we check if two or more neck tubes fired *or* a single neck tube + * fired with a high charge and 70 ns before the average time of the normal + * PMTs. */ + n = 0; + high_charge = 0; + for (i = 0; i < MAX_PMTS; i++) { + if (pmts[i].pmt_type == PMT_NECK && ev->pmt_hits[i].hit) { + n += 1; + + /* FIXME: The SNOMAN companion says "Early if defined by the neck + * tube having an ECA time 70ns or more before the average ECA time + * of the PSUP PMTS with z les than 0." When does this change take + * place? */ + if ((ev->pmt_hits[i].ehs > 70 || ev->pmt_hits[i].ehs < -110) && + (ev->pmt_hits[i].ept < avg_ept_psup - 70.0)) high_charge = 1; + } + } + + if (n >= 2 || high_charge) return 1; + + return 0; +} /* Returns 1 if the event is a flasher and 0 if it isn't. The definition of * a flasher event comes from |