BmnScWallRaw2Digit.cxx

Go to the documentation of this file.

#include "BmnScWallRaw2Digit.h"
 
#include "BmnADCDigit.h"
#include "BmnScWallDigi.h"
#include "BmnSyncDigit.h"
#include "TMath.h"
#include "TSystem.h"
#include "UniDetectorParameter.h"
 
#include <TStopwatch.h>
#include <bitset>
 
void BmnScWallRaw2Digit::print()
{
    LOG(info) << Form("%s", GetName().Data());
}
 
BmnScWallRaw2Digit::BmnScWallRaw2Digit()
    : WfmProcessor()
 
    , fworkTime_real(0.0)
    , fworkTime_cpu(0.0)
{}
 
BmnScWallRaw2Digit::BmnScWallRaw2Digit(int period,
                                       int run,
                                       const std::string& mappingFile,
                                       const std::string& calibrationFile)
    : WfmProcessor()
 
    , fworkTime_real(0.0)
    , fworkTime_cpu(0.0)
{
    auto logIfError = [this](const char* what, BmnStatus st) {
        if (st == kBMNERROR)
            LOG(error) << GetName() << ": " << what << " was not loaded, decoder may misbehave!";
    };
 
    LOG(info) << Form("%s : Config file: %s", GetName().Data(), mappingFile.c_str());
    logIfError("mapping", ParseConfig(mappingFile));
 
    BmnStatus calibStatus;
    if (!calibrationFile.empty()) {
        LOG(info) << Form("%s : Hand-set Calibration file: %s", GetName().Data(), calibrationFile.c_str());
        calibStatus = ParseCalibration(calibrationFile);
    } else {
        calibStatus = LoadCalibration(period, run);
    }
    logIfError("calibration", calibStatus);
}
 
BmnStatus BmnScWallRaw2Digit::ParseConfig(const std::string& file)
{
    std::string fullPath = std::string(getenv("VMCWORKDIR")) + "/parameters/scwall/" + file;
    std::ifstream config_file(fullPath);
    if (!config_file.is_open()) {
        LOG(error) << GetName() << "::ParseConfig: Cannot open config file: " << fullPath;
        return kBMNERROR;
    }
 
    nlohmann::json j;
    try {
        config_file >> j;
    } catch (const std::exception& e) {
        LOG(error) << GetName() << "::ParseConfig: Failed to parse JSON: " << e.what();
        return kBMNERROR;
    }
 
    LOG(debug) << GetName() << "::ParseConfig. Version: " << j.value("version", "unknown");
    LOG(debug) << GetName() << "::ParseConfig. Comment: " << j.value("comment", "none");
 
    fuoChannelMap.clear();
    for (const auto& entry : j["mapping"]) {
        std::string serial = entry["serial"];
        int channel = entry["channel"];
        int cell = entry["cell"];
 
        auto key = std::make_pair(std::stoul(serial, nullptr, 16), channel);
        fuoChannelMap[key] = BmnScWallAddress::GetAddress(cell);
    }
 
    return kBMNSUCCESS;
}
 
BmnStatus BmnScWallRaw2Digit::LoadCalibration(int period, int run)
{
    nlohmann::json j;
    bool gotFromDb = false;
 
    UniDetectorParameter* pParam =
        UniDetectorParameter::GetDetectorParameter("SCWALL", "scwall_calibration", period, run);
    if (pParam != nullptr) {
        JsonValue* pValue = (JsonValue*)pParam->GetValue();
        if (pValue != nullptr) {
            j = pValue->value;
            gotFromDb = true;
            delete pValue;
        }
        delete pParam;
    }
 
    if (gotFromDb) {
        LOG(info) << GetName() << " : Loaded calibration from DB for period " << period << ", run " << run;
 
        return ParseCalibrationJson(j);
    } else {
        const std::string fallbackFile = Form("SCWALL_calibration_period%d.json", period);
 
        LOG(warn) << GetName() << " : No calibration in DB for period " << period << ", run " << run
                  << ". Falling back to file " << fallbackFile;
 
        return ParseCalibration(fallbackFile);
    }
}
 
BmnStatus BmnScWallRaw2Digit::ParseCalibration(const std::string& file)
{
    std::string fullPath = std::string(getenv("VMCWORKDIR")) + "/parameters/scwall/" + file;
    std::ifstream config_file(fullPath);
    if (!config_file.is_open()) {
        LOG(error) << GetName() << "::ParseCalibration: Cannot open config file: " << fullPath;
        return kBMNERROR;
    }
 
    nlohmann::json j;
    try {
        config_file >> j;
    } catch (const std::exception& e) {
        LOG(error) << GetName() << "::ParseCalibration: Failed to parse JSON: " << e.what();
        return kBMNERROR;
    }
 
    return ParseCalibrationJson(j);
}
 
BmnStatus BmnScWallRaw2Digit::ParseCalibrationJson(const nlohmann::json& j)
{
    LOG(debug) << GetName() << "::ParseCalibrationJson. Version: " << j.value("version", "unknown");
    LOG(debug) << GetName() << "::ParseCalibrationJson. Comment: " << j.value("comment", "none");
 
    try {
 
        const auto& par = j.at("parameters");
 
        fdigiPars.isWriteWfm = par.value("isWriteWfm", false);
        fdigiPars.gateBegin = par.at("gateBegin");
        fdigiPars.gateEnd = par.at("gateEnd");
        fdigiPars.threshold = par.at("threshold");
        fdigiPars.doInvert = par.value("doInvert", false);
 
        fuoCalibMap.clear();
        for (const auto& ch : j.at("calibration")) {
            int cell = ch.at("cell");
            float calib = ch.at("ampl");
 
            fuoCalibMap[BmnScWallAddress::GetAddress(cell)] = std::make_pair(calib, 0.0);
        }
 
        if (j.contains("fitParameters")) {
            const auto& fitpar = j.at("fitParameters");
            fdigiPars.isfit = fitpar.value("isfit", false);
            fdigiPars.harmonics.clear();
            for (const auto& pair : fitpar.at("harmonics")) {
                float re = pair.at(0).get<float>();
                float im = pair.at(1).get<float>();
                fdigiPars.harmonics.emplace_back(re, im);
            }
 
            if (fdigiPars.isfit) {
                int model_order = fdigiPars.harmonics.size() + 1;
                fSignalLength = fdigiPars.gateEnd - fdigiPars.gateBegin + 1;
                auto maxElement = std::max_element(fdigiPars.harmonics.begin(), fdigiPars.harmonics.end(),
                                                   [](const std::complex<float>& a, const std::complex<float>& b) {
                                                       return std::abs(a) < std::abs(b);
                                                   });
                if (maxElement != fdigiPars.harmonics.end()) {
                    // set signal length as 5 * harmonic length
                    int maxElementLength = floor(5. / (-log(real(*maxElement))));
                    if (maxElementLength < fSignalLength)
                        fSignalLength = maxElementLength;
                }
 
                fAZik = new std::complex<float>*[model_order];
                for (int i = 0; i < model_order; i++) {
                    fAZik[i] = new std::complex<float>[model_order];
                    for (int k = 0; k < model_order; k++)
                        fAZik[i][k] = {0., 0.};
                }
                PsdSignalFitting::PronyFitter Pfitter;
                Pfitter.Initialize(fdigiPars.harmonics.size(), fdigiPars.harmonics.size(), fdigiPars.gateBegin,
                                   fdigiPars.gateEnd);
                Pfitter.SetExternalHarmonics(fdigiPars.harmonics);
                Pfitter.MakeInvHarmoMatrix(fSignalLength, fAZik);
            }
        }
    } catch (const std::exception& e) {
        LOG(error) << GetName() << "::ParseCalibrationJson: error while parsing: " << e.what();
        return kBMNERROR;
    }
 
    return kBMNSUCCESS;
}
 
std::optional<uint32_t> BmnScWallRaw2Digit::GetAddressFromBoard(std::pair<size_t, size_t> key)
{
    try {
        return fuoChannelMap.at(key);
    } catch (const std::out_of_range& e) {
        LOG(debug) << "BmnScWallRaw2Digit:: GetAddressFromBoard. Board serial " << std::hex << key.first << std::dec
                   << " channel " << key.second << " is not connected. Skipping it";
        return std::nullopt;
    }
}
 
std::optional<std::pair<float, float>> BmnScWallRaw2Digit::GetCalibPairFromAddress(uint32_t address)
{
    try {
        return fuoCalibMap.at(address);
    } catch (const std::out_of_range& e) {
        LOG(debug) << "BmnScWallRaw2Digit:: GetCalibPairFromAddress " << std::bitset<32>(address)
                   << " is marked badly in calibration file.";
        return std::nullopt;
    }
}
void BmnScWallRaw2Digit::FillEvent(TClonesArray* data, TClonesArray* ScWalldigit)
{
 
    LOG(debug) << "BmnScWallRaw2Digit::fillEvent";
    TStopwatch timer;
    timer.Start();
 
    for (int i = 0; i < data->GetEntriesFast(); i++) {
        BmnADCDigit<int16_t>* digit = (BmnADCDigit<int16_t>*)data->At(i);
        auto key = std::make_pair(digit->GetSerial(), digit->GetChannel());
        auto catch_address = GetAddressFromBoard(key);
        if (!catch_address.has_value())
            continue;   // Not connected lines
        auto address = catch_address.value();
        if (address == 0)
            continue;   // Not connected lines
 
        auto catch_calib = GetCalibPairFromAddress(address);
        if (!catch_calib.has_value())
            continue;   // not involved in the analysis
        auto calib_pair = catch_calib.value();
 
        std::vector<float> wfm(digit->GetValue().begin(), digit->GetValue().end());
        BmnScWallDigi ThisDigi;
        ThisDigi.reset();
        ThisDigi.SetAddress(address);
        ProcessWfm(wfm, &ThisDigi);
 
        // Apply calibration
        LOG(debug4) << "BmnScWallRaw2Digit::ProcessWfm  Calibration";
        if (fdigiPars.signalType == 0)
            ThisDigi.SetSignal((float)ThisDigi.fAmpl * calib_pair.first);
        if (fdigiPars.signalType == 1)
            ThisDigi.SetSignal((float)ThisDigi.fIntegral * calib_pair.first);
        if (ThisDigi.GetSignal() < fdigiPars.threshold)
            continue;
 
        TClonesArray& ar_ScWall = *ScWalldigit;
        new (ar_ScWall[ScWalldigit->GetEntriesFast()]) BmnScWallDigi(ThisDigi);
    }
 
    timer.Stop();
    fworkTime_cpu += timer.RealTime();
    fworkTime_real += timer.CpuTime();
}
 
BmnScWallRaw2Digit::~BmnScWallRaw2Digit()
{
    LOG(info) << Form("%10s: Real time %4.3f s, CPU time %4.3f s", GetName().Data(), fworkTime_real, fworkTime_cpu);
}