readWD.cc

Go to the documentation of this file.

 
#include "readWD.hh"
 
using namespace std;
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ FUNCTIONS                                                               │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
ostream &operator<<(ostream &o, const TAG &tag) // cout << TAG
{
    if (tag.tag[0] == 'B')
    {
        o << tag.tag[0] << tag.tag[1] << *(short *)(tag.tag + 2);
    }
    else if (tag.tag[0] == 'T' and tag.tag[1] == '#')
    {
        o << tag.tag[0] << tag.tag[1] << *(short *)(tag.tag + 2);
    }
    else
    {
        o << tag.tag;
    }
    return o;
}
 
ostream &operator<<(ostream &o, const EventHeader &eh) // cout << EventHeader
{
    o << eh.tag << endl;
    o << "Date: " << eh.year << "/" << eh.month << "/" << eh.day << endl;
    o << "Hour: " << eh.hour << ":" << eh.min << ":" << eh.sec << "." << eh.ms << endl;
    o << "Range: " << eh.rangeCenter << " Serial number: " << eh.serialNumber;
    return o;
}
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ CLASSES : DAQEvent                                                      │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
DAQEvent::DAQEvent()
{
    is_getch_ = false;
    is_init_ = false;
    is_iw_ = false;
    is_ped_ = false;
    is_peak_ = false;
    ped_interval_ = {0, 100};
    iw_ = {0, SAMPLES_PER_WAVEFORM - 1};
    peak_threshold_ = 1.;
    routine_ = {false, false, false};
    ch_old_ = {-1, -1};
}
 
DAQEvent &DAQEvent::GetChannel(const int &board, const int &channel)
{
    if (board < 0 or channel < 0)
    {
        cerr << "!! Error: board and channel ID number(s) must be positive integers" << endl;
        exit(0);
    }
 
    else if (ch_.first == board and ch_.second == channel)
    {
        is_getch_ = true;
        return *this;
    }
 
    else if (volts_.find(board) != volts_.end())
    {
        if (volts_[board].find(channel) != volts_[board].end())
        {
            is_getch_ = true;
            ch_ = {board, channel};
            routine_ = {false, false, false};
            return *this;
        }
        cerr << "!! Error: invalid channel, max channel ID number for this board is " << volts_[board].size() - 1 << endl;
        exit(0);
    }
 
    else if (is_init_ == false && config_.is_makeconfig_ == true)
    {
        is_getch_ = true;
        ch_ = {board, channel};
        return *this;
    }
    cerr << "!! Error: invalid board, max board ID number is " << volts_.size() - 1 << endl;
    exit(0);
}
 
void DAQEvent::SetPedInterval(int a, int b)
{
    if (is_getch_)
    {
        config_.SetPedInterval(pair<int, int>{a, b}, ch_.first, ch_.second);
    }
    else
    {
        config_.SetPedInterval(pair<int, int>{a, b});
    }
    is_getch_ = false;
    return;
}
 
void DAQEvent::SetPeakThr(float thr)
{
    if (is_getch_)
    {
        config_.SetPeakThr(thr, ch_.first, ch_.second);
    }
    else
    {
        config_.SetPeakThr(thr);
    }
    is_getch_ = false;
    return;
}
 
void DAQEvent::SetIntWindow(int a, int b)
{
    if (is_getch_)
    {
        config_.SetIntWindow({a, b}, ch_.first, ch_.second);
    }
    else
    {
        config_.SetIntWindow({a, b});
    }
    is_getch_ = false;
    return;
}
 
void DAQEvent::SetIntWindow(float a, float b)
{
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    auto &times = times_[ch_.first][ch_.second];
    if (a < times[0] or a > b or b > times[SAMPLES_PER_WAVEFORM - 1])
    {
        cerr << "!! Error: invalid times passed as integration window" << endl;
        exit(0);
    }
 
    config_.intWindow_[ch_.first][ch_.second].first = distance(times.begin(), lower_bound(times.begin(), times.end(), a));
 
    auto &iw_first = config_.intWindow_[ch_.first][ch_.second].first;
    config_.intWindow_[ch_.first][ch_.second].second = distance(times.begin(), lower_bound(times.begin() + iw_first, times.end(), b));
    return;
}
 
