path: root/src/zdab_utils.c

/* Copyright (c) 2019, Anthony Latorre <tlatorre at uchicago>
 *
 * This program is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation, either version 3 of the License, or (at your option)
 * any later version.

 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details.

 * You should have received a copy of the GNU General Public License along with
 * this program. If not, see <https://www.gnu.org/licenses/>.
 */

#include "Record_Info.h"
#include <stdint.h>
#include "zdab_utils.h"
#include "pack2b.h"
#include <stdlib.h> /* for size_t */
#include <stdio.h> /* for fprintf() */
#include "event.h"
#include "zebra.h"
#include "misc.h"

size_t get_nhit(event *ev)
{
    /* Returns the number of PMT hits in event `ev`.
     *
     * Note: Only hits on normal PMTs which aren't flagged are counted. */
    size_t i, nhit;

    nhit = 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) continue;

        nhit++;
    }

    return nhit;
}

int get_event(zebraFile *f, event *ev, zebraBank *bev)
{
    /* Read all the PMT banks from the zebra file and update `ev`.
     *
     * Returns 0 on success, -1 on error. */
    int i, rv;
    PMTBank bpmt;
    zebraBank b;
    int id, crate, card, channel;
    static int pmt_links[] = {KEV_PMT,KEV_OWL,KEV_LG,KEV_FECD,KEV_BUTT,KEV_NECK};
    static int pmt_types[] = {PMT_NORMAL,PMT_OWL,PMT_LG,PMT_CALIBRATION,PMT_BUTT,PMT_NECK};
    static char *pmt_names[] = {"PMT","OWL","LG","FECD","BUTT","NECK"};
    char pmt_type_string[256];

    for (i = 0; i < MAX_PMTS; i++) {
        ev->pmt_hits[i].hit = 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);
            card = bpmt.pin/1024;
            crate = (bpmt.pin % 1024)/32;
            channel = bpmt.pin % 32;
            id = crate*512 + card*32 + channel;

            if (id >= MAX_PMTS) {
                fprintf(stderr, "PMT hit from %i/%i/%i\n", crate, card, channel);
                return -1;
            }

            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;
            ev->pmt_hits[id].pt1 = bpmt.pt1;

            /* Clear the PMT_FLAG_DIS bit. */
            ev->pmt_hits[id].flags &= ~PMT_FLAG_DIS;

            /* Make sure we have good calibrations. Technically, we should only
             * have to check the KPF_DIS bit since according to the SNOMAN
             * docs:
             *
             *     The discard bit KPF_DIS of the 1-bit flags KPMT_PF is also
             *     set if any other failure bit, such as KPF_NO_CAL or
             *     KPF_BAD_CAL is set. So only this bit need be tested for PMT
             *     rejection.
             *
             * But this isn't true! In run 10,000 GTID 140730 channel 11/6/17
             * has the KPF_BAD_CAL bit set but not the KPF_DIS bit set. I also
             * looked at the SNOMAN code and confirmed that this isn't the
             * case. */
            if (bpmt.pf & (KPF_DIS | KPF_NO_CAL | KPF_BAD_CAL))
                ev->pmt_hits[id].flags |= PMT_FLAG_DIS;

            if (bpmt.qms)
                ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;

            /* Currently, the charge model only deals with QHS, so we flag any hits which have a bad or railed QHS value.
             *
             * FIXME: In the future, it would be nice to use the best charge
             * word (either QHS or QLX) depending on the if QHS is railed or
             * not, but I need to do more work to see how the QLX values are
             * normalized and if the existing charge model is good enough. */
            if (bpmt.pihs >= 4095 || bpmt.pihs < 300)
                ev->pmt_hits[id].flags |= PMT_FLAG_CHARGE;

            if (pmts[id].pmt_type != pmt_types[i]) {
                get_pmt_type_string(pmts[id].pmt_type,pmt_type_string);
                fprintf(stderr, "%i/%i/%i has PMT type %s but expected %s based on bank\n", crate, card, channel, pmt_type_string, pmt_names[i]);
            }

            /* The above check prints out a warning on a handful of channels:
             *
             *     3/15/9 has PMT type NECK but expected OWL based on bank
             *     18/0/14 has PMT type INVALID but expected BUTT based on bank
             *     (+ a handful of other INVALID/BUTT discrepancies in crates 18 and 15)
             *
             * I sent an email to the SNO mailing list and Stan responded:
             *
             * > That PMT ID jogged a neuron somewhere, and I went hunting
             * > through my email archive.  Here is an email I sent Neil 14 years
             * > ago (?!) after we added the new neck tubes which might explain
             * > how a neck tube ended up classified as an OWL.  When new neck
             * > tubes were added for the NCD phase, they broke the range-based
             * > logical PMT numbering:
             *
             * > On Tue, 9 Nov 2004 at 09:40:06 AM,  Neil wrote:
             * > > Hi,
             * > >
             * > >     A new file has been released as should be in the database that has
             * > > the updated neck tube positions. This file will also be in the new
             * > > version of snoman to be released today or tomorrow.
             * > >
             * > >    Neil.
             * > >
             * > > Stan Seibert wrote:
             * > > >
             * > > > Hi Neil,
             * > > >
             * > > > I've noticed several issues in how neck tubes are handled in SNOMAN.
             * > > >
             * > > > * The 3 new tubes are plugged into channels 18/15/0, 18/15/10, and
             * > > > 3/15/9.  According to map_ccc_tube.dat, the last two are tube numbers
             * > > > 9635 and 9655.  The first, 18/15/0, currently has the entry -999999,
             * > > > indicating it is a spare.  Is there a more up-to-date version of
             * > > > map_ccc_tube.dat I should be looking at?  I tried reading through
             * > > > Phil's pmt database for xsnoed, but was unable to find anything that
             * > > > resembled tube numbers as SNOMAN defines them.
             * > > >
             * > > >
             * > > > * Once that is corrected, I will need to fix ccc_type.for which
             * > > > decides which PMT list (normal, owl, neck, etc.) each PMT hit goes
             * > > > into.  This will uglify the function a bit because the neck tubes are
             * > > > no longer in a contiguous block of tube numbers.  Numbers 9635 and
             * > > > 9655 are actually in the OWL range.  How should backward compatibility
             * > > > be handled here?  Do I hardcode a run number and check for it?
             * > > >
             * > > >
             * > > > * With that fix, upk_unpacker_zdab_pmt.for will then file the neck
             * > > > tube in its appropriate list, right?
             * > > >
             * > > >
             * > > > Please let me know if this makes sense and where I can get updated
             * > > > information.   Thanks.
             * > > >
             * > > > ---
             * > > > Stan Seibert
             *
             * I also looked at the run in XSNOED and it looks like 3/15/9
             * should be an OWL PMT and not a neck PMT.
             *
             * Based on the email and the fact that 3/15/9 *should* be an OWL
             * PMT, I am going to assume that the database that SNOMAN uses to
             * figure out the PMT types is correct and so reassign the PMT type
             * based on that. */
            pmts[id].pmt_type = pmt_types[i];

            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;
            }
        }
    }

    ev->nhit = get_nhit(ev);

    return 0;
}

