BmnTrigRaw2Digit.cxx

Go to the documentation of this file.

#include "BmnTrigRaw2Digit.h"
// C++
#include <climits>
// Boost
#include <boost/program_options.hpp>
// ROOT
#include <TStopwatch.h>
// FairRoot
#include <FairLogger.h>
// BmnRoot
#include "BmnTDCDigit.h"
#include "BmnTQDCADCDigit.h"
#include "BmnTrigDigit.h"
#include "BmnTrigWaveDigit.h"
 
namespace po = boost::program_options;
using namespace std;
using namespace TMath;
 
BmnTrigParameters::BmnTrigParameters()
{
    for (UInt_t i = 0; i < CHANNEL_COUNT_MAX; i++) {
        for (UInt_t j = 0; j < TDC_BIN_COUNT; j++)
            INL[i][j] = 0u;
        ChannelMap[i] = 0u;
        NegativeMap[i] = kFALSE;
        branchArrayPtr[i] = NULL;
        t[i] = 0.0;
    }
    IsT0 = kFALSE;
}
 
BmnTrigRaw2Digit::BmnTrigRaw2Digit(TString PlacementMapFile, TString StripMapFile, TTree* digiTree)
    : fSetup(kBMNSETUP)
    //    , fWorkArray{}
    , fTrcSerial(0)
{
    //    fWorkHist->SetDirectory(0);
    ReadPlacementMap(PlacementMapFile);
    if (ReadChannelMap(StripMapFile) == kBMNERROR) {
        printf("Error! Unable to read the trigger channel map!\n");
        return;
    }
    if (digiTree == NULL)
        return;
    // Create corresponding branches for each trigger.       //
    // different channels of the one trigger will be stored //
    // in one branch with different mod ID  //
    map<TString, TClonesArray*> brName2ptr;
    for (BmnTrigChannelData& record : fMap) {
        TString& detName = record.name;
        TClass* cl = detName.Contains("TQDC") ? BmnTrigWaveDigit::Class() : BmnTrigDigit::Class();
        auto it = brName2ptr.find(detName);
        if (it == brName2ptr.end()) {
            TClonesArray* ar = new TClonesArray(cl);
            ar->SetName(detName.Data());
            digiTree->Branch(detName.Data(), &ar);
            trigArrays.push_back(ar);
            record.branchArrayPtr = ar;
            brName2ptr.insert(make_pair(detName, ar));
        } else {
            record.branchArrayPtr = it->second;
        }
    }
    // Fill elements of placement map with channel->(strip, mod, branchRef) map //
    for (BmnTrigChannelData& record : fMap) {
        BmnTrigParameters* par = nullptr;
        map<PlMapKey, BmnTrigParameters*>::iterator itPar = fPlacementMap.find(PlMapKey(record.serial, record.slot));
        if (itPar == fPlacementMap.end()) {
            printf("CrateSeral %08X slot %u not found in the placement map!\n", record.serial, record.slot);
            par = new BmnTrigParameters();
            par->BoardSerial = record.serial;
            par->name = record.name;
            par->ChannelCount = CHANNEL_COUNT_MAX;   // ChanCntByName(record.name);
            fPlacementMap.insert(pair<PlMapKey, BmnTrigParameters*>(PlMapKey(par->BoardSerial, par->slot), par));
        } else
            par = itPar->second;
        par->ChannelMap[record.channel] = record.module;
        par->branchArrayPtr[record.channel] = record.branchArrayPtr;
        par->NegativeMap[record.channel] = record.isNegative;
    }
    for (auto& el : fPlacementMap) {
        BmnTrigParameters* par = el.second;
        ReadINLFromFile(par);
    }
}
 
// BmnTrigRaw2Digit::BmnTrigRaw2Digit(TString PlacementMapFile, TString StripMapFile, FairRootManager* frm, bool
// saveOutput) :
 
