diff options
Diffstat (limited to 'src/dc.c')
| -rw-r--r-- | src/dc.c | 201 |
1 files changed, 201 insertions, 0 deletions
diff --git a/src/dc.c b/src/dc.c new file mode 100644 index 0000000..897a645 --- /dev/null +++ b/src/dc.c @@ -0,0 +1,201 @@ +#include "event.h" +#include "sort.h" +#include <stddef.h> /* for size_t */ +#include "misc.h" +#include "pmt.h" +#include "vector.h" + +/* 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, + * with a few small changes. + * + * Here we define a flasher as an event with the following criteria: + * + * - 31 <= nhit <= 1000 + * + * - Look for the paddle card with the most PMTs hit. This paddle card should + * have 4 or more hits. + * + * - If such a paddle card exists, look for one channel in the card that hosts + * this paddle card which has either QHS or QHL values (uncalib'd) > 4000, + * whereas the QHS and QHL values for the other channels in the card are < + * 800. Note that QHS and QHL values below 300 are sent above 4096. + * + * Note: This is changed from the original definition. We look for all the + * other channels to have a QHL and QHS below 800 instead of 700 since the + * latter was failing to tag many obvious flashers. + * + * - At least 70% of the (regular) pmts in the event which are not in the + * flasher slot should have fired more than 50 ns after the high charge + * channel. + * + * - At least 70% of the (regular) pmts in the event which are not in the + * flasher slot should be at a distance >= 12.0 m from the high charge pmt. + * + * The one criteria we *don't* use that is on the webpage is the AMB cut: + * + * - AMB Int >= 200 + * + * since I can't seem to figure out where that is stored in the data structure. + */ +int is_flasher(event *ev) +{ + int hits_pc[1280] = {0}; + int qhs, qhl; + int qhl_pc[8] = {0}; + int qhs_pc[8] = {0}; + double t_pc[8] = {0}; + size_t crate, card, channel, id; + int i, j; + size_t index[1280], index_qhs[8], index_qhl[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; + + /* Flasher event must have an nhit between 31 and 1000. */ + if (ev->nhit < 31 || 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; + + crate = i/512; + card = (i % 512)/32; + channel = i % 32; + + id = crate*64 + card*4 + channel/8; + + hits_pc[id] += 1; + } + + argsort(hits_pc, LEN(hits_pc), index); + + /* The paddle card with the most hits must have >= 4 hits. */ + if (hits_pc[index[1279]] < 4) return 0; + + /* Loop over the most hit paddle cards in order from most hits to least + * hits. We do this because it's possible that more than one paddle card + * has the same number of channels hit, and so we need to test all of them. */ + 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; + + for (i = 0; i < MAX_PMTS; i++) { + if (!ev->pmt_hits[i].hit) continue; + + crate = i/512; + card = (i % 512)/32; + channel = i % 32; + + id = crate*64 + card*4 + channel/8; + + if (id != index[j]) continue; + + /* Get the uncalibrated QHS and QHL charges. */ + qhs = ev->pmt_hits[i].qihs; + qhl = ev->pmt_hits[i].qihl; + + if (qhs < 300) { + /* QHS values below 300 are set to 4095. */ + qhs_pc[channel % 8] = 4095; + } else { + qhs_pc[channel % 8] = qhs; + } + + if (qhl < 300) { + /* QHL values below 300 are set to 4095. */ + qhl_pc[channel % 8] = 4095; + } else { + qhl_pc[channel % 8] = qhl; + } + + t_pc[channel % 8] = ev->pmt_hits[i].t; + } + + /* 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; + } + } + } + + /* 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); + + /* 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; + + /* Only count regular PMTs without any flags. */ + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; + + crate = i/512; + card = (i % 512)/32; + channel = i % 32; + + id = crate*64 + card*4 + channel/8; + + /* Skip PMTs in the same paddle card as the high charge channel. */ + if (id == flasher_pc) continue; + + nhit += 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 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 fired within 50 ns of the high + * charge channel, then it's not a flasher. */ + if (nhit_late < nhit*0.7) return 0; + + /* 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; + + return 1; +} |