int isOrphan(aPmtEventRecord *pmtRecord)
{
    /* Returns non-zero if the specified event is an orphan. */
    int i;

    uint32_t *mtc_data = (uint32_t *) &pmtRecord->TriggerCardData;

    for (i = 0; i < 6; ++i) {
        if (*mtc_data != 0) return 0;
        ++mtc_data;
    }

    return 1;
}

// PH 04/23/98
// Swap 4-byte integer/floats between native and external format
void swap_int32(int32_t *val_pt, int count)
{
    int32_t *last = val_pt + count;
    while (val_pt < last) {
        *val_pt = ((*val_pt << 24) & 0xff000000) |
                  ((*val_pt <<  8) & 0x00ff0000) |
                  ((*val_pt >>  8) & 0x0000ff00) |
                  ((*val_pt >> 24) & 0x000000ff);
        ++val_pt;
    }
    return;
}

// Swap 2-byte integers between native and external format
void swap_int16(int16_t *val_pt, int count)
{
    char tmp;
    int i;
    for (i=0; i<count; ++i) {
        tmp = ((char *) val_pt)[0];
        ((char *) val_pt)[0] = ((char *) val_pt)[1];
        ((char *) val_pt)[1] = tmp;
        ++val_pt;
    }
    return;
}

void unpack_mcgn(uint32_t *data, MCGNBank *b)
{
    unpack((uint8_t *) data,      "l",&b->id);
    unpack((uint8_t *) (data+1),  "l",&b->num);
    unpack((uint8_t *) (data+2),  "l",&b->radcor_proc);
    unpack((uint8_t *) (data+3),  "l",&b->radcor_made_gamma);
    unpack((uint8_t *) (data+4),  "l",&b->spare5);
    unpack((uint8_t *) (data+5),  "l",&b->spare6);
    unpack((uint8_t *) (data+7),  "l",&b->spare7);
    unpack((uint8_t *) (data+8),  "l",&b->spare8);
    unpack((uint8_t *) (data+9),  "l",&b->spare9);
    unpack((uint8_t *) (data+10), "l",&b->spare10);
    unpack((uint8_t *) (data+12), "f",&b->radcor_xtot);
    unpack((uint8_t *) (data+13), "f",&b->radcor_xdif);
    unpack((uint8_t *) (data+14), "f",&b->gentim_xpar);
    unpack((uint8_t *) (data+15), "f",&b->spare14);
    unpack((uint8_t *) (data+16), "f",&b->spare15);
    unpack((uint8_t *) (data+17), "f",&b->spare16);
    unpack((uint8_t *) (data+18), "f",&b->spare17);
    unpack((uint8_t *) (data+19), "f",&b->spare18);
    unpack((uint8_t *) (data+20), "f",&b->spare19);
    unpack((uint8_t *) (data+21), "f",&b->spare20);
}

