BmnNdetReconstructor.cxx

Go to the documentation of this file.

/** @file BmnNdetReconstructor.cxx
 ** @author Nikolay Karpushkin <karpushkin@inr.ru>
 ** @date 03.04.2021
 **
 ** Code for BmnNdetEvent reconstruction (creation) from BmnNdetDigi {Data} or BmnNdetDigit {Sim}
 **/
 
#include "BmnNdetReconstructor.h"
 
#include "BmnNdetDigi.h"
#include "BmnNdetDigit.h"
#include "BmnNdetGeo.h"
#include "BmnNdetPoint.h"
#include "BmnNdetSurfacePoint.h"
#include "FairLogger.h"
#include "TTreeFormula.h"
 
#include <TStopwatch.h>
 
BmnNdetReconstructor::BmnNdetReconstructor(bool isExp, bool isGlobalCoordinates)
    : fBmnNdetEvent(nullptr)
    , fCuts(nullptr)
    , fIsExp(isExp)
    , fIsGlobal(isGlobalCoordinates)
    , fworkTime(0.)
{
    LOG(detail) << Form("BmnNdetReconstructor: isExp %i; isGlobal %i", fIsExp, fIsGlobal);
    fDoMatchMC = (fIsExp) ? false : true;
}
 
void BmnNdetReconstructor::SetRecoCutsFile(const std::string& cuts)
{
    std::string fullPath = std::string(getenv("VMCWORKDIR")) + "/parameters/ndet/" + cuts;
    std::ifstream test(fullPath);
    if (!test.is_open()) {
        LOG(warn) << Form("BmnNdetReconstructor: cuts file not found: %s. Using default (no cuts).", fullPath.c_str());
        fCuts = nullptr;
        return;
    }
 
    fCuts = new BmnNdetCuts(BmnNdetCuts::Load(fullPath, fIsExp));
    LOG(info) << Form("BmnNdetReconstructor: cuts loaded from %s (isExp=%d)", fullPath.c_str(), fIsExp);
}
 
BmnNdetReconstructor::~BmnNdetReconstructor()
{
    if (fBmnNdetEvent)
        delete fBmnNdetEvent;
    if (fGeoHandler)
        delete fGeoHandler;
    if (fCuts)
        delete fCuts;
}
 
InitStatus BmnNdetReconstructor::Init()
{
    fworkTime = 0.;
    fpFairRootMgr = FairRootManager::Instance();
 
    if (fIsExp) {
        if (auto* a = (TClonesArray*)fpFairRootMgr->GetObject("HgndDigi"))
            fDigits.push_back(a);
        if (auto* a = (TClonesArray*)fpFairRootMgr->GetObject("NdetDigi"))
            fDigits.push_back(a);
 
        if (fDigits.empty()) {
            LOG(error) << "BmnNdetReconstructor::Init(): no experimental digi branches (HgndDigi/NdetDigi) found! "
                          "Task will be deactivated";
            SetActive(kFALSE);
            return kERROR;
        }
        LOG(info) << Form("BmnNdetReconstructor::Init(): using %zu experimental branch(es)", fDigits.size());
    } else {
        if (auto* a = (TClonesArray*)fpFairRootMgr->GetObject("NdetDigit"))
            fDigits.push_back(a);
        else {
            LOG(error) << "BmnNdetReconstructor::Init(): branch with Digits not found! Task will be deactivated";
            SetActive(kFALSE);
            return kERROR;
        }
 
        fArrayOfPoints = (TClonesArray*)fpFairRootMgr->GetObject("NdetPoint");
        if (!fArrayOfPoints) {
            LOG(error) << "BmnNdetReconstructor::Init(): branch with NdetPoints not found! Task will be deactivated";
            SetActive(kFALSE);
            return kERROR;
        }
 
        fArrayOfSurfacePoints = (TClonesArray*)fpFairRootMgr->GetObject("NdetSurfacePoint");
        if (!fArrayOfSurfacePoints) {
            LOG(error)
                << "BmnNdetReconstructor::Init(): branch with NdetSurfacePoint not found! Task will be deactivated";
            SetActive(kFALSE);
            return kERROR;
        }
        fpFairRootMgr->Register("NdetSurfacePoint", "NDET", fArrayOfSurfacePoints, kTRUE);
    }
 
    fBmnNdetEvent = new BmnNdetEvent();
    fBmnNdetEvent->reset();
    ParseGeometry();
    fpFairRootMgr->RegisterAny("NdetEvent", fBmnNdetEvent, kTRUE);
    fpFairRootMgr->Register("NdetEvent", "NDET", fBmnNdetEvent, kFALSE);
 
    LOG(debug) << "Ndet Reconstructor ready";
    return kSUCCESS;
}
 
void BmnNdetReconstructor::ParseGeometry()
{
    fGeoHandler = new BmnNdetGeo();
    fGeoHandler->ReadGeometryFromGeoManager(fIsGlobal);
}
 
void BmnNdetReconstructor::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
    TStopwatch sw;
    sw.Start();
    fBmnNdetEvent->reset();
 
    auto digiVector = GetSelectedDigis();
    for (auto* ThisDigi : digiVector) {
        auto address = ThisDigi->GetAddress();
        const auto& posMap = fGeoHandler->GetPositionMap();
        auto pos = posMap.find(address);
        if (pos == posMap.end()) {
            LOG(debug2) << Form("BmnNdetReconstructor : Exec. Digi address %d not found in geometry. Skipping",
                                address);
            if (FairLogger::GetLogger()->IsLogNeeded(fair::Severity::debug2))
                ThisDigi->Print();
            continue;
        }
        BmnNdetCell* this_cell = new BmnNdetCell();
        this_cell->SetAddress(address);
        this_cell->SetPosition(pos->second.first);
        this_cell->SetPositionError(pos->second.second);
        this_cell->SetSignal(ThisDigi->GetSignal());
        this_cell->SetTime(ThisDigi->GetTime());
        fBmnNdetEvent->AddCell(this_cell);
    }
 
    if (!fIsExp && fDoMatchMC)
        MatchMCTracks();
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
std::vector<BmnNdetDigit*> BmnNdetReconstructor::GetSelectedDigis()
{
    std::vector<BmnNdetDigit*> denoised;
 
    size_t total = 0;
    for (auto* arr : fDigits)
        total += arr->GetEntriesFast();
    denoised.reserve(total);
 
    for (auto* arr : fDigits) {
        for (int i = 0, n = arr->GetEntriesFast(); i < n; ++i) {
            TObject* obj = arr->At(i);
            if (!obj)
                continue;
 
            if (fIsExp) {
                if (auto* dexp = dynamic_cast<BmnHgndDigi*>(obj)) {
                    if (!fCuts || fCuts->PassExp(dexp))
                        denoised.push_back(dexp);
                } else if (auto* dprt = dynamic_cast<BmnNdetDigi*>(obj)) {
                    if (!fCuts || fCuts->PassPrototype(dprt))
                        denoised.push_back(dprt);
                }
            } else {
                auto* d = static_cast<BmnNdetDigit*>(obj);
                if (!fCuts || fCuts->PassSim(d))
                    denoised.push_back(d);
            }
        }
    }
 
    LOG(debug4) << Form("BmnNdetReconstructor::GetSelectedDigis: kept %zu / %zu", denoised.size(), total);
    return denoised;
}
 
void BmnNdetReconstructor::MatchMCTracks()
{
    // Enable the use of FairLinks in FairRootManager
    fpFairRootMgr->SetUseFairLinks(kTRUE);
 
    // Iterate over all cells in the event
    auto* cells = fBmnNdetEvent->GetModifiableCells();
    for (int ice = 0, nce = cells->GetEntriesFast(); ice < nce; ++ice) {
        auto* cell = static_cast<BmnNdetCell*>(cells->At(ice));
        if (!cell)
            continue;
 
        // Iterate over all points
        for (int iPt = 0, nPt = fArrayOfPoints->GetEntriesFast(); iPt < nPt; ++iPt) {
            auto* point = static_cast<BmnNdetPoint*>(fArrayOfPoints->At(iPt));
            if (!point)
                continue;
            if (point->GetAddress() != cell->GetAddress())
                continue;   // Skip if point is from another cell
            cell->AddLink(FairLink(static_cast<int>(BmnNdetCell::LinkType::MCPoint), iPt, point->GetEnergyLoss()));
            for (auto link : point->GetLinks()) {
                link.SetWeight(point->GetEnergyLoss());
                cell->AddLink(link);
            }
            LOG(debug) << Form("BmnNdetReconstructor::MatchMCTracks. Matched MC point and MC track. Updated cell with "
                               "address %d. Now it has %d links",
                               cell->GetAddress(), cell->GetNLinks());
        }
    }
 
    // Disable the use of FairLinks in FairRootManager
    fpFairRootMgr->SetUseFairLinks(kFALSE);
}
 
void BmnNdetReconstructor::OnlineWrite(const std::unique_ptr<TTree>& resultTree)
{
    if (!IsActive())
        return;
 
    resultTree->Branch("NdetEvent", &fBmnNdetEvent);
    if (!fIsExp)
        resultTree->Branch("NdetSurfacePoint", &fArrayOfSurfacePoints);
    resultTree->Fill();
}
 
void BmnNdetReconstructor::Finish()
{
    printf("Work time of BmnNdetReconstructor: %4.2f sec.\n", fworkTime);
}