BmnFHCalDigitizer.cxx

Go to the documentation of this file.

#include "BmnFHCalDigitizer.h"
 
#include "BmnFHCalAddress.h"
#include "FairRootManager.h"
#include "TStopwatch.h"
#include "TVirtualMC.h"
 
#include <FairRun.h>
#include <FairRunSim.h>
#include <cstdlib>
#include <fstream>
#include <nlohmann/json.hpp>
#include <stdexcept>
#include <string>
 
BmnFHCalDigitizer::BmnFHCalDigitizer(int period)
    : BmnFHCalDigitizer(Form("FHCAL_digitizer_period%d.json", period))
{}
 
BmnFHCalDigitizer::BmnFHCalDigitizer(const char* config)
    : FairTask("BmnFHCalDigitizer")
    , fScale(1.0 / 0.0048)
    , fThreshold(0.0)
    , fTimeCut(1000.0)
    , fSiPM(90000, 40, 0.3, 0.0048)
    , fPointArray(nullptr)
    , fDigiArray(nullptr)
    , fworkTime(0.0f)
{
    if (config && config[0] != '\0') {
        LoadConfig(config);
    }
}
 
InitStatus BmnFHCalDigitizer::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/fhcal/" + 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("sipm") && j["sipm"].is_object()) {
        const auto& jsipm = j["sipm"];
 
        int nPixels = fSiPM.Npixels;
        int mipPixels = fSiPM.pixPerMIP;
        double noise = fSiPM.noiseMIP;
        double mipEnergy = fSiPM.gevPerMIP;
 
        if (jsipm.contains("n_pixels"))
            nPixels = jsipm["n_pixels"].get<int>();
 
        if (jsipm.contains("mip_pixels"))
            mipPixels = jsipm["mip_pixels"].get<int>();
 
        if (jsipm.contains("noise"))
            noise = jsipm["noise"].get<double>();
 
        if (jsipm.contains("mip_energy"))
            mipEnergy = jsipm["mip_energy"].get<double>();
 
        fSiPM = SiPM(nPixels, mipPixels, noise, mipEnergy);
    }
 
    LOG(info) << Form("%s::LoadConfig: loaded %s (version=%s, comment=%s)", GetName(), fullPath.c_str(),
                      version.c_str(), comment.c_str());
 
    return kSUCCESS;
}
 
// Destructor
BmnFHCalDigitizer::~BmnFHCalDigitizer()
{
    if (fDigiArray) {
        fDigiArray->Delete();
        delete fDigiArray;
    }
}
 
// Initialization method
InitStatus BmnFHCalDigitizer::Init()
{
    fworkTime = 0.0;
    LOG(detail) << "-I- BmnFHCalDigitizer: Init started..." << std::endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << "-E- BmnFHCalDigitizer::Init: RootManager not instantiated!" << std::endl;
        return kFATAL;
    }
 
    fPointArray = static_cast<TClonesArray*>(ioman->GetObject("FHCalPoint"));
    if (!fPointArray) {
        LOG(error) << "-E- BmnFHCalDigitizer::Init: No FHCalPoint array!" << std::endl;
        return kERROR;
    }
 
    fDigiArray = new TClonesArray("BmnFHCalDigit");
    ioman->Register("FHCalDigit", "FHCal", fDigiArray, kTRUE);
 
    LOG(detail) << "-I- BmnFHCalDigitizer: Initialization successful" << std::endl;
    return kSUCCESS;
}
 
void BmnFHCalDigitizer::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
    TStopwatch sw;
    sw.Start();
 
    LOG(debug2) << "BmnFHCalDigitizer::Exec() started..." << std::endl;
    if (!fDigiArray)
        Fatal("Exec", "No DigiArray");
 
    fDigiArray->Clear();
    FillHitMap();   // One address (key) is now associated with a sorted vector of tracks which fired it
 
    for (auto& it : fuoHitMap) {
        uint32_t address = it.first;
        const auto& tracks_contributions = it.second;
 
        double totalEnergy = 0.0;
        double crossTime = 0.0;   // time when threshold is crossed
        bool reachedThreshold = false;
        for (const auto* point : tracks_contributions) {
            if (point->GetTime() > fTimeCut)
                break;
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug2))
                point->Print();
 
            totalEnergy += point->GetEnergyLoss();
            if (!reachedThreshold && totalEnergy * fScale > fThreshold) {
                crossTime = point->GetTime();
                reachedThreshold = true;
                LOG(debug2) << Form("Reached threshold. Energy %.4f Weighted time %.4f", totalEnergy, crossTime);
            }
        }
        totalEnergy = fSiPM.ModelResponse(totalEnergy);
        totalEnergy *= fScale;
 
        if (!reachedThreshold || totalEnergy < std::numeric_limits<double>::epsilon()) {
            LOG(debug2) << "BmnFHCalDigitizer::Skip Digit " << BmnFHCalAddress::GetInfoString(address).c_str();
            continue;
        } else {
            TClonesArray& ar = *fDigiArray;
            long entries = fDigiArray->GetEntriesFast();
            new (ar[entries]) BmnFHCalDigit(address, crossTime, totalEnergy);
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug2))
                static_cast<BmnFHCalDigit*>(ar.At(entries))->Print();
        }
    }
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
void BmnFHCalDigitizer::FillHitMap()
{
    fuoHitMap.clear();
 
    // Fill the hit map with points from the TClonesArray
    for (int iPoint = 0, nPoint = fPointArray->GetEntriesFast(); iPoint < nPoint; ++iPoint) {
        auto* point = static_cast<BmnFHCalPoint*>(fPointArray->At(iPoint));
        uint32_t address = BmnFHCalAddress::GetPhysicalAddress(point->GetAddress());
        fuoHitMap[address].push_back(point);
    }
 
    // Sort each vector of points by their time
    for (auto& pair : fuoHitMap) {
        std::vector<BmnFHCalPoint*>& points = pair.second;
        std::sort(points.begin(), points.end(),
                  [](BmnFHCalPoint* a, BmnFHCalPoint* b) { return a->GetTime() < b->GetTime(); });
    }
 
    LOG(debug2) << "BmnFHCalDigitizer::FillHitMap() event " << gMC->CurrentEvent() << " fired sections "
                << fuoHitMap.size();
}
 
void BmnFHCalDigitizer::Finish()
{
    printf("Work time of the FHCal digitizer: %.4f sec.\n", fworkTime);
}