BmnTrigRaw2Digit::BmnTrigRaw2Digit(TString PlacementMapFile,
                                   TString StripMapFile,
                                   std::function<void(TString, TClonesArray*)> BranchRegFun)
    : fSetup(kBMNSETUP)
    //    , fWorkHist(make_unique<TH1I>(Form("workhist_%p", this), "wh", SHRT_MAX - SHRT_MIN, SHRT_MIN, SHRT_MAX))
    , fTrcSerial(0)
{
    //    fWorkHist->SetDirectory(0);
    ReadPlacementMap(PlacementMapFile);
    if (ReadChannelMap(StripMapFile) == kBMNERROR) {
        printf("Error! Unable to read the trigger channel map!\n");
        return;
    }
    // Create corresponding branches for each trigger.       //
    // different channels of the one trigger will be stored //
    // in one branch with different mod ID  //
    map<TString, TClonesArray*> brName2ptr;
    for (BmnTrigChannelData& record : fMap) {
        TString& detName = record.name;
        TClass* cl = detName.Contains("TQDC") ? BmnTrigWaveDigit::Class() : BmnTrigDigit::Class();
        auto it = brName2ptr.find(detName);
        if (it == brName2ptr.end()) {
            TClonesArray* ar = new TClonesArray(cl);
            ar->SetName(detName.Data());
            BranchRegFun(detName, ar);
            //            frm->Register(detName, "Trig", ar, saveOutput);
            trigArrays.push_back(ar);
            record.branchArrayPtr = ar;
            brName2ptr.insert(make_pair(detName, ar));
        } else
            record.branchArrayPtr = it->second;
    }
    // Fill elements of placement map with channel->(strip, mod, branchRef) map //
    for (BmnTrigChannelData& record : fMap) {
        BmnTrigParameters* par = nullptr;
        map<PlMapKey, BmnTrigParameters*>::iterator itPar = fPlacementMap.find(PlMapKey(record.serial, record.slot));
        if (itPar == fPlacementMap.end()) {
            printf("CrateSeral %08X slot %u not found in the placement map!\n", record.serial, record.slot);
            par = new BmnTrigParameters();
            par->BoardSerial = record.serial;
            par->name = record.name;
            par->ChannelCount = CHANNEL_COUNT_MAX;   // ChanCntByName(record.name);
            fPlacementMap.insert(pair<PlMapKey, BmnTrigParameters*>(PlMapKey(par->BoardSerial, par->slot), par));
        } else
            par = itPar->second;
        par->ChannelMap[record.channel] = record.module;
        par->branchArrayPtr[record.channel] = record.branchArrayPtr;
        par->NegativeMap[record.channel] = record.isNegative;
    }
    for (auto& el : fPlacementMap) {
        BmnTrigParameters* par = el.second;
        ReadINLFromFile(par);
    }
}
 
BmnStatus BmnTrigRaw2Digit::ReadPlacementMap(TString mappingFile)
{
    TString PlMapFileName = TString(getenv("VMCWORKDIR")) + TString("/input/") + mappingFile;
    LOGF(info, "Reading Triggers placement mapping file %s...", PlMapFileName.Data());
    ifstream pmFile;
    pmFile.open(PlMapFileName.Data());
    if (!pmFile.is_open()) {
        cout << "Error opening map-file (" << PlMapFileName << ")!" << endl;
        return kBMNERROR;
    }
    string dummy;
    string name;
    UInt_t crateSerial, boardSerial;
    uint16_t slot;
    uint16_t isT0;
 
    pmFile >> dummy >> dummy >> dummy >> dummy >> dummy;
    pmFile >> dummy;
    while (pmFile >> name >> hex >> crateSerial >> dec >> slot >> hex >> boardSerial >> dec >> isT0) {
        BmnTrigParameters* par = new BmnTrigParameters();
        par->BoardSerial = boardSerial;
        par->CrateSerial = crateSerial;
        par->slot = slot;
        par->name = name;
        par->IsT0 = isT0;
        par->ChannelCount = CHANNEL_COUNT_MAX;   // ChanCntByName(name);
        fPlacementMap.insert(pair<PlMapKey, BmnTrigParameters*>(PlMapKey(par->CrateSerial, par->slot), par));
    }
    pmFile.close();
    return kBMNSUCCESS;
}
 
