diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/dc.c | 259 | ||||
| -rw-r--r-- | src/dc.h | 8 | ||||
| -rw-r--r-- | src/event.h | 5 | ||||
| -rw-r--r-- | src/zdab_utils.c | 10 | ||||
| -rw-r--r-- | src/zdab_utils.h | 17 |
5 files changed, 299 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 @@ -18,7 +18,15 @@ #define DC_H #include "event.h" +#include "zebra.h" +/* QvNHIT ratio threshold. */ +#define QRATIO_THRESHOLD 0.25 + +int junk_cut(zebraFile *f, zebraBank *bmast, zebraBank *bev); +int crate_isotropy(event *ev); +int qvnhit(event *ev); +int is_neck_event(event *ev); int is_flasher(event *ev); #endif diff --git a/src/event.h b/src/event.h index a561a1f..0f1efce 100644 --- a/src/event.h +++ b/src/event.h @@ -28,7 +28,11 @@ typedef struct pmt_hit { int hit; + float ept; float t; + float ehl; + float ehs; + float elx; float qhl; float qhs; float qlx; @@ -36,6 +40,7 @@ typedef struct pmt_hit { uint16_t qihs; uint16_t qilx; int flags; + int pf; } pmt_hit; typedef struct event { diff --git a/src/zdab_utils.c b/src/zdab_utils.c index f0651d5..ba76b63 100644 --- a/src/zdab_utils.c +++ b/src/zdab_utils.c @@ -77,14 +77,24 @@ int get_event(zebraFile *f, event *ev, zebraBank *bev) crate = (bpmt.pin % 1024)/32; channel = bpmt.pin % 32; id = crate*512 + card*32 + channel; + + if (ev->pmt_hits[id].hit) { + fprintf(stderr, "%i/%i/%i is in the PMT bank twice!\n", crate, card, channel); + } + ev->pmt_hits[id].hit = 1; + ev->pmt_hits[id].ept = bpmt.ept; ev->pmt_hits[id].t = bpmt.pt; ev->pmt_hits[id].qihl = bpmt.pihl; ev->pmt_hits[id].qihs = bpmt.pihs; ev->pmt_hits[id].qilx = bpmt.pilx; + ev->pmt_hits[id].ehl = bpmt.ehl; + ev->pmt_hits[id].ehs = bpmt.ehs; + ev->pmt_hits[id].elx = bpmt.elx; ev->pmt_hits[id].qhl = bpmt.phl; ev->pmt_hits[id].qhs = bpmt.phs; ev->pmt_hits[id].qlx = bpmt.plx; + ev->pmt_hits[id].pf = bpmt.pf; /* Clear the PMT_FLAG_DIS bit. */ ev->pmt_hits[id].flags &= ~PMT_FLAG_DIS; if (bpmt.pf & KPF_DIS) diff --git a/src/zdab_utils.h b/src/zdab_utils.h index a2804e1..1c138b9 100644 --- a/src/zdab_utils.h +++ b/src/zdab_utils.h @@ -50,6 +50,12 @@ #define KMC_MCVX 2 #define KMC_MCGN 1 +/* Status bits for the EV bank. */ +#define KEV_ORPHAN 0 /* Event contains orphan data. */ +#define KEV_TSLH 1 /* Event contains TSLH information. */ +#define KEV_HCA 2 /* Event contains HCA information. */ +#define KEV_MC_GE_ERR 3 /* Event produced from MC that had geometry errors suppressed. */ + /* Links for the EV bank. */ #define KEV_NPFA 32 #define KEV_NECL 31 @@ -84,6 +90,17 @@ #define KEV_PT 2 #define KEV_FT 1 +/* Status bits for the MAST bank. */ + +#define KMAST_PMT 1 /* PMT data. */ +#define KMAST_PERM 2 /* Permanent data of any kind e.g. Run Header. */ +#define KMAST_NCD 3 /* NCD data. */ +#define KMAST_CMA 4 /* CMA data. */ +#define KMAST_ECA 5 /* ECA data. */ +#define KMAST_SUBRUN 6 /* Sub-run data, i.e. permanent bank data that is duplicated + at the start of the second and subsequent sub runs to + ensure the context is maintained across file boundaries. */ + /* Links for the MAST bank. */ #define KMAST_RLOG 25 /* Runlog data. */ #define KMAST_RLAI 24 /* Runlog analysis data. */ |