bool DAQEvent::IsSaturated()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    auto &volts = volts_[ch_.first][ch_.second];
    return any_of(volts.begin() + 2, volts.end() - 2, [](float val)
                  { return (val < -0.499) || (val > +0.499); });
}
 
float DAQEvent::GetCharge()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    (*this).EvalPedestal();
    (*this).FindPeaks();
    (*this).EvalIntegrationBounds();
 
    const auto &volts = volts_[ch_.first][ch_.second];
    const auto &times = times_[ch_.first][ch_.second];
 
    is_getch_ = false;
    iw_ = config_.intWindow_[ch_.first][ch_.second];
    float charge = 0;
 
    for (int i = iw_.first; i < iw_.second; ++i)
    {
        charge += (volts[i + 1] + volts[i] - 2 * ped_.first) / (2 * (times[i + 1] - times[i]));
    }
 
    return abs(charge);
}
 
float DAQEvent::GetAmplitude()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    (*this).EvalPedestal();
    (*this).FindPeaks();
 
    is_getch_ = false;
 
    return peak_.first - ped_.first;
}
 
float DAQEvent::GetTime(float thr)
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    (*this).EvalPedestal();
 
    auto &volts = volts_[ch_.first][ch_.second];
    auto &times = times_[ch_.first][ch_.second];
    int i = 10;
 
    if (thr < ped_.first)
    {
        while (volts[i] > thr && i < SAMPLES_PER_WAVEFORM - 10)
        {
            ++i;
        }
    }
    else
    {
        while (volts[i] < thr && i < SAMPLES_PER_WAVEFORM - 10)
        {
            ++i;
        }
    }
 
    if (i == SAMPLES_PER_WAVEFORM)
    {
        return 0;
    }
 
    auto time = times[i] + (thr - volts[i]) * (times[i + 1] - times[i]) / (volts[i + 1] - volts[i]);
 
    is_getch_ = false;
 
    return time;
}
 
float DAQEvent::GetTimeCF(float CF)
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (CF <= 0 or CF > 1)
    {
        cerr << "!! Error: CF value must be in range (0, 1)" << endl;
        exit(0);
    }
 
    (*this).EvalPedestal();
    (*this).FindPeaks();
 
    float thr = ped_.first + (peak_.first - ped_.first) * CF;
 
    return (*this).GetTime(thr);
}
 
float DAQEvent::GetRiseTime()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    return (*this).GetTimeCF(0.9) - (*this).GetTimeCF(0.1);
}
 
const pair<float, float> &DAQEvent::GetPedestal()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel with DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    (*this).EvalPedestal();
    is_getch_ = false;
    return ped_;
}
 
const vector<float> &DAQEvent::GetVolts()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel with DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    is_getch_ = false;
    return volts_[ch_.first][ch_.second];
}
 
const vector<float> &DAQEvent::GetTimes()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel with DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    is_getch_ = false;
    return times_[ch_.first][ch_.second];
}
 
const vector<int> &DAQEvent::GetPeakIndices()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
    (*this).EvalPedestal();
    (*this).FindPeaks();
 
    is_getch_ = false;
 
    return indexMin_;
}
 
const pair<int, int> &DAQEvent::GetIntegrationBounds()
{
    if (!is_init_)
    {
        cerr << "!! Error: no event read yet" << endl;
        exit(0);
    }
 
    if (!config_.user_iw_[ch_.first][ch_.second])
    {
        (*this).EvalPedestal();
        (*this).FindPeaks();
        (*this).EvalIntegrationBounds();
    }
 
    is_getch_ = false;
 
    return config_.intWindow_[ch_.first][ch_.second];
}
 
DAQEvent &DAQEvent::TimeCalibration(const unsigned short &tCell, const std::vector<float> &times, int i, int j)
{
    times_[i][j] = times;
    vector<float> &times_ij = times_[i][j];
    rotate(times_ij.begin(), times_ij.begin() + tCell, times_ij.end());
    partial_sum(times_ij.begin(), times_ij.end(), times_ij.begin());
 
    return *this;
}
 