BmnStatus BmnTrigRaw2Digit::ReadChannelMap(TString mappingFile)
{
    fMapFileName = TString(getenv("VMCWORKDIR")) + TString("/input/") + mappingFile;
    LOGF(info, "Reading Triggers strip mapping file %s...", fMapFileName.Data());
 
    vector<string> trc_trigger_bits;
    vector<string> scalers;
    vector<string> configuration;
    vector<string> t0_raw;
 
    // Setup options.
    po::options_description desc("Options");
    desc.add_options()("INPUTS.serial", po::value<string>(), "TRC serial")(
        "INPUTS.channels", po::value<vector<string>>(&trc_trigger_bits)->multitoken(), "TRC trigger input signal bits")(
        "SCALERS.channels", po::value<vector<string>>(&scalers)->multitoken(),
        "Scaler channels")("TDC-CONFIGURATION.channels", po::value<vector<string>>(&configuration)->multitoken(),
                           "TDC configuration")("T0-RAW.virt_serial", po::value<string>(), "T0 virtual serial")(
        "T0-RAW.channels", po::value<vector<string>>(&t0_raw)->multitoken(), "TDC configuration");
 
    // Load config file.
    po::variables_map vm;
    ifstream config_file(fMapFileName.Data(), ifstream::in);
    if (!config_file.is_open()) {
        printf("%s - file open error!\n", mappingFile.Data());
        return kBMNERROR;
    }
    po::store(po::parse_config_file(config_file, desc), vm);
    config_file.close();
    po::notify(vm);
 
    stringstream hs(vm["INPUTS.serial"].as<string>());
    hs >> std::hex >> fTrcSerial;
    LOGF(info, "TRC serial %8X", fTrcSerial);
 
    //========== read TRC trig bits map ==========//
    string tr_name;
    uint16_t bit;
    for (auto it : trc_trigger_bits) {
        istringstream ss(it);
        ss >> tr_name >> bit;
        fTrcBitMap[tr_name] = bit;
        //        printf("trig[%8s] = %u\n", tr_name.data(), fTrigBitsMap[tr_name]);
    }
    //========== read scalers channel map ==========//
    UInt_t ser;
    uint16_t ch;
    uint16_t scalers_ar_size(0);   // MSC list size
    for (auto it : scalers) {
        istringstream ss(it);
        ss >> tr_name >> hex >> ser >> dec >> ch;
        std::pair p_ser_ch = make_pair(ser, ch);
        scalers2name_map.insert(make_pair(p_ser_ch, tr_name));
        scalers2index_map.insert(make_pair(p_ser_ch, scalers_ar_size));
        scalers_name2index_map.insert(make_pair(tr_name, scalers_ar_size));
        scalers_index2name_map.push_back(tr_name);
        //        scalers_index2name_map.insert(make_pair(scalers_ar_size, tr_name));
        name2scaler_map[tr_name] = p_ser_ch;
        auto it_trc = fTrcBitMap.find(tr_name);
        if (it_trc != fTrcBitMap.end())
            trcIdx2scalerIdx.insert(make_pair(it_trc->second, scalers_ar_size));
        scalers_ar_size++;
        LOGF(debug, "msc[%15s] = (0x%08X : %2u)", tr_name.data(), name2scaler_map[tr_name].first,
             name2scaler_map[tr_name].second);
    }
    //========== read TDC channel map ==========//
    TString name;
    uint16_t slot, mod;
    uint16_t neg;
    for (auto it : configuration) {
        istringstream ss(it);
        ss >> name >> mod >> hex >> ser >> dec >> slot >> ch >> neg;
        if (ch < 0)
            continue;
        BmnTrigChannelData record;
        record.branchArrayPtr = NULL;
        record.name = name;
        record.serial = ser;
        record.module = mod;
        record.slot = slot;
        record.channel = ch;
        record.isNegative = neg > 0 ? kTRUE : kFALSE;
        fMap.push_back(record);
    }
    //========== read T0 raw signals from Sergeev's txt ==========//
    stringstream vs(vm["T0-RAW.virt_serial"].as<string>());
    UInt_t t0_virt_serial;
    vs >> std::hex >> t0_virt_serial;
    LOGF(info, "t0_virt_serial %0X", t0_virt_serial);
    for (auto it : t0_raw) {
        istringstream ss(it);
        ss >> tr_name >> mod >> ch;
        fT0Map[make_pair(tr_name, mod)] = ch;
        LOGF(debug, "trig[%10s] mod %2u  ch %3u", tr_name.data(), mod, ch);
    }
    return kBMNSUCCESS;
}
 
