add online monitoring scripts

include/MonitorDisplay.h

+#pragma once
+
+R__LOAD_LIBRARY(libScandalizer.so)
+#include "scandalizer/PostProcessors.h"
+#include "scandalizer/ScriptHelpers.h"
+#include "scandalizer/SpectrometerMonitor.h"
+#include "monitor/DetectorDisplay.h"
+#include "monitor/EventDisplays.h"
+#include "monitor/ExperimentMonitor.h"
+#include "TStyle.h"
+
+hallc::MonitoringDisplay *CreateMonitorDisplay(const std::string &root_dir, int port, TTree *tree)
+{
+    gStyle->SetHistFillStyle(3002);
+    gStyle->SetPalette(kBird, 0, 0.5);
+    gStyle->SetOptStat(10);
+    gStyle->SetTitleSize(0.04, "");
+    gStyle->SetTitleSize(0.04, "xyz");
+    gStyle->SetLabelSize(0.03, "xyz");
+    // gStyle->SetTitleW(1.0);
+    gStyle->SetTitleStyle(0);
+    gStyle->SetTitleBorderSize(0);
+    gROOT->SetBatch();
+
+    auto dply = new hallc::MonitoringDisplay(port);
+    dply->SetRootFolder(root_dir);
+
+    // cherenkov channels
+    dply->CreateDisplayPlot("cher", "cher",
+        [] (hallc::DisplayPlot& plt) {
+            plt._plot_data._canvas = new TCanvas("cher", "MaPMT Channels", 1280, 720); // 720p HD
+            plt._plot_data._canvas->Divide(4, 4);
+            for (int i = 1; i <= 4; ++i) {
+                for (int j = 1; j <= 4; ++j) {
+                    plt._plot_data._canvas->cd((i - 1)*4 + j);
+                    auto chn = fmt::format("Cer{}{}_5_Ppeak", i, j);
+                    auto h1f = new TH1F(("h" + chn).c_str(), chn.c_str(), 400, 0, 4096);
+                    plt._plot_data._hists1.push_back(h1f);
+                    h1f->Draw("PFC");
+                    gPad->SetLogy();
+                    h1f->GetYaxis()->SetRangeUser(1, 1e5);
+                }
+            }
+            return 0;
+        },
+        [tree] (hallc::DisplayPlot& plt) {
+            // std::cout << "update canvas" << std::endl;
+            for (size_t i = 0; i < plt._plot_data._hists1.size(); ++i) {
+                TH1F *h1f = plt._plot_data._hists1[i];
+                plt._plot_data._canvas->cd(i + 1);
+                std::string name = h1f->GetName();
+                std::string title = h1f->GetTitle();
+                tree->Draw(fmt::format("{} >> {}", title, name).c_str());
+                // std::cout << name << ": " << h1f->Integral(1, 400) << std::endl;
+            }
+            
+            return 0;
+        }
+    );
+
+    return dply;
+}
+

online/monitor.cxx

+#include "MonitorDisplay.h"
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <csignal>
+#include "utils.h"
+#include "ConfigParser.h"
+#include "ConfigObject.h"
+#include "ConfigOption.h"
+#include "TSpectrum.h"
+#include "TTree.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "simpleLib.h"
+#include "Fadc250Decoder.h"
+#include "ETChannel.h"
+
+
+#define FADC_BANK 3
+#define PROGRESS_COUNT 100
+
+
+volatile sig_atomic_t term = 0;
+
+void handle_sig(int signum)
+{
+    std::cout << "Received SIGINT, stopping the monitor..." << std::endl;
+    term = 1;
+}
+
+
+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;
+}
+
+uint32_t parseBlock(const uint32_t *buf)
+{
+    auto header = (ETChannel::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;
+}
+
+// parse an evio event
+bool parseEvent(const uint32_t *buf, bool verbose = false)
+{
+    auto header = (ETChannel::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 ETChannel::CODA_PHY1:
+    case ETChannel::CODA_PHY2:
+        break;
+    case ETChannel::CODA_PRST:
+    case ETChannel::CODA_GO:
+    case ETChannel::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;
+    */
+}
+
+#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 processEvent(const uint32_t *buf, int &count, TTree *tree, bool &init_tree, int &data_mode,
+                  const std::vector<Module> &modules, std::unordered_map<std::string, BranchData> &brdata)
+{
+    if (!parseEvent(buf, false)) {
+        return;
+    }
+
+    // get block level
+    int blvl = 1;
+    simpleGetRocBlockLevel(modules.front().crate, FADC_BANK, &blvl);
+
+    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 (data_mode < 0 && slot_data.GetMode() > 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 ++;
+    }
+}
+
+void monitor(const std::string &addr = "clrlpc", int port = 11111, const std::string &etfile = "/tmp/et_test",
+             const std::string &opath = "processed_data/online.root",
+             const std::string &module_conf = "config/esb_module.conf",
+             int nev = -1)
+{
+    auto modules = read_modules(module_conf);
+    ETChannel et_chan(100, 100);
+    if (et_chan.Connect(addr.c_str(), port, etfile.c_str(), "TCD_MONITOR") != ET_OK) {
+        return;
+    }
+
+    auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT test results");
+    auto tree = new TTree("EvTree", "Cherenkov Test Events");
+    // 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();
+        }
+    }
+
+    auto dply = CreateMonitorDisplay("/online/latest/", 9100, tree);
+    dply->InitAll();
+
+    int count = 0;
+    bool init_tree = false;
+    int data_mode = -1;
+    int status = ETChannel::READ_EMPTY;
+    et_chan.Open();
+    signal (SIGINT, handle_sig);
+
+    while (true) {
+        if (term > 0 || (nev > 0 && nev < count)) { break; }
+        if (count % PROGRESS_COUNT == 0 && count > 0) {
+            std::cout << "Read and processed events - " << count << std::endl;
+            dply->Process();
+            dply->UpdateAll();
+        }
+
+        status = et_chan.ReadEvent();
+
+        switch (status) {
+        case ETChannel::READ_ERROR:
+        case ETChannel::READ_EOF:
+            term = 1;
+            break;
+        case ETChannel::READ_EMPTY:
+            std::cout << "ET station is empty, wait 2 secs..." << std::endl;
+            std::this_thread::sleep_for(std::chrono::seconds(2));
+            break;
+        case ETChannel::READ_OK:
+            processEvent(et_chan.GetEvBuffer(), count, tree, init_tree, data_mode, modules, brdata);
+            break;
+        }
+    }       
+
+    dply->Process();
+    dply->UpdateAll();
+    std::cout << "Read and processed events - " << count << std::endl;
+
+    hfile->Write();
+    hfile->Close();
+}
+