DAQEvent &DAQEvent::EvalPedestal()
{
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    auto same_ch = (ch_old_.first == ch_.first) && (ch_old_.second == ch_.second);
    auto same_evt = (evtserial_old_ == eh_.serialNumber);
 
    if (routine_[0] && same_ch && same_evt)
    {
        return *this;
    }
    else
    {
        routine_[0] = true;
        ch_old_ = ch_;
        evtserial_old_ = eh_.serialNumber;
    }
 
    ped_interval_ = config_.pedInterval_[ch_.first][ch_.second];
    int ped_interval_dist = ped_interval_.second - ped_interval_.first;
    const vector<float> &volts = volts_[ch_.first][ch_.second];
    ped_ = {0., 0.};
    ped_.first = accumulate(volts.begin() + ped_interval_.first, volts.begin() + ped_interval_.second, 0.) / ped_interval_dist;
    for (int i = ped_interval_.first; i < ped_interval_.second; ++i)
    {
        ped_.second += pow(volts[i] - ped_.first, 2);
    }
    ped_.second = sqrt(ped_.second / ped_interval_dist);
 
    return *this;
}
 
DAQEvent &DAQEvent::EvalIntegrationBounds()
{
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    if (config_.user_iw_[ch_.first][ch_.second])
    {
        return *this;
    }
 
    auto same_ch = (ch_old_.first == ch_.first) && (ch_old_.second == ch_.second);
    auto same_evt = (evtserial_old_ == eh_.serialNumber);
 
    if (routine_[2] && same_ch && same_evt)
    {
        return *this;
    }
    else
    {
        routine_[2] = true;
        ch_old_ = ch_;
        evtserial_old_ = eh_.serialNumber;
    }
 
    iw_ = {indexMin_[0], indexMin_[0]};
 
    const vector<float> &volts = volts_[ch_.first][ch_.second];
    auto lower_bound = ped_.first - 5 * ped_.second;
 
    if (peak_.first < lower_bound)
    {
        while (volts[iw_.first] < lower_bound and iw_.first > 10)
        {
            --iw_.first;
        }
 
        while (volts[iw_.second] < lower_bound and iw_.second < SAMPLES_PER_WAVEFORM - 10)
        {
            ++iw_.second;
        }
    }
 
    config_.intWindow_[ch_.first][ch_.second] = {iw_.first, iw_.second};
    return *this;
}
 
DAQEvent &DAQEvent::FindPeaks()
{
    if (!is_getch_)
    {
        cerr << "!! Error: select a channel using DAQEvent::GetChannel()" << endl;
        exit(0);
    }
 
    auto same_ch = (ch_old_.first == ch_.first) && (ch_old_.second == ch_.second);
    auto same_evt = (evtserial_old_ == eh_.serialNumber);
 
    if (routine_[1] && same_ch && same_evt)
    {
        return *this;
    }
    else
    {
        routine_[1] = true;
        ch_old_ = ch_;
        evtserial_old_ = eh_.serialNumber;
    }
 
    long index_min;
    auto &volts = volts_[ch_.first][ch_.second];
    auto &times = times_[ch_.first][ch_.second];
    indexMin_ = {};
 
    iw_ = config_.intWindow_[ch_.first][ch_.second];
    peak_threshold_ = config_.peakThr_[ch_.first][ch_.second];
 
    if (config_.user_iw_[ch_.first][ch_.second]) // Integration window set by the user
    {
        index_min = distance(volts.begin(), min_element(volts.begin() + iw_.first, volts.begin() + iw_.second));
        indexMin_.push_back(index_min);
    }
    else // No user integration window set
    {
        index_min = distance(volts.begin() + 10, min_element(volts.begin() + 10, volts.end() - 10)) + 10;
        bool signal, min_left, min_right, at_least;
        for (int i = 10; i < SAMPLES_PER_WAVEFORM - 10; ++i)
        {
            signal = abs(volts[i] - ped_.first) > 5 * ped_.second;
            min_left = abs(volts[i]) > abs(volts[i - 1]) + ped_.second;
            min_right = abs(volts[i]) > abs(volts[i + 1]) + ped_.second;
            if (config_.peakThr_[ch_.first][ch_.second] == 0.5) // Default threshold level
            {
                at_least = abs(volts[i] - ped_.first) > abs(volts[index_min] - ped_.first) * 0.5;
            }
            else // Custom threshold level
            {
                at_least = volts[i] < peak_threshold_;
            }
 
            if (signal and min_left and min_right and at_least)
            {
                indexMin_.push_back(i);
            }
        }
 
        if (find(indexMin_.begin(), indexMin_.end(), index_min) == indexMin_.end()) // Sorted insertion of min element's index
        {
            auto pos = find_if(indexMin_.begin(), indexMin_.end(), [index_min](auto i)
                               { return i > index_min; });
            indexMin_.insert(pos, index_min);
        }
    }
 
    if (indexMin_.size() == 0) // Assure that at least global minimum is inserted in indexMin_
    {
        index_min = distance(volts.begin() + 10, min_element(volts.begin() + 10, volts.end() - 10)) + 10;
        indexMin_.push_back(index_min);
    }
 
    peak_ = {volts[indexMin_[0]], times[indexMin_[0]]}; // First local minimum taken as peak
 
    return *this;
}
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ CLASSES : DRSEvent                                                      │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
DRSEvent::DRSEvent()
{
    evtserial_old_ = 1;
}
const string DRSEvent::type_ = "DRS";
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ CLASSES : WDBEvent                                                      │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
WDBEvent::WDBEvent()
{
    evtserial_old_ = 0;
}
const string WDBEvent::type_ = "WDB";
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ CLASSES : DAQConfig                                                     │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
DAQConfig::DAQConfig()
{
    is_makeconfig_ = false;
}
 