void unpack_mcvx(uint32_t *data, MCVXBank *b)
{
    unpack((uint8_t *) data,      "l",&b->cls);
    unpack((uint8_t *) (data+1),  "l",&b->inc);
    unpack((uint8_t *) (data+2),  "f",&b->x);
    unpack((uint8_t *) (data+3),  "f",&b->y);
    unpack((uint8_t *) (data+4),  "f",&b->z);
    unpack((uint8_t *) (data+5),  "F",&b->tim);
    unpack((uint8_t *) (data+7),  "l",&b->rgn);
    unpack((uint8_t *) (data+8),  "l",&b->idm);
    unpack((uint8_t *) (data+9),  "l",&b->rg2);
    unpack((uint8_t *) (data+10), "l",&b->im2);
    unpack((uint8_t *) (data+12), "f",&b->bnx);
    unpack((uint8_t *) (data+13), "f",&b->bny);
    unpack((uint8_t *) (data+14), "f",&b->bnz);
    unpack((uint8_t *) (data+15), "l",&b->cer);
}

void unpack_mctk(uint32_t *data, MCTKBank *b)
{
    unpack((uint8_t *) data,      "l",&b->idp);
    unpack((uint8_t *) (data+1),  "f",&b->drx);
    unpack((uint8_t *) (data+2),  "f",&b->dry);
    unpack((uint8_t *) (data+3),  "f",&b->drz);
    unpack((uint8_t *) (data+4),  "f",&b->ene);
    unpack((uint8_t *) (data+5),  "l",&b->rgn);
    unpack((uint8_t *) (data+6),  "l",&b->idm);
    unpack((uint8_t *) (data+7),  "f",&b->plx);
    unpack((uint8_t *) (data+8),  "f",&b->ply);
    unpack((uint8_t *) (data+9),  "f",&b->plz);
    unpack((uint8_t *) (data+10), "f",&b->stp);
    unpack((uint8_t *) (data+11), "f",&b->near);
}

