esb_analyze.cpp

#include <iostream>
#include <iomanip>
#include <fstream>
#include "utils.h"
#include "ConfigParser.h"
#include "ConfigObject.h"
#include "ConfigOption.h"
#include "THaCodaFile.h"
#include "THaEvData.h"
#include "CodaDecoder.h"
#include "evio.h"
#include "TSpectrum.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "simpleLib.h"
#include "Fadc250Decoder.h"


#define FADC_BANK 3
#define PROGRESS_COUNT 1000


void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules);
bool parseEvent(const uint32_t *buf, bool verbose);
uint32_t parseBlock(const uint32_t *buf);

/* 32 bit event header structure
 * --------------------------------
 * |             length           |
 * --------------------------------
 * |    etype     | dtype |  num  |
 * --------------------------------
 */
struct CodaEvHeader
{
    uint32_t length;
    uint8_t num, dtype;
    uint16_t etype;
};

enum EvType {
    CODA_PRST = 0xffd1,
    CODA_GO = 0xffd2,
    CODA_END = 0xffd4,
    CODA_PHY1 = 0xff50,
    CODA_PHY2 = 0xff70,
};

uint32_t swap_endian32(uint32_t num)
{
    uint32_t b0 = (num & 0x000000ff) << 24u;
    uint32_t b1 = (num & 0x0000ff00) << 8u;
    uint32_t b2 = (num & 0x00ff0000) >> 8u;
    uint32_t b3 = (num & 0xff000000) >> 24u;
    return b0 | b1 | b2 | b3;
}

// parse an evio event
bool parseEvent(const uint32_t *buf, bool verbose = false)
{
    auto header = (CodaEvHeader*)buf;
    const uint32_t buf_size = header->length + 1;
    if (verbose) {
        std::cout << "Event header: " << std::dec << buf_size << "\n"
                  << (int) header->etype << ", " << (int) header->dtype << ", " << (int) header->num
                  << std::endl;
        std::cout << std::hex;
        for (uint32_t i = 0; i < 2; ++i) {
            std::cout << "0x" << std::setw(8) << std::setfill('0') << buf[i] << "\n";
        }
    }

    switch(header->etype) {
    case CODA_PHY1:
    case CODA_PHY2:
        break;
    case CODA_PRST:
    case CODA_GO:
    case CODA_END:
    default:
        return false;
    }

    simpleScan((volatile uint32_t*)buf, buf_size);
    return true;

    /*
    // parse ROC data
    uint32_t index = 2;
    while(index < buf_size) {
        // skip header size and data size 2 + (length - 1)
        index += parseBlock(&buf[index]);
    }

    // return parsed event type
    return buf_size;
    */
}


uint32_t parseBlock(const uint32_t *buf)
{
    auto header = (CodaEvHeader*) buf;
    const uint32_t buf_size = header->length + 1;
    std::cout << "ROC header: " << std::dec << buf_size << "\n"
              << (int) header->etype << ", " << (int) header->dtype << ", " << (int) header->num
              << std::endl;
    std::cout << std::hex;
    for (uint32_t i = 0; i < buf_size; ++i) {
        std::cout << "0x" << std::setw(8) << std::setfill('0') << buf[i] << "\n";
    }
    return buf_size;
}


int main(int argc, char* argv[])
{
    // setup input arguments
    ConfigOption conf_opt;
    conf_opt.AddOpts(ConfigOption::help_message, 'h', "help");
    conf_opt.AddOpt(ConfigOption::arg_require, 'n');
    conf_opt.AddLongOpt(ConfigOption::arg_require, "config-module", 'm');
    conf_opt.AddLongOpt(ConfigOption::arg_require, "waveform-output", 'w');

    conf_opt.SetDesc("usage: %0 <data_file> <out_file>");
    conf_opt.SetDesc('n', "number of events to process (< 0 means all).");
    conf_opt.SetDesc('m', "configuration file for modules to be read-in, default \"config/mapmt_module.conf\"");

    if (!conf_opt.ParseArgs(argc, argv) || conf_opt.NbofArgs() != 2) {
        std::cout << conf_opt.GetInstruction() << std::endl;
        return -1;
    }

    int nev = -1;
    std::string mconf = "config/esb_module.conf";

    for (auto &opt : conf_opt.GetOptions()) {
        switch (opt.mark) {
        case 'n':
            nev = opt.var.Int();
            break;
        case 'm':
            mconf = opt.var.String();
            break;
        default :
            std::cout << conf_opt.GetInstruction() << std::endl;
            return -1;
        }
    }

    auto modules = read_modules(mconf);

    write_raw_data(conf_opt.GetArgument(0), conf_opt.GetArgument(1), nev, modules);
}