BmnStatus BmnTrigRaw2Digit::ReadINLFromFile(BmnTrigParameters* par)
{
    fstream ff;
    TString fINLFileName =
        TString(getenv("VMCWORKDIR")) + TString("/input/") + Serial2FileName(par->name, par->BoardSerial);
    ff.open((fINLFileName).Data(), ios::in);
    if (!ff.is_open()) {
        cout << "Error opening INL-file (" << fINLFileName << ")!" << endl;
        return kBMNERROR;
    }
    LOGF(info, "Open INL file %s", fINLFileName.Data());
    TPRegexp reInlHdr("\\[.*(inl_corr).*\\]");   // INL block header
    //    regex reInlHdr("\\‍[.*(inl_corr).*\\‍]"); // INL block header
    //    regex reInlChannel("\\s*(\\d+)=(.+)"); // chID=c0, c1, ...
    Bool_t isInlHdr = kFALSE;
    while (!ff.eof()) {
        string line;
        std::getline(ff, line, '\n');
        if (reInlHdr.MatchB(line)) {
            //        if (regex_match(line, reInlHdr)) {
            isInlHdr = kTRUE;
            break;
        }
    }
    if (!isInlHdr) {
        printf("Incorrect INL file format!\n");
        return kBMNERROR;
    }
    uint16_t channelID = 0;
    while (!ff.eof()) {
        string line;
        std::getline(ff, line, '\n');
        //        printf("Read %lu \n", line.length());
        //        if (!regex_match(line, reInlChannel))
        //            continue;
        //        line = regex_replace(line, reInlChannel, "$1 $2");
        //        printf("%s\n", line.c_str());
        istringstream ss(line);
        ss >> channelID;
        // printf("Channel ID = %u\n", channelID);
        uint16_t i_bin = 0;
        while (ss.tellg() != -1) {
            if (i_bin > TDC_BIN_COUNT) {
                printf("INL File contains too many bins in channel.\n");
                ff.close();
                return kBMNERROR;
            }
            if (ss.peek() == ',' || ss.peek() == '=') {
                ss.ignore();
            }
            ss >> par->INL[channelID][i_bin];
            i_bin++;
        }
    }
    ff.close();
    return kBMNSUCCESS;
}
 
void BmnTrigRaw2Digit::ProcessWave(vector<int16_t>& iValue, Bool_t& isNeg)
{
    //    Double_t rtime;
    //    Double_t ctime;
    //    TStopwatch timer;
    //    timer.Start();
    int16_t baseLine = GetBaselineByMap(iValue);
    //    timer.Stop();
    //    rtime = timer.RealTime();
    //    ctime = timer.CpuTime();
    //    printf("Real time %f s, CPU time %f s baseline\n", rtime, ctime);
    //    timer.Start();
    for (auto& val : iValue)
        val -= baseLine;
    if (isNeg)
        for (auto& val : iValue)
            val = -val;
    //    timer.Stop();
    //    rtime = timer.RealTime();
    //    ctime = timer.CpuTime();
    //    printf("Real time %f s, CPU time %f s sub\n", rtime, ctime);
}
 
// int16_t BmnTrigRaw2Digit::GetBaseline(vector<int16_t>& iValue)
//{
//     fWorkArray = {};
//     for (auto& val : iValue)
//         fWorkArray[val - SHRT_MIN]++;
//     uint16_t maxCount = 0;
//     size_t mode = 0;
//     for (size_t i = 0; i < ADC_RANGE; ++i)
//         if (fWorkArray[i] > maxCount) {
//             maxCount = fWorkArray[i];
//             mode = i;
//         }
//     return (mode + SHRT_MIN);
// }
 
int16_t BmnTrigRaw2Digit::GetBaselineByMap(vector<int16_t>& iValue)
{
    fValsMap.clear();
    for (auto& val : iValue)
        fValsMap[val]++;
    uint16_t maxCount = 0;
    int16_t mode = 0;
    for (auto& [val, cnt] : fValsMap)
        if (cnt > maxCount) {
            maxCount = cnt;
            mode = val;
        }
    return mode;
}
 