void DAQConfig::MakeConfig(DAQFile &file)
{
    if (file.initialization_ == false)
    {
        cerr << "!! Error : File was not initialised, use DAQFile::Open()" << endl;
        exit(0);
    }
 
    is_makeconfig_ = true;
    for (auto &[bKey, bVal] : file.times_)
    {
        for (auto &[cKey, cVal] : bVal)
        {
            intWindow_[bKey][cKey] = {0, SAMPLES_PER_WAVEFORM - 1};
            pedInterval_[bKey][cKey] = {0, 100};
            peakThr_[bKey][cKey] = +0.5;
            user_iw_[bKey][cKey] = false;
        }
    }
}
 
void DAQConfig::ShowConfig()
{
    cout << "----- CONFIGURATION SETTINGS -----" << endl;
    for (auto &[bKey, bVal] : peakThr_)
    {
        for (auto &[cKey, cVal] : bVal)
        {
            cout << " - Board/Channel ID : " << bKey << "/" << cKey << endl
                 << "       - Integration window : (" << intWindow_[bKey][cKey].first << ", " << intWindow_[bKey][cKey].second << ")" << endl
                 << "       - Pedestal interval : (" << pedInterval_[bKey][cKey].first << ", " << pedInterval_[bKey][cKey].second << ")" << endl
                 << "       - Peak threshold : " << cVal << " V" << endl;
        }
    }
}
 
void DAQConfig::SetIntWindow(pair<int, int> intWindow, int b, int c)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (intWindow.first < 0 || intWindow.first > intWindow.second || intWindow.second > SAMPLES_PER_WAVEFORM - 1)
    {
        cerr << "!! Error : Integration window has invalid value" << endl
             << " Values must be in interval (0, " << SAMPLES_PER_WAVEFORM << "), passed values are ( " << intWindow.first << ", " << intWindow.second << ")" << endl;
        exit(0);
    }
 
    if (intWindow_.find(b) != intWindow_.end())
    {
        if (intWindow_[b].find(c) != intWindow_[b].end())
        {
            intWindow_[b][c] = intWindow;
            user_iw_[b][c] = true;
            return;
        }
    }
    cerr << "!! Error : Couldn't find board-channel of ID (" << b << ", " << c << ")" << endl;
    ;
    exit(0);
}
 
void DAQConfig::SetIntWindow(pair<int, int> intWindow)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (intWindow.first < 0 || intWindow.first > intWindow.second || intWindow.second > SAMPLES_PER_WAVEFORM - 1)
    {
        cerr << "!! Error : Integration window has invalid value" << endl
             << " Values must be in interval (0, " << SAMPLES_PER_WAVEFORM << "), passed values are ( " << intWindow.first << ", " << intWindow.second << ")" << endl;
        exit(0);
    }
 
    for (auto &[bKey, bVal] : intWindow_)
    {
        for (auto &[cKey, cVal] : bVal)
        {
            cVal = intWindow;
            user_iw_[bKey][cKey] = true;
        }
    }
}
 