#define MAX_NPEAKS 20
#define MAX_RAW 300
struct BranchData
{
    int npul, nraw;
    float integral[MAX_NPEAKS], peak[MAX_NPEAKS], time[MAX_NPEAKS];
    int raw[MAX_RAW];
    float ped_mean, ped_err;
};

#define BUF_SIZE 1000
static double buffer[BUF_SIZE], wfbuf[BUF_SIZE], bkbuf[BUF_SIZE];
void refine_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
{
    int count = 0;
    float mean = 0.;
    for (auto &val : raw) {
        if (std::abs(val - ped_mean) < 2.0 * ped_err) {
            buffer[count] = val;
            mean += val;
            count ++;
        }
    }
    if (count == 0) {
        return;
    }

    mean /= count;
    float change = std::abs(ped_mean - mean);

    // no outliers
    if (change < 0.05 * ped_err) {
        return;
    }

    // recursively refine
    float err = 0.;
    for (int i = 0; i < count; ++i) {
        err += (buffer[i] - mean)*(buffer[i] - mean);
    }
    ped_mean = mean;
    ped_err = std::sqrt(err/count);
    refine_pedestal(raw, ped_mean, ped_err);
}

// a function for a simple constant baseline fit
void find_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
{
    ped_mean = 0;
    ped_err = 0;

    for (auto &val : raw) {
        ped_mean += val;
    }
    ped_mean /= raw.size();

    for (auto &val : raw) {
        ped_err += (val - ped_mean)*(val - ped_mean);
    }
    ped_err = std::sqrt(ped_err/raw.size());

    refine_pedestal(raw, ped_mean, ped_err);
}

// analyze the waveform data and fill the result in a branch data
void waveform_analysis(const std::vector<uint32_t> &raw, BranchData &res)
{
    // no need to analyze
    if (raw.empty()) {
        return;
    }

    // fill in the raw data
    res.nraw = raw.size();
    for (size_t i = 0; i < raw.size(); ++i) {
        res.raw[i] = raw[i];
    }

    // get pedestals
    find_pedestal(raw, res.ped_mean, res.ped_err);

    // fill in spectrum buffer
    for (size_t i = 0; i < raw.size(); ++i) {
        wfbuf[i] = raw[i] - res.ped_mean;
    }

    // find peaks
    TSpectrum s;
    s.SetResolution(0.5);
    int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 3.0, 20, false, 5, true, 3);

    // fill branch data
    double *pos = s.GetPositionX();
    res.npul = 0;
    for (int i = 0; i < npeaks; ++i) {
        int j = pos[i];
        if (wfbuf[j] < 5.*res.ped_err) { continue; }
        res.time[res.npul] = j;
        res.peak[res.npul] = wfbuf[j];
        /*
        res.integral[res.npul] = wfbuf[j];
        j = pos[i] - 1;
        while ( (j > 0) && (wfbuf[j] - wfbuf[j - 1])*wfbuf[j] > 0. ) {
            res.integral[res.npul] += wfbuf[j--];
        }
        j = pos[i] + 1;
        while ( (j < res.nraw - 1) && (wfbuf[j] - wfbuf[j + 1])*wfbuf[j] > 0. ) {
            res.integral[res.npul] += wfbuf[j++];
        }
        */
        res.npul ++;
    }
}


// fill decoded FADC250 data into the container (branch data)
void fill_branch(const Fadc250Data &slot_data, const Module &mod, std::unordered_map<std::string, BranchData> &brdata)
{
    for (auto &event : slot_data.events) {
        for (auto &ch : mod.channels) {
            auto &channel = event.channels[ch.id];
            auto &bd = brdata[ch.name];
            bd.npul = channel.integral.size();
            for (uint32_t i = 0; i < channel.integral.size(); ++i) {
                bd.integral[i] = channel.integral[i];
                bd.time[i] = channel.time[i];
            }
            waveform_analysis(channel.raw, bd);
        }
    }
}

