path: root/src/zdab-cat.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 <stdio.h>
#include "zebra.h"
#include "event.h"
#include "zdab_utils.h"
#include "pmt.h"
#include "db.h"
#include "dqxx.h"
#include <inttypes.h> /* for PRIu32 macro */
#include <string.h> /* for memcpy() */
#include <errno.h> /* for errno */
#include "release.h"
#include "dc.h"
#include "hdf5.h"
#include <math.h>
#include "hdf5_utils.h"

/* Maximum number of events, primary seed particle tracks, and fits.
 *
 * FIXME: I should write out the events, tracks, and fits as they come up
 * instead of storing them until the end. */
#define MAX_NEVENTS 1000000

#define EV_RECORD      0x45562020 // 'EV  '
#define MCTK_RECORD    0x4d43544b // 'MCTK'
#define MCVX_RECORD    0x4d435658 // 'MCVX'

char *GitSHA1(void);
char *GitDirty(void);

void usage(void)
{
    fprintf(stderr,"Usage: ./zdab-cat [options] FILENAME\n");
    fprintf(stderr,"  -o      output file (default: stdout)\n");
    fprintf(stderr,"  --skip-second-event  only fit the first event after a MAST bank\n");
    fprintf(stderr,"  -v                   verbose\n");
    fprintf(stderr,"  -h      display this help message\n");
    exit(1);
}

