BmnNdetDigitizer.cxx

Go to the documentation of this file.

/*************************************************************************************
 *            BmnNdetDigitizer
 *    Class to create digital data taken from BmnNdet detector
 *  Author:   Elena Litvinenko
 *  e-mail:   litvin@nf.jinr.ru
 *  Version:  18-11-2015
 *  Modified by M.Golubeva July 2022
 ************************************************************************************/
 
#include "BmnNdetDigitizer.h"
 
#include "BmnNdetCell.h"
#include "FairRootManager.h"
#include "FairRun.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "TStopwatch.h"
#include "TVirtualMC.h"
 
#include <nlohmann/json.hpp>
 
BmnNdetDigitizer::BmnNdetDigitizer(int period)
    : BmnNdetDigitizer(Form("NDET_digitizer_period%d.json", period))
{}
 
BmnNdetDigitizer::BmnNdetDigitizer(const char* config)
    : FairTask("BmnNdetDigitizer")
    , fScale(1.0)
    , fThreshold(0.0025)
    , fTimeCut(100.0)
    , fSiPM(160000, 165, 0.003, 0.005)
    , fPointArray(nullptr)
    , fDigiArray(nullptr)
    , fworkTime(0.0f)
{
    if (config && config[0] != '\0') {
        LoadConfig(config);
    }
}
 
InitStatus BmnNdetDigitizer::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/ndet/" + 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;
}
 
BmnNdetDigitizer::~BmnNdetDigitizer()
{
    if (fDigiArray) {
        fDigiArray->Delete();
        delete fDigiArray;
    }
}
 
InitStatus BmnNdetDigitizer::Init()
{
    fworkTime = 0.0;
    LOG(detail) << "-I- BmnNdetDigitizer: Init started..." << std::endl;
 
    // Get RootManager
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << "-E- BmnNdetDigitizer::Init: RootManager not instantiated!" << std::endl;
        return kFATAL;
    }
 
    // Get input array
    fPointArray = static_cast<TClonesArray*>(ioman->GetObject("NdetPoint"));
    if (!fPointArray) {
        LOG(error) << "-W- BmnNdetDigitizer::Init: No NdetPoint array!" << std::endl;
        return kERROR;
    }
 
    // Create and register output array
    fDigiArray = new TClonesArray("BmnNdetDigit");
    ioman->Register("NdetDigit", "Ndet", fDigiArray, kTRUE);
 
    LOG(detail) << "-I- BmnNdetDigitizer: Initialization successful" << std::endl;
    return kSUCCESS;
}
 
void BmnNdetDigitizer::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
    TStopwatch sw;
    sw.Start();
 
    LOG(debug) << "BmnNdetDigitizer::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 (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;
                point->AddLink(FairLink(static_cast<int>(BmnNdetCell::LinkType::SurfaceMCPointCrossThr),
                                        point->GetLinksWithType(static_cast<int>(BmnNdetCell::LinkType::SurfaceMCPoint))
                                            .GetLink(0)
                                            .GetIndex(),
                                        point->GetEnergyLoss()));
                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) << "BmnNdetDigitizer::Skip Digit " << BmnNdetAddress::GetInfoString(address).c_str();
            continue;
        } else {
            TClonesArray& ar = *fDigiArray;
            long entries = fDigiArray->GetEntriesFast();
            new (ar[entries]) BmnNdetDigit(address, crossTime, totalEnergy);
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug2))
                static_cast<BmnNdetDigit*>(ar.At(entries))->Print();
        }
    }
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
void BmnNdetDigitizer::FillHitMap()
{
    fuoHitMap.clear();
 
    // Fill the hit map with points from the TClonesArray
    for (int iPoint = 0, nPoints = fPointArray->GetEntriesFast(); iPoint < nPoints; ++iPoint) {
        auto* point = static_cast<BmnNdetPoint*>(fPointArray->At(iPoint));
        uint32_t address = point->GetAddress();
        fuoHitMap[address].push_back(point);
    }
 
    // Sort each vector of points by their time
    for (auto& pair : fuoHitMap) {
        std::vector<BmnNdetPoint*>& points = pair.second;
        std::sort(points.begin(), points.end(),
                  [](BmnNdetPoint* a, BmnNdetPoint* b) { return a->GetTime() < b->GetTime(); });
    }
 
    LOG(debug) << "BmnNdetDigitizer::FillHitMap() event " << gMC->CurrentEvent() << " fired cells " << fuoHitMap.size();
}
 
void BmnNdetDigitizer::Finish()
{
    printf("Work time of the Ndet digitizer: %.4f sec.\n", fworkTime);
}