BmnScWallDigitizer.cxx

Go to the documentation of this file.

/*
 * File:   BmnScWallDigitizer.cxx
 * Author: Sergey Morozov
 *
 * Created on 16.09.2021, 12:00
 */
 
#include "BmnScWallDigitizer.h"
 
#include "BmnScWallAddress.h"
#include "FairRootManager.h"
#include "TStopwatch.h"
 
#include <FairRun.h>
#include <FairRunSim.h>
#include <nlohmann/json.hpp>
 
BmnScWallDigitizer::BmnScWallDigitizer(int period)
    : BmnScWallDigitizer(Form("SCWALL_digitizer_period%d.json", period))
{}
 
BmnScWallDigitizer::BmnScWallDigitizer(const char* config)
    : FairTask("BmnScWallDigitizer")
    , fScale(1.0)
    , fThreshold(0.0)
    , fTimeCut(1000.0)
    , fSiPM(2668, 55, 0.1, 0.0021)
    , fPointArray(nullptr)
    , fDigiArray(nullptr)
    , fworkTime(0.0f)
{
    if (config && config[0] != '\0') {
        LoadConfig(config);
    }
}
 
InitStatus BmnScWallDigitizer::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/scwall/" + 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
BmnScWallDigitizer::~BmnScWallDigitizer()
{
    if (fDigiArray) {
        fDigiArray->Delete();
        delete fDigiArray;
    }
}
 
// Initialization method
InitStatus BmnScWallDigitizer::Init()
{
    fworkTime = 0.0;
    LOG(detail) << "-I- BmnScWallDigitizer: Init started..." << std::endl;
 
    FairRootManager* ioman = FairRootManager::Instance();
    if (!ioman) {
        LOG(error) << "-E- BmnScWallDigitizer::Init: RootManager not instantiated!" << std::endl;
        return kFATAL;
    }
 
    fPointArray = static_cast<TClonesArray*>(ioman->GetObject("ScWallPoint"));
    if (!fPointArray) {
        LOG(error) << "-E- BmnScWallDigitizer::Init: No ScWallPoint array!" << std::endl;
        return kERROR;
    }
 
    fDigiArray = new TClonesArray("BmnScWallDigit");
    ioman->Register("ScWallDigit", "ScWall", fDigiArray, kTRUE);
 
    LOG(detail) << "-I- BmnScWallDigitizer: Initialization successful" << std::endl;
    return kSUCCESS;
}
 
void BmnScWallDigitizer::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
    TStopwatch sw;
    sw.Start();
 
    LOG(debug2) << "BmnScWallDigitizer::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;
 
            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()) {
            continue;
        } else {
            TClonesArray& ar = *fDigiArray;
            long entries = fDigiArray->GetEntriesFast();
            new (ar[entries]) BmnScWallDigit(address, crossTime, totalEnergy);
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug2))
                static_cast<BmnScWallDigit*>(ar.At(entries))->Print();
        }
    }
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
void BmnScWallDigitizer::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<BmnScWallPoint*>(fPointArray->At(iPoint));
        uint32_t address = BmnScWallAddress::GetAddress(point->GetCopyMother());
        fuoHitMap[address].push_back(point);
    }
 
    // Sort each vector of points by their time
    for (auto& pair : fuoHitMap) {
        std::vector<BmnScWallPoint*>& points = pair.second;
        std::sort(points.begin(), points.end(),
                  [](BmnScWallPoint* a, BmnScWallPoint* b) { return a->GetTime() < b->GetTime(); });
    }
 
    LOG(debug2) << "BmnScWallDigitizer::FillHitMap() event " << gMC->CurrentEvent() << " fired sections "
                << fuoHitMap.size();
}
 
void BmnScWallDigitizer::Finish()
{
    printf("Work time of the ScWall digitizer: %.4f sec.\n", fworkTime);
}