// fill branch data into the root tree
void fill_tree(TTree *tree, std::unordered_map<std::string, BranchData> &brdata, bool &init, int mode)
{
    if ((mode != 1) && (mode != 3)) {
        std::cout << "Warning: unsupported mode " << mode << ", data won't be recorded." << std::endl;
        return;
    }

    // initialize
    if (!init) {
        for (auto &it : brdata) {
            auto n = it.first;
            tree->Branch((n + "_Npulse").c_str(), &brdata[n].npul, (n + "_N/I").c_str());
            tree->Branch((n + "_Pint").c_str(), &brdata[n].integral[0], (n + "_Pint[" + n + "_N]/F").c_str());
            tree->Branch((n + "_Ptime").c_str(), &brdata[n].time[0], (n + "_Ptime[" + n + "_N]/F").c_str());
            // raw waveform provides more information
            if (mode == 1) {
                tree->Branch((n + "_Ppeak").c_str(), &brdata[n].peak[0], (n + "_Ppeak[" + n + "_N]/F").c_str());
                tree->Branch((n + "_Nraw").c_str(), &brdata[n].nraw, (n + "_Nraw/I").c_str());
                tree->Branch((n + "_raw").c_str(), &brdata[n].raw[0], (n + "_raw[" + n + "_Nraw]/I").c_str());
                tree->Branch((n + "_ped_mean").c_str(), &brdata[n].ped_mean, (n + "_ped_mean/F").c_str());
                tree->Branch((n + "_ped_err").c_str(), &brdata[n].ped_err, (n + "_ped_err/F").c_str());
            }
        }
        init = true;
        std::cout << "Initialized root file for mode " << mode << " data." << std::endl;
    }

    tree->Fill();
}


void write_raw_data(const std::string &dpath, const std::string &opath, int nev, const std::vector<Module> &modules)
{
    Decoder::THaCodaFile datafile(dpath.c_str());
    THaEvData *evdata = new Decoder::CodaDecoder();

    auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
    auto tree = new TTree("EvTree", "Cherenkov Test Events");
    int evtype;
    tree->Branch("Event Type", &evtype, "EventType/I");
    // set up branches and a container to be combined with the branches
    std::unordered_map<std::string, BranchData> brdata;
    for (auto &mod : modules) {
        for (auto &ch : mod.channels) {
            auto n = ch.name;
            if (brdata.find(n) != brdata.end()) {
                std::cout << "WARNING: Duplicated channel names found! -- " << n << std::endl;
            }
            brdata[n] = BranchData();
        }
    }

    int count = 0;
    bool init_tree = false;
    while ((datafile.codaRead() == S_SUCCESS) && (nev-- != 0)) {
        // read-in data
        evdata->LoadEvent(datafile.getEvBuffer());
        evtype = evdata->GetEvType();

        if((count % PROGRESS_COUNT) == 0) {
            std::cout << "Processed events - " << count << "\r" << std::flush;
        }

        if (!parseEvent(evdata->GetRawDataBuffer())) {
            continue;
        }

        // get block level
        int blvl = 1;
        simpleGetRocBlockLevel(modules.front().crate, FADC_BANK, &blvl);

        int data_mode = -1;
        for (int ii = 0; ii < blvl; ++ii) {
            // clear data buffer
            for (auto &br : brdata) {
                br.second.npul = 0;
                br.second.nraw = 0;
            }

            // parse module data
            for (auto &mod : modules) {
                uint32_t header = 0;
                auto status = simpleGetSlotBlockHeader(mod.crate, FADC_BANK, mod.slot, &header);
                if (status <= 0) {
                    std::cout << "Error getting header for crate = "
                              << mod.crate << ", slot = " << mod.slot << std::endl;
                    continue;
                }

                uint32_t *dbuf;
                auto len = simpleGetSlotEventData(mod.crate, FADC_BANK, mod.slot, ii, &dbuf);
                if (len <= 0) {
                    std::cout << "No data for crate = " << mod.crate << ", slot = " << mod.slot
                              << ", block_level = " << ii << std::endl;
                    continue;
                }
                auto slot_data = fadc250_decoder(header, dbuf, len);

                // check data mode
                if (slot_data.GetMode() > 0 && data_mode < 0) {
                    data_mode = slot_data.GetMode();
                }
                // fill branch data
                fill_branch(slot_data, mod, brdata);
            }

            // fill event
            fill_tree(tree, brdata, init_tree, data_mode);
            count ++;
        }

    }
    std::cout << "Processed events - " << count << std::endl;

    hfile->Write();
    hfile->Close();
}