BmnStatus BmnTrigRaw2Digit::FillEvent(TClonesArray* tdc, TClonesArray* adc, unordered_map<UInt_t, Long64_t>& tsMap)
{
    std::vector<Double_t> times;
    std::vector<Double_t> diff;
    for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); iAdc++) {
        times.clear();
        diff.clear();
        BmnTQDCADCDigit<int16_t>* adcDig = (BmnTQDCADCDigit<int16_t>*)adc->At(iAdc);
        uint16_t iChannel = adcDig->GetChannel();
        auto plIter = fPlacementMap.find(PlMapKey(adcDig->GetSerial(), adcDig->GetSlot()));
        if (plIter == fPlacementMap.end())
            continue;
        BmnTrigParameters* par = plIter->second;
        TClonesArray* trigAr = par->branchArrayPtr[iChannel];
        if (trigAr) {
            uint16_t iMod = par->ChannelMap[iChannel];
            Double_t adcTimestamp = adcDig->GetAdcTimestamp() * ADC_CLOCK_TQDC16VS;
            Double_t trgTimestamp = adcDig->GetTrigTimestamp() * ADC_CLOCK_TQDC16VS;
            Double_t boardShift = par->IsT0 ? 0 : tsMap[adcDig->GetSerial()];
 
            for (Int_t iTdc = 0; iTdc < tdc->GetEntriesFast(); ++iTdc) {
                BmnTDCDigit* tdcDig = (BmnTDCDigit*)tdc->At(iTdc);
                if (tdcDig->GetSerial() != adcDig->GetSerial() || tdcDig->GetSlot() != adcDig->GetSlot())
                    continue;
                if (tdcDig->GetChannel() != iChannel)
                    continue;
                Double_t time = (tdcDig->GetValue() + par->INL[iChannel][tdcDig->GetValue() % TDC_BIN_COUNT])
                                    * TDC_CLOCK / TDC_BIN_COUNT
                                + boardShift;
                times.push_back(time);
            }
            BmnTrigWaveDigit* dig = new ((*trigAr)[trigAr->GetEntriesFast()])
                BmnTrigWaveDigit(iMod, move(adcDig->GetValue()), trgTimestamp, adcTimestamp);
            ProcessWave(dig->GetValue(), par->NegativeMap[iChannel]);
            dig->TdcVector() = std::move(times);
        }
    }
    return kBMNSUCCESS;
}
 
BmnStatus BmnTrigRaw2Digit::FillEvent(TClonesArray* tdc, unordered_map<UInt_t, Long64_t>& tsMap)
{
    for (auto& el : fPlacementMap)
        for (UInt_t i = 0; i < CHANNEL_COUNT_MAX; i++)
            el.second->t[i] = -1.0;
    for (Int_t iTdc = 0; iTdc < tdc->GetEntriesFast(); ++iTdc) {
        BmnTDCDigit* tdcDig = (BmnTDCDigit*)tdc->At(iTdc);
        auto plIter = fPlacementMap.find(PlMapKey(tdcDig->GetSerial(), tdcDig->GetSlot()));
        if (plIter == fPlacementMap.end())
            continue;
        BmnTrigParameters* par = plIter->second;
        Double_t boardShift = par->IsT0 ? 0 : tsMap[tdcDig->GetSerial()];
        uint16_t rChannel = tdcDig->GetHptdcId() * kNCHANNELS + tdcDig->GetChannel();
        Double_t time =
            (tdcDig->GetValue() + par->INL[rChannel][tdcDig->GetValue() % TDC_BIN_COUNT]) * TDC_CLOCK / TDC_BIN_COUNT
            + boardShift;
        //        printf("\tCrateSeral %08X slot %02u channel %02u  time %+2.2f  leading %d neg %d\n",
        //        tdcDig->GetSerial(), tdcDig->GetSlot(), rChannel, time,
        //        tdcDig->GetLeading(),par->NegativeMap[rChannel]);
        if (tdcDig->GetLeading() ^ par->NegativeMap[rChannel]) {
            par->t[rChannel] = time;
        } else {
            if (time < par->t[rChannel])
                continue;
            if (par->t[rChannel] < 0)
                continue;
            uint16_t iMod = par->ChannelMap[rChannel];
            TClonesArray* trigAr = par->branchArrayPtr[rChannel];
            if (trigAr == NULL)
                continue;
            Double_t tL = par->t[rChannel];
            Double_t tT = time;
            par->t[rChannel] = -1.0;
            new ((*trigAr)[trigAr->GetEntriesFast()]) BmnTrigDigit(iMod, tL, tT - tL);
        }
    }
    return kBMNSUCCESS;
}
 
BmnStatus BmnTrigRaw2Digit::ClearArrays()
{
    for (TClonesArray* ar : trigArrays)
        ar->Delete();
    return kBMNSUCCESS;
}