int main(int argc, char **argv)
{
    int i;
    zebraFile *f;
    zebraBank bmast, mc, bmcgn, mctk, b;
    int rv;
    EVBank bev;
    FTPVBank bftpv;
    FTXKBank bftxk;
    RSPBank bftxr;
    MCBank bmc;
    MCTKBank bmctk;
    MCVXBank bmcvx;
    event ev = {0};
    char *filename = NULL;
    char *output = NULL;
    int skip_second_event = 0;
    int last_run;
    char dqxx_file[256];
    int nevents = 0;
    static HDF5Event hdf5_events[MAX_NEVENTS];
    int nmcgn = 0;
    static HDF5MCGN hdf5_mcgn[MAX_NEVENTS];
    int verbose = 0;

    for (i = 1; i < argc; i++) {
        if (strlen(argv[i]) >= 2 && !strncmp(argv[i], "--", 2)) {
            if (!strcmp(argv[i]+2,"skip-second-event")) {
                skip_second_event = 1;
                continue;
            }
        } else if (argv[i][0] == '-') {
            switch (argv[i][1]) {
            case 'o':
                output = argv[++i];
                break;
            case 'v':
                verbose = 1;
                break;
            case 'h':
                usage();
            default:
                fprintf(stderr, "unrecognized option '%s'\n", argv[i]);
                exit(1);
            }
        } else {
            filename = argv[i];
        }
    }

    if (!filename)
        usage();

    f = zebra_open(filename);

    if (!f) {
        fprintf(stderr, "%s\n", zebra_err);
        return 1;
    }

    if (load_pmt_info()) {
        zebra_close(f);
        return 1;
    }

    for (i = 0; i < MAX_PMTS; i++) {
        ev.pmt_hits[i].hit = 0;
        ev.pmt_hits[i].flags = 0;
    }

    ev.run = -1;

    dict *db = db_init();

    last_run = 10000;

    if (load_file(db, "DQXX_0000010000.dat", 0)) {
        fprintf(stderr, "failed to load DQXX_0000010000.dat: %s\n", db_err);
        goto err;
    }

    if (dqxx_init(db, &ev)) {
        fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err);
        goto err;
    }

    while (1) {
        rv = zebra_read_next_logical_record(f);

        if (rv == -1) {
            fprintf(stderr, "error getting logical record: %s\n", zebra_err);
            goto err;
        } else if (rv == 1) {
            /* EOF */
            break;
        }

        rv = zebra_get_bank(f, &bmast, f->first_bank);

        if (rv) {
            fprintf(stderr, "error getting MAST bank: %s\n", zebra_err);
            goto err;
        }

        if (bmast.links[KMAST_EV-1] == 0) {
            /* First logical record in SNOCR files doesn't have an EV bank. */
            continue;
        }

        if (bmast.links[KMAST_MC-1] == 0) goto skip_mc;

        rv = zebra_get_bank(f,&mc,bmast.links[KMAST_MC-1]);

        if (rv) {
            fprintf(stderr, "error getting MC bank: %s\n", zebra_err);
            goto err;
        }

        if (bmast.links[KMC_MCGN-1] == 0) {
            fprintf(stderr, "MCGN link is zero!\n");
            goto err;
        }

        unpack_mc(mc.data, &bmc);

        rv = zebra_get_bank(f,&bmcgn,mc.links[KMC_MCGN-1]);

        if (rv) {
            fprintf(stderr, "error getting MCGN bank: %s\n", zebra_err);
            goto err;
        }

        while (1) {
            if (bmcgn.links[KMCGN_MCTK-1] == 0) {
                fprintf(stderr, "MCTK link is zero!\n");
                goto err;
            }

            rv = zebra_get_bank(f,&mctk,bmcgn.links[KMCGN_MCTK-1]);

            if (rv) {
                fprintf(stderr, "error getting MCTK bank: %s\n", zebra_err);
                goto err;
            }

            if (mctk.orig == mctk.up - KMCVX_MCTK) {
                /* This is the first MCTK bank. */
                unpack_mctk(mctk.data, &bmctk);
            } else {
                /* For some reason SNOMAN sometimes links to the second MCTK
                 * from the MCGN bank. */
                rv = zebra_get_bank(f,&b,mctk.orig);

                if (b.idh != MCTK_RECORD) {
                    fprintf(stderr, "error following origin link from MCTK bank!\n");
                    goto err;
                }

                if (rv) {
                    fprintf(stderr, "error getting MCTK bank: %s\n", zebra_err);
                    goto err;
                }

                unpack_mctk(b.data, &bmctk);
            }

            if (mctk.up == 0) {
                fprintf(stderr, "MCVX link is zero!\n");
                goto err;
            }

            rv = zebra_get_bank(f,&b,mctk.up);

            if (rv) {
                fprintf(stderr, "error getting MCVX bank: %s\n", zebra_err);
                goto err;
            }

            unpack_mcvx(b.data, &bmcvx);

            if (output) {
                hdf5_mcgn[nmcgn].evn = bmc.evn;
                hdf5_mcgn[nmcgn].id = bmctk.idp;
                hdf5_mcgn[nmcgn].energy = bmctk.ene;
                hdf5_mcgn[nmcgn].x = bmcvx.x;
                hdf5_mcgn[nmcgn].y = bmcvx.y;
                hdf5_mcgn[nmcgn].z = bmcvx.z;
                hdf5_mcgn[nmcgn].dirx = bmctk.drx;
                hdf5_mcgn[nmcgn].diry = bmctk.dry;
                hdf5_mcgn[nmcgn].dirz = bmctk.drz;
                hdf5_mcgn[nmcgn].time = bmcvx.tim;
                nmcgn++;
            }

            if (bmcgn.next) {
                rv = zebra_get_bank(f,&bmcgn,bmcgn.next);

                if (rv) {
                    fprintf(stderr, "error getting MCGN bank: %s\n", zebra_err);
                    goto err;
                }
            } else {
                break;
            }
        }

skip_mc:

        rv = zebra_get_bank(f,&b,bmast.links[KMAST_EV-1]);

        if (rv) {
            fprintf(stderr, "error getting EV bank: %s\n", zebra_err);
            goto err;
        }

        /* Skip to the last event so we can traverse them in reverse order. The
         * reason for this is that for some reason SNOMAN puts the events in
         * reverse order within each logical record. */
        while (b.next) {
            rv = zebra_get_bank(f,&b,b.next);

            if (rv) {
                fprintf(stderr, "error getting EV bank: %s\n", zebra_err);
                goto err;
            }
        }

        while (1) {
            unpack_ev(b.data, &bev);

            if (bev.run != last_run) {
                fprintf(stderr, "loading DQXX file for run %010i\n", bev.run);

                sprintf(dqxx_file, "DQXX_%010i.dat", bev.run);
                if (load_file(db, dqxx_file, 1)) {
                    fprintf(stderr, "failed to load %s: %s\n", dqxx_file, db_err);
                    goto err;
                }

                if (dqxx_init(db, &ev)) {
                    fprintf(stderr, "failed to initialize DQXX bank: %s\n", dqxx_err);
                    goto err;
                }

                last_run = bev.run;
            }

            rv = get_event(f,&ev,&b);

            if (output) {
                hdf5_events[nevents].run = ev.run;
                hdf5_events[nevents].evn = bmc.evn;
                hdf5_events[nevents].gtr = ev.trigger_time;
                hdf5_events[nevents].nhit = ev.nhit;
                hdf5_events[nevents].gtid = ev.gtid;
                hdf5_events[nevents].trg_type = ev.trigger_type;
                hdf5_events[nevents].dc = get_dc_word(&ev, f, &bmast, &b);

                if (get_ftpv(f,&b,&bftpv)) {
                    if (verbose) fprintf(stderr, "%s\n", zdab_err);
                    hdf5_events[nevents].ftp_x = NAN;
                    hdf5_events[nevents].ftp_y = NAN;
                    hdf5_events[nevents].ftp_z = NAN;
                } else {
                    hdf5_events[nevents].ftp_x = bftpv.x;
                    hdf5_events[nevents].ftp_y = bftpv.y;
                    hdf5_events[nevents].ftp_z = bftpv.z;
                }

                if (get_ftxk(f,&b,&bftxk)) {
                    if (verbose) fprintf(stderr, "%s\n", zdab_err);
                    hdf5_events[nevents].ftk_energy = NAN;
                } else {
                    hdf5_events[nevents].ftk_energy = bftxk.energy;
                }

                if (get_rsp(f,&b,&bftxr)) {
                    if (verbose) fprintf(stderr, "%s\n", zdab_err);
                    hdf5_events[nevents].rsp_energy = NAN;
                } else {
                    hdf5_events[nevents].rsp_energy = bftxr.ene;
                }

                nevents++;
            }

            /* Note the origin link for the first EV bank points back to the
             * structural link location in the MAST bank. These links are super
             * confusing! */
            if ((b.orig == f->first_bank - KMAST_EV) || skip_second_event) break;

            rv = zebra_get_bank(f,&b,b.orig);

            if (rv) {
                fprintf(stderr, "error getting EV bank: %s\n", zebra_err);
                goto err;
            }
        }
    }

    if (output) {
        save_output(output, hdf5_events, nevents, hdf5_mcgn, nmcgn, NULL, 0);
    }

    db_free(db);

    zebra_close(f);

    return 0;

err:
    db_free(db);

    zebra_close(f);

    return 1;
}