void unpack_ev(uint32_t *data, EVBank *b)
{
    unpack((uint8_t *) data,     "l",&b->run);
    unpack((uint8_t *) (data+1), "l",&b->evn);
    unpack((uint8_t *) (data+2), "l",&b->dtp);
    unpack((uint8_t *) (data+3), "l",&b->jdy);
    unpack((uint8_t *) (data+4), "l",&b->ut1);
    unpack((uint8_t *) (data+5), "l",&b->ut2);
    unpack((uint8_t *) (data+6), "l",&b->dte);
    unpack((uint8_t *) (data+7), "l",&b->hmsc);
    unpack((uint8_t *) (data+8), "F",&b->gtr);
    unpack((uint8_t *) (data+10),"l",&b->npm);
    unpack((uint8_t *) (data+11),"f",&b->nph);
    unpack((uint8_t *) (data+12),"l",&b->sub_run);
    unpack((uint8_t *) (data+13),"l",&b->mc_pck);
    unpack((uint8_t *) (data+14),"l",&b->rec);
    unpack((uint8_t *) (data+15),"l",&b->vpck);
    unpack((uint8_t *) (data+16),"l",&b->gtr_id);
    unpack((uint8_t *) (data+17),"l",&b->trg_type);
    unpack((uint8_t *) (data+18),"l",&b->peak);
    unpack((uint8_t *) (data+19),"l",&b->diff);
    unpack((uint8_t *) (data+20),"l",&b->integral);
    unpack((uint8_t *) (data+21),"l",&b->err);
    unpack((uint8_t *) (data+22),"l",&b->data_set);
    unpack((uint8_t *) (data+22),"lll",&b->spare1[0],
                                       &b->spare1[1],
                                       &b->spare1[2]);
    unpack((uint8_t *) (data+26),"l",&b->ncd_status);
    unpack((uint8_t *) (data+27),"l",&b->num_muxg);
    unpack((uint8_t *) (data+29),"l",&b->num_mux);
    unpack((uint8_t *) (data+29),"l",&b->num_scope);
    unpack((uint8_t *) (data+30),"lllll",&b->spare2[0],
                                         &b->spare2[1],
                                         &b->spare2[2],
                                         &b->spare2[3],
                                         &b->spare2[4]);
    unpack((uint8_t *) (data+35),"l",&b->ncd_clk_up);
    unpack((uint8_t *) (data+36),"l",&b->ncd_clk_lw);
    unpack((uint8_t *) (data+37),"l",&b->ncd_reg);
    unpack((uint8_t *) (data+38),"l",&b->ncd_gtid);
    unpack((uint8_t *) (data+39),"l",&b->ncd_sync);
    unpack((uint8_t *) (data+40),"l",&b->spare3[0],
                                     &b->spare3[1],
                                     &b->spare3[2],
                                     &b->spare3[3],
                                     &b->spare3[4],
                                     &b->spare3[5],
                                     &b->spare3[6],
                                     &b->spare3[7],
                                     &b->spare3[8],
                                     &b->spare3[9]);
}

void unpack_pmt(uint32_t *data, PMTBank *b)
{
    unpack((uint8_t *) data,"l",&b->pn);
    unpack((uint8_t *) (data+1),"l",&b->pf);
    unpack((uint8_t *) (data+2),"f",&b->pt);
    unpack((uint8_t *) (data+3),"f",&b->phl);
    unpack((uint8_t *) (data+4),"f",&b->phs);
    unpack((uint8_t *) (data+5),"f",&b->plx);
    unpack((uint8_t *) (data+6),"f",&b->pt0);
    unpack((uint8_t *) (data+7),"l",&b->pif);
    unpack((uint8_t *) (data+8),"f",&b->pit);
    unpack((uint8_t *) (data+9),"f",&b->pihl);
    unpack((uint8_t *) (data+10),"f",&b->pihs);
    unpack((uint8_t *) (data+11),"f",&b->pilx);
    unpack((uint8_t *) (data+12),"f",&b->pit0);
    unpack((uint8_t *) (data+13),"l",&b->cell);
    unpack((uint8_t *) (data+14),"l",&b->pin);
    unpack((uint8_t *) (data+15),"f",&b->tslh);
    unpack((uint8_t *) (data+16),"f",&b->hca);
    unpack((uint8_t *) (data+17),"l",&b->eca_val);
    unpack((uint8_t *) (data+18),"l",&b->pca_val);
    unpack((uint8_t *) (data+19),"l",&b->anxx);
    unpack((uint8_t *) (data+20),"f",&b->ept);
    unpack((uint8_t *) (data+21),"f",&b->ehl);
    unpack((uint8_t *) (data+22),"f",&b->ehs);
    unpack((uint8_t *) (data+23),"f",&b->elx);
    unpack((uint8_t *) (data+24),"f",&b->pt1);
    unpack((uint8_t *) (data+25),"f",&b->ptm);
    unpack((uint8_t *) (data+26),"f",&b->ptms);
    unpack((uint8_t *) (data+27),"f",&b->qm);
    unpack((uint8_t *) (data+28),"l",&b->qms);
    unpack((uint8_t *) (data+29),"f",&b->qrc);
}

int swap_PmtRecord(aPmtEventRecord *aPmtRecord, size_t size)
{
    /* Swap a Pmt Event Record. This function swaps both the Pmt event record
     * and the PMT hits and sub fields. Returns -1 if the PMT record has too
     * many hits. */
    SWAP_INT32(aPmtRecord, sizeof(aPmtEventRecord)/sizeof(uint32_t));

    int npmt = aPmtRecord->NPmtHit;

    if (npmt > MAX_NHIT) {
        fprintf(stderr, "Read error: Bad ZDAB -- %d pmt hit!", npmt);
        return -1;
    } else {
        if (size < sizeof(aPmtEventRecord) + 3*npmt*4) {
            fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, but there are %i PMT hits", size, npmt);
            return -1;
        }
        // swap the hit data
        SWAP_INT32(aPmtRecord + 1, 3*npmt);
        // swap the sub-fields
        uint32_t *sub_header = &aPmtRecord->CalPckType;
        while (*sub_header & SUB_NOT_LAST) {
            if (size < (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4 + 4) {
                fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, "
                    "but sub-field requires %lu bytes",
                    size,
                    (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4 + 4);
                return -1;
            }
            sub_header += (*sub_header & SUB_LENGTH_MASK);
            SWAP_INT32(sub_header, 1); // swap the sub-field header
            // get number of data words (-1 because we don't want to include header size)
            uint32_t data_words = (*sub_header & SUB_LENGTH_MASK) - 1;
            if (size < (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4) {
                fprintf(stderr, "swap_PmtRecord: size of record is %zu bytes, "
                    "but sub-field requires %lu bytes",
                    size,
                    (sub_header - (uint32_t *) aPmtRecord)*4 + (*sub_header & SUB_LENGTH_MASK)*4);
                return -1;
            }
            SWAP_INT32(sub_header+1, data_words);
        }
    }

    return 0;
}

void swap_TrigRecord(struct TriggerInfo *aTrigRecord)
{
    /* Byte swap a Trigger Record. */
    SWAP_INT32(aTrigRecord, sizeof(struct TriggerInfo)/sizeof(uint32_t));
}

void swap_RunRecord(struct RunRecord *aRunRecord)
{
    /* Byte swap a Run Record. */
    SWAP_INT32(aRunRecord, sizeof(struct RunRecord)/sizeof(uint32_t));
}