diff --git a/include/MonitorDisplay.h b/include/MonitorDisplay.h new file mode 100644 index 0000000000000000000000000000000000000000..dcee20ad39889689a35fee208dfd3e7b3c9d5120 --- /dev/null +++ b/include/MonitorDisplay.h @@ -0,0 +1,63 @@ +#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; +} + diff --git a/online/monitor.cxx b/online/monitor.cxx new file mode 100644 index 0000000000000000000000000000000000000000..46af8c3c319636c18e962d713dc629c244573510 --- /dev/null +++ b/online/monitor.cxx @@ -0,0 +1,379 @@ +#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(); +} +