void DAQConfig::SetPedInterval(pair<int, int> pedInterval, int b, int c)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (pedInterval.first < 0 || pedInterval.first > pedInterval.second || pedInterval.second > SAMPLES_PER_WAVEFORM - 1)
    {
        cerr << "!! Error : Pedestal interval has invalid value" << endl
             << " Values must be in interval (0, " << SAMPLES_PER_WAVEFORM << "), passed values are ( " << pedInterval.first << ", " << pedInterval.second << ")" << endl;
        exit(0);
    }
 
    if (pedInterval_.find(b) != pedInterval_.end())
    {
        if (pedInterval_[b].find(c) != pedInterval_[b].end())
        {
            pedInterval_[b][c] = pedInterval;
            return;
        }
    }
    cerr << "!! Error : Couldn't find board-channel of ID (" << b << ", " << c << ")" << endl;
    exit(0);
}
 
void DAQConfig::SetPedInterval(pair<int, int> pedInterval)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (pedInterval.first < 0 || pedInterval.first > pedInterval.second || pedInterval.second > SAMPLES_PER_WAVEFORM - 1)
    {
        cerr << "!! Error : Pedestal interval has invalid value" << endl
             << " Values must be in interval (0, " << SAMPLES_PER_WAVEFORM << "), passed values are ( " << pedInterval.first << ", " << pedInterval.second << ")";
        exit(0);
    }
 
    for (auto &[bKey, bVal] : pedInterval_)
    {
        for (auto &[cKey, cVal] : bVal)
        {
            cVal = pedInterval;
        }
    }
}
 
void DAQConfig::SetPeakThr(float thr, int b, int c)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (thr < -0.5 || thr > +0.5)
    {
        cerr << "!! Error : Peak threshold has invalid value" << endl
             << " Values must be in interval (-0.5, +0.5), passed values is " << thr << endl;
        exit(0);
    }
 
    if (peakThr_.find(b) != peakThr_.end())
    {
        if (peakThr_[b].find(c) != peakThr_[b].end())
        {
            peakThr_[b][c] = thr;
            return;
        }
    }
    cerr << "!! Error : Couldn't find board-channel of ID (" << b << ", " << c << ")" << endl;
    exit(0);
}
 
void DAQConfig::SetPeakThr(float thr)
{
    if (is_makeconfig_ == false)
    {
        cerr << "!! Error : Configuration class not initialised, use DAQEvent::MakeConfig()" << endl;
        exit(0);
    }
 
    if (thr < -0.5 || thr > +0.5)
    {
        cerr << "!! Error : Peak threshold has invalid value" << endl
             << " Values must be in interval (-0.5, +0.5), passed values is " << thr << endl;
        exit(0);
    }
 
    for (auto &[bKey, bVal] : peakThr_)
    {
        for (auto &[cKey, cVal] : bVal)
        {
            cVal = thr;
        }
    }
}
 
/*
  ┌─────────────────────────────────────────────────────────────────────────┐
  │ CLASSES : DAQFile                                                       │
  └─────────────────────────────────────────────────────────────────────────┘
 */
 
DAQFile::DAQFile()
{
    std::cout << "Created DAQFile, open a file using DAQFile::Open()" << endl;
    is_lab_ = 0;
    initialization_ = 0;
}
 
DAQFile::DAQFile(const string &fname)
{
    filename_ = fname;
    in_.open(fname, std::ios::in | std::ios::binary);
    std::cout << "Created DAQFile, opened file " << fname << std::endl;
    initialization_ = 0;
    is_lab_ = 0;
    (*this).Initialise();
}
 
DAQFile &DAQFile::Initialise()
{
    DAQFile &file = *this;
 
    if (!in_.is_open())
    {
        cerr << "!! Error: file not open --> use DAQFile(filename)" << endl;
        return file;
    }
 
    if (initialization_ == 1)
    {
        return file;
    }
 
    TAG bTag, cTag;
    vector<float> times(SAMPLES_PER_WAVEFORM);
 
    cout << "Initializing file " << filename_ << endl;
 
    file >> bTag; // DRSx (TIME for Lab's DRS boards)
    if (bTag.tag[0] == 'D' && bTag.tag[1] == 'R' && bTag.tag[2] == 'S')
    {
        cout << bTag;
        if (bTag.tag[3] == '8')
        {
            cout << " --> WaveDREAM Board" << endl;
            type_ = "WDB";
        }
        else
        {
            cout << " --> DRS Evaluation Board" << endl;
            type_ = "DRS";
        }
        file >> bTag; // TIME
    }
    else if (strcmp(bTag.tag, "TIME") == 0)
    {
        cout << "LAB-DRS" << endl;
        is_lab_ = 1;
        type_ = "DRS";
    }
    else
    {
        cerr << "!! Error: invalid file header --> expected \"DRS\", found " << bTag << endl;
        cerr << "Initialisation failed" << endl;
        return file;
    }
    if (strcmp(bTag.tag, "TIME") != 0)
    {
        cerr << "!! Error: invalid time header --> expected \"TIME\", found " << cTag << endl;
        cerr << "Initialisation failed" << endl;
        return file;
    }
 
    o_ = 'B';
    int i = 0, j = 0;
    while (file >> bTag)
    {
        cout << bTag << ":" << endl;
        j = 0;
        while (file >> cTag)
        {
            cout << " --> " << cTag << endl;
            file.Read(times);
            times_[i][j] = times;
            ++j;
        }
        ++i;
        file.ResetTag();
    }
 
    initialization_ = true;
    file.ResetTag();
    first_evt_pos_ = in_.tellg();
 
    return file;
}
 
