BmnHodoDigitizer.cxx

Go to the documentation of this file.

#include "BmnHodoDigitizer.h"
 
#include "BmnHodoAddress.h"
#include "FairRootManager.h"
#include "TMath.h"
#include "TStopwatch.h"
 
#include <FairRun.h>
#include <FairRunSim.h>
#include <cstdlib>
#include <fstream>
#include <nlohmann/json.hpp>
#include <stdexcept>
#include <string>
 
BmnHodoDigitizer::BmnHodoDigitizer(int period)
    : BmnHodoDigitizer(Form("HODO_digitizer_period%d.json", period))
{
    LoadCalibration(Form("HODO_Q_calibration_period%d.json", period));
}
 
BmnHodoDigitizer::BmnHodoDigitizer(const char* config)
    : FairTask("BmnHodoDigitizer")
    , fScale(1337.8)    // [Z^2/GeV] set Xe peak to 54*54
    , fThreshold(0.0)   // [Z^2] we see noise below this in data
    , fTimeCut(1000.0)
    , fSaturation(std::numeric_limits<double>::max())
    , fPointArray(nullptr)
    , fDigiArray(nullptr)
{
    if (config && config[0] != '\0') {
        LoadConfig(config);
    }
}
 
InitStatus BmnHodoDigitizer::LoadConfig(const char* config)
{
    const char* vmc = getenv("VMCWORKDIR");
    if (!vmc || !vmc[0]) {
        LOG(error) << Form("%s::LoadConfig: VMCWORKDIR is not set", GetName());
        SetActive(kFALSE);
        return kERROR;
    }
 
    const std::string fullPath = std::string(vmc) + "/parameters/hodo/" + config;
    std::ifstream input(fullPath);
    if (!input.is_open()) {
        LOG(error) << Form("%s::LoadConfig: cannot open file %s", GetName(), fullPath.c_str());
        SetActive(kFALSE);
        return kERROR;
    }
 
    nlohmann::json j;
    try {
        input >> j;
    } catch (const std::exception& e) {
        LOG(error) << Form("%s::LoadConfig: failed to parse JSON %s: %s", GetName(), fullPath.c_str(), e.what());
        SetActive(kFALSE);
        return kERROR;
    }
 
    const std::string version = j.value("version", "unknown");
    const std::string comment = j.value("comment", "none");
 
    if (j.contains("scale"))
        fScale = j["scale"].get<double>();
 
    if (j.contains("threshold"))
        fThreshold = j["threshold"].get<double>();
 
    if (j.contains("timecut"))
        fTimeCut = j["timecut"].get<double>();
 
    LOG(info) << Form("%s::LoadConfig: loaded %s (version=%s, comment=%s)", GetName(), fullPath.c_str(),
                      version.c_str(), comment.c_str());
 
    return kSUCCESS;
}
 
// Destructor
BmnHodoDigitizer::~BmnHodoDigitizer()
{
    if (fDigiArray) {
        fDigiArray->Delete();
        delete fDigiArray;
    }
}
 
// Initialization method
InitStatus BmnHodoDigitizer::Init()
{
    fworkTime = 0.0;
    LOG(detail) << "-I- BmnHodoDigitizer: Init started..." << std::endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << "-E- BmnHodoDigitizer::Init: RootManager not instantiated!" << std::endl;
        return kFATAL;
    }
 
    fPointArray = static_cast<TClonesArray*>(ioman->GetObject("HodoPoint"));
    if (!fPointArray) {
        LOG(error) << "-W- BmnHodoDigitizer::Init: No HodoPoint array!" << std::endl;
        return kERROR;
    }
 
    fDigiArray = new TClonesArray("BmnHodoDigit");
    ioman->Register("HodoDigit", "Hodo", fDigiArray, kTRUE);
 
    LOG(detail) << "-I- BmnHodoDigitizer: Initialization successful" << std::endl;
    return kSUCCESS;
}
 
// Setter for energy resolution configuration
void BmnHodoDigitizer::LoadCalibration(const char* inputFile)
{
    std::string fullPath = std::string(getenv("VMCWORKDIR")) + "/parameters/hodo/" + inputFile;
    std::ifstream input(fullPath);
    if (!input.is_open()) {
        LOG(error) << "BmnHodoDigitizer::LoadCalibration: Cannot open file " << fullPath;
        return;
    }
 
    nlohmann::json j;
    try {
        input >> j;
    } catch (const std::exception& e) {
        LOG(error) << "BmnHodoDigitizer::LoadCalibration: Failed to parse JSON: " << e.what();
        return;
    }
 
    std::string version = j.value("version", "unknown");
    std::string comment = j.value("comment", "none");
 
    LOG(debug) << "BmnHodoDigitizer::LoadCalibration. Version: " << version;
    LOG(debug) << "BmnHodoDigitizer::LoadCalibration. Comment: " << comment;
 
    fSaturation = j.at("parameters").at("saturation");
 
    fuoResolutionMap.clear();
    for (const auto& entry : j["calibration"]) {
        int strip = entry["strip"];
        int side = entry["side"];
        int gain = entry["gain"];
        double calib = entry["calib"];
        double stochastic = entry["stochastic"];
        double constant = entry["constant"];
 
        uint32_t address = BmnHodoAddress::GetAddress(strip, side, gain);
        Kernel resolution(stochastic, constant);
        fuoResolutionMap.emplace(address, std::make_pair(calib, resolution));
    }
}
 
// Execute method (digitization process)
void BmnHodoDigitizer::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
    TStopwatch sw;
    sw.Start();
 
    fDigiArray->Clear();
    FillHitMap();
 
    for (const auto& it : fuoHitMap) {
        uint32_t strip = it.first;
        const auto& points = it.second;
 
        double energyLoss = 0.0;
        double time = 0.0;
        for (const auto* point : points) {
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug))
                point->Print("");
            energyLoss += point->GetEnergyLoss();
            time += point->GetTime();
        }
        time /= points.size();
 
        if (fuoResolutionMap.empty()) {
            uint32_t address = BmnHodoAddress::GetAddress(strip, 0, 0);
            double z2 = energyLoss * fScale;
            LOG(debug) << "BmnHodoDigitizer::Exec Empty resoultion map. Smoothed energy: " << z2;
 
            if (z2 < fThreshold)
                continue;
            else {
                TClonesArray& ar = *fDigiArray;
                new (ar[ar.GetEntriesFast()]) BmnHodoDigit(address, time, z2);
            }
        } else {
            for (const auto& jt : fuoResolutionMap) {
                uint32_t address = jt.first;
                if (BmnHodoAddress::GetStripId(address) != strip)
                    continue;
                double z2 = jt.second.second(energyLoss * fScale);
                double adc = z2 / jt.second.first;
                LOG(debug) << Form(
                    "BmnHodoDigitizer::Exec Resoultion map for strid %d side %d gain %d. Smoothed energy: %f",
                    BmnHodoAddress::GetStripId(address), BmnHodoAddress::GetStripSide(address),
                    BmnHodoAddress::GetGain(address), z2);
 
                if (z2 < fThreshold)
                    continue;
                else {
                    TClonesArray& ar = *fDigiArray;
                    BmnHodoDigit* digit = new (ar[ar.GetEntriesFast()]) BmnHodoDigit(address, time, z2);
                    if (adc >= fSaturation) {
                        digit->SetSignal(fSaturation * jt.second.first);
                        digit->SetIsSaturated(true);
                    }
                }
            }
        }
    }
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
// Fill the hit map with points
void BmnHodoDigitizer::FillHitMap()
{
    fuoHitMap.clear();
 
    for (int i = 0; i < fPointArray->GetEntriesFast(); ++i) {
        auto* point = static_cast<BmnHodoPoint*>(fPointArray->At(i));
        int iStrip = point->GetCopyMother();
        fuoHitMap[iStrip].emplace_back(point);
    }
 
    LOG(debug) << "BmnHodoDigitizer::FillHitMap() - Number of fired cells: " << fuoHitMap.size();
}
 
// Finish method
void BmnHodoDigitizer::Finish()
{
    printf("Work time of the Hodo digitizer: %.4f sec.\n", fworkTime);
}