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))
{}
 
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())
    , fKernel()
    , 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>();
 
    if (j.contains("stochastic") && j.contains("constant"))
        fKernel = BmnHodoResolution::Kernel(j["stochastic"].get<double>(), j["constant"].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;
}
 
// 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();
 
        uint32_t address = BmnHodoAddress::GetAddress(strip, 0, 0);
        double z2 = fKernel(energyLoss * fScale);
        LOG(debug) << "BmnHodoDigitizer::Exec. Smoothed energy: " << z2;
 
        if (z2 < fThreshold)
            continue;
        else {
            TClonesArray& ar = *fDigiArray;
            new (ar[ar.GetEntriesFast()]) BmnHodoDigit(address, time, z2);
        }
    }
 
    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);
    }
 
    // Sort each vector of points by their time
    for (auto& pair : fuoHitMap) {
        std::vector<BmnHodoPoint*>& points = pair.second;
        std::sort(points.begin(), points.end(),
                  [](BmnHodoPoint* a, BmnHodoPoint* b) { return a->GetTime() < b->GetTime(); });
    }
 
    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);
}