DAQFile &DAQFile::Close()
{
    if (in_.is_open())
    {
        cout << "Closing file " << filename_ << "..." << endl;
        in_.close();
    }
    else
    {
        cout << "File is already closed" << endl;
    }
    return *this;
}
 
DAQFile &DAQFile::Open(const string &fname)
{
    if (!in_.is_open())
    {
        filename_ = fname;
        in_.open(fname, std::ios::in | std::ios::binary);
        std::cout << std::endl
                  << "Created DAQFile, opened file " << fname << std::endl;
        (*this).Initialise();
        return *this;
    }
    else
    {
        std::cerr << "!! Error: File is already opened --> " << filename_ << std::endl;
        return *this;
    }
}
 
DAQFile &DAQFile::Reset()
{
    if (initialization_ == 0)
    {
        cerr << "Warning: file is not initialised. Nothing to reset..." << endl;
        return *this;
    }
 
    in_.seekg(first_evt_pos_);
    return *this;
}
 
DAQFile &DAQFile::GetEvent(int evt_id)
{
    DAQFile &file = *this;
    int next_evt_pos;
    int evt_id_pos;
    int evt_size;
    int file_size;
    TAG tag;
    EventHeader eh;
    file.Initialise();
 
    // Evaluate event size
    file.Read(tag); // Reading EHDR
    file.Read(tag); // Reading new line to reset tag.tag
    while (strcmp(tag.tag, "EHDR") != 0)
    {
        file.Read(tag);
    }
    file.ResetTag();
    next_evt_pos = in_.tellg();
    evt_size = next_evt_pos - first_evt_pos_;
    cout << "Event size: " << evt_size << endl;
 
    // Evaluate file size
    in_.seekg(0, in_.end);
    file_size = in_.tellg();
 
    // Reset position
    in_.seekg(first_evt_pos_);
 
    // Check to stay into boundaries of file
    if (file_size - (first_evt_pos_ + evt_id * evt_size) < evt_size)
    {
        cerr << "!! Error : Invalid position reached, out of bounds of file" << endl
             << "Reset position to first event header..." << endl;
    }
    else
    {
        // Moving to requested event
        cout << "Moving to event: " << evt_id << endl;
        in_.seekg(first_evt_pos_ + evt_size * evt_id);
        evt_id_pos = in_.tellg();
        file.Read(eh);
        cout << eh << endl;
        in_.seekg(evt_id_pos);
    }
 
    return file;
}
 
bool DAQFile::operator>>(TAG &t) // DAQFile >> TAG
{
    if (!in_.good())
    {
        return 0;
    }
    this->Read(t);
    return *this;
}
 
bool DAQFile::operator>>(EventHeader &eh) // DAQFile >> EventHeader
{
    if (!in_.good())
    {
        return 0;
    }
    this->Read(eh);
    return *this;
}
 
bool DAQFile::operator>>(DRSEvent &event) // DAQFile >> DRSEvent
{
    if (!in_.good())
    {
        return 0;
    }
 
    if (type_ != event.type_)
    {
        cerr << "!! Error: Invalid type of class used" << endl
             << "Type expected: " << type_ << endl
             << "Event given: " << event.type_ << endl;
        exit(0);
    }
 
    if (!event.config_.is_makeconfig_)
    {
        cout << "Autocall to: DAQEvent::MakeConfig()...";
        event.MakeConfig(*this);
        cout << " Done!" << endl;
    }
 
    DAQFile &file = *this;
    TAG bTag, cTag, tag;
    vector<float> volts(SAMPLES_PER_WAVEFORM);
    int i = 0, j = 0;
 
    // Read only one event
    file >> event.eh_;
    if (event.eh_.serialNumber % 100 == 0)
    {
        cout << "Event serial number: " << event.eh_.serialNumber << endl;
    }
 
    event.is_init_ = true;
    event.routine_ = {false, false, false};
 
    while (file >> bTag)
    {
        file >> tag; // Trigger cell
        auto tCell = *(unsigned short *)(tag.tag + 2);
        j = 0;
        while (file >> cTag)
        {
            if (!is_lab_)
            {
                file >> tag; // Time scaler, LAB-DRS don't have time scaler
            }
            file.Read(volts, event.eh_.rangeCenter);
            event.volts_[i][j] = volts;
            event.TimeCalibration(tCell, times_[i][j], i, j);
            ++j;
        }
        ++i;
        file.ResetTag();
    }
    file.ResetTag();
    return 1;
}
 
bool DAQFile::operator>>(WDBEvent &event) // DAQFile >> WDBEvent
{
    if (!in_.good())
    {
        return 0;
    }
 
    if (type_ != event.type_)
    {
        cerr << "!! Error: Invalid type of class used" << endl
             << "Type expected: " << type_ << endl
             << "Event given: " << event.type_ << endl;
        exit(0);
    }
 
    if (!event.config_.is_makeconfig_)
    {
        cout << "Autocall to: DAQEvent::MakeConfig()...";
        event.MakeConfig(*this);
        cout << " Done!" << endl;
    }
 
    DAQFile &file = *this;
    TAG bTag, cTag, tag;
    vector<float> volts(SAMPLES_PER_WAVEFORM);
    int i = 0, j = 0;
 
    // Read only one event
    file >> event.eh_;
    if (event.eh_.serialNumber % 100 == 0 and event.eh_.serialNumber > 0)
    {
        cout << "Event serial number: " << event.eh_.serialNumber << endl;
    }
 
    event.is_init_ = true;
    event.routine_ = {false, false, false};
 
    while (file >> bTag)
    {
        j = 0;
        while (file >> cTag)
        {
            file >> tag; // Time scaler
            file >> tag; // Trigger cell
            auto tCell = *(unsigned short *)(tag.tag + 2);
            file.Read(volts, event.eh_.rangeCenter);
            event.volts_[i][j] = volts;
            event.TimeCalibration(tCell, times_[i][j], i, j);
            ++j;
        }
        file.ResetTag();
        ++i;
    }
    file.ResetTag();
    return 1;
}
 
void DAQFile::Read(TAG &t)
{
    in_.read(t.tag, 4);
    n_ = t.tag[0];
    return;
}
 
void DAQFile::Read(EventHeader &eh)
{
    in_.read((char *)&eh, sizeof(eh));
    n_ = eh.tag[0];
}
 
void DAQFile::Read(vector<float> &vec)
{
    for (int i = 0; i < vec.size(); ++i)
    {
        in_.read((char *)&vec.at(i), sizeof(float));
    }
    return;
}
 
void DAQFile::Read(vector<float> &vec, const unsigned short &range_center)
{
    unsigned short val;
    for (int i = 0; i < vec.size(); ++i)
    {
        in_.read((char *)&val, sizeof(unsigned short));
        vec[i] = val / 65536. + range_center / 1000. - 0.5;
    }
    return;
}
 
DAQFile::operator bool()
{
    map<char, char> header{{'E', 'B'}, {'B', 'C'}, {'C', 'B'}};
    if (!in_.good())
    {
        cout << "End of file reached" << endl;
        return 0;
    }
    else if (n_ == 'T' or n_ == 'D') // Ignores DRSx and TIME
        return 0;
    else if (n_ == o_) // C --> C, B --> B
        return 1;
    else if (n_ == header[o_]) // E --> B, B --> C, C --> B
    {
        o_ = n_;
        return 1;
    }
    else if (n_ == 'E' and initialization_ == 0) // End of initialization, first EHDR
    {
        cout << "Initialization done --> EHDR next" << endl;
        initialization_ = 1;
        return 0;
    }
    return 0;
}