BmnHodoReconstructor.cxx

Go to the documentation of this file.

/** @file BmnHodoReconstructor.cxx
 ** @author Nikolay Karpushkin <karpushkin@inr.ru>
 ** @date 03.04.2021
 **
 ** Code for BmnHodoEvent reconstruction (creation) from BmnHodoDigi {Data} or BmnHodoDigit {Sim}
 **/
 
#include "BmnHodoReconstructor.h"
 
#include "BmnHodoGeo.h"
 
BmnHodoReconstructor::BmnHodoReconstructor(bool isExp)
    : fArrayOfDigits(nullptr)
    , fBmnHodoEvent(nullptr)
    , fCuts(nullptr)
    , fIsExp(isExp)
    , fworkTime(0.)
{
    LOG(detail) << Form("BmnHodoReconstructor: isExp %i", fIsExp);
}
 
void BmnHodoReconstructor::SetRecoCutsFile(const std::string& cuts)
{
    std::string fullPath = std::string(getenv("VMCWORKDIR")) + "/parameters/hodo/" + cuts;
    std::ifstream test(fullPath);
    if (!test.is_open()) {
        LOG(warn) << Form("BmnHodoReconstructor: cuts file not found: %s. Using default (no cuts).", fullPath.c_str());
        fCuts = nullptr;
        return;
    }
 
    fCuts = new BmnHodoCuts(BmnHodoCuts::Load(fullPath, fIsExp));
    LOG(info) << Form("BmnHodoReconstructor: cuts loaded from %s (isExp=%d)", fullPath.c_str(), fIsExp);
}
 
void BmnHodoReconstructor::ParseGeometry()
{
    BmnHodoGeo* GeoHandler = new BmnHodoGeo();
    GeoHandler->ReadGeometryFromGeoManager(true);
    const auto& positionMap = GeoHandler->GetPositionMap();
    for (auto it : positionMap) {
        auto unique_address = it.first;
        if (unique_address == 0)
            continue;
        auto strip_id = BmnHodoAddress::GetStripId(unique_address);
        auto* strip = fBmnHodoEvent->GetOrAddStrip(strip_id);
        strip->SetStripId(strip_id);
        strip->SetPosition(it.second.first);
        strip->SetPositionError(it.second.second);
    }
    delete GeoHandler;
}
 
BmnHodoReconstructor::~BmnHodoReconstructor()
{
    if (fBmnHodoEvent)
        delete fBmnHodoEvent;
    if (fCuts)
        delete fCuts;
}
 
InitStatus BmnHodoReconstructor::Init()
{
    fworkTime = 0.;
    fpFairRootMgr = FairRootManager::Instance();
    if (!fpFairRootMgr) {
        LOG(error) << "BmnHodoReconstructor::Init(): FairRootManager instance is nullptr! Task will be deactivated";
        SetActive(kFALSE);
        return kERROR;
    }
 
    fArrayOfDigits = (TClonesArray*)fpFairRootMgr->GetObject((fIsExp) ? "HodoDigi" : "HodoDigit");
    if (!fArrayOfDigits) {
        LOG(error) << "BmnHodoReconstructor::Init(): branch with Digits not found! Task will be deactivated";
        SetActive(kFALSE);
        return kERROR;
    }
 
    fBmnHodoEvent = new BmnHodoEvent();
    fBmnHodoEvent->reset();
    ParseGeometry();
    fpFairRootMgr->RegisterAny("HodoEvent", fBmnHodoEvent, kTRUE);
    fpFairRootMgr->Register("HodoEvent", "Hodo", fBmnHodoEvent, kFALSE);
 
    LOG(debug) << "Hodo Reconstructor ready";
    return kSUCCESS;
}
 
void BmnHodoReconstructor::ProcessStrip(BmnHodoStrip& strip)
{
    BmnHodoDigit* downDigit = nullptr;
    BmnHodoDigit* upDigit = nullptr;
 
    int n = 0;
    double logSum = 0.0;
    for (auto* digit : strip.GetDigits()) {
        if (!digit)
            continue;
 
        const double s = digit->GetSignal();
        if (s < 0.0)
            continue;
 
        logSum += std::log(s);
        ++n;
 
        if (digit->GetStripSide() == 0 && !downDigit)
            downDigit = digit;
        else if (digit->GetStripSide() == 1 && !upDigit)
            upDigit = digit;
    }
 
    if (n == 0) {
        strip.SetSignal(0.f);
        strip.SetIntegral(0.f);
        return;
    }
 
    const double signal = std::exp(logSum / n);
    strip.SetSignal(signal);
 
    if (fIsExp) {
        double integral = 0.0, dn = 0.0, up = 0.0;
 
        if (downDigit)
            dn = static_cast<BmnHodoDigi*>(downDigit)->fSignalIntegral;
        if (upDigit)
            up = static_cast<BmnHodoDigi*>(upDigit)->fSignalIntegral;
 
        if (dn > 0.0 && up > 0.0)
            integral = std::sqrt(dn * up);
        else if (dn > 0.0)
            integral = dn;
        else if (up > 0.0)
            integral = up;
 
        strip.SetIntegral(integral);
    } else {
        strip.SetIntegral(signal);
    }
}
 
void BmnHodoReconstructor::Exec(Option_t* opt)
{
    if (!IsActive())
        return;
 
    TStopwatch sw;
    sw.Start();
 
    fBmnHodoEvent->ResetStrips();
    BuildFilteredStripMap();
 
    auto* strips = fBmnHodoEvent->GetStrips();
    if (!strips)
        return;
 
    for (int i = 0, n = strips->GetEntriesFast(); i < n; ++i) {
        auto* strip = static_cast<BmnHodoStrip*>(strips->At(i));
        if (!strip)
            continue;
 
        ProcessStrip(*strip);
    }
 
    fBmnHodoEvent->SummarizeEvent();
 
    sw.Stop();
    fworkTime += sw.RealTime();
}
 
std::vector<BmnHodoDigit*> BmnHodoReconstructor::GetDenoisedDigis()
{
    std::vector<BmnHodoDigit*> denoised;
    denoised.reserve(fArrayOfDigits->GetEntriesFast());
 
    for (int i = 0, n = fArrayOfDigits->GetEntriesFast(); i < n; ++i) {
        if (fIsExp) {
            auto* d = static_cast<BmnHodoDigi*>(fArrayOfDigits->At(i));
            if (!d)
                continue;
            if (d->GetGain() == 1)
                continue;
            if (!fCuts || fCuts->PassExp(d))
                denoised.push_back(d);
        } else {
            auto* d = static_cast<BmnHodoDigit*>(fArrayOfDigits->At(i));
            if (!d)
                continue;
            if (!fCuts || fCuts->PassSim(d))
                denoised.push_back(d);
        }
    }
 
    LOG(debug4) << Form("Hodo denoise: kept %zu / %d (%.1f%%)%s", denoised.size(), fArrayOfDigits->GetEntriesFast(),
                        100.0 * denoised.size() / std::max(1, fArrayOfDigits->GetEntriesFast()),
                        fCuts ? "" : " [no cuts applied]");
    return denoised;
}
 
void BmnHodoReconstructor::MarkPileup(const std::vector<BmnHodoDigit*>& input)
{
    std::unordered_map<uint32_t, std::vector<BmnHodoDigi*>> byAddr;
    byAddr.reserve(input.size());
 
    for (auto* d : input)
        byAddr[d->GetAddress()].push_back(static_cast<BmnHodoDigi*>(d));
 
    for (auto& [addr, vec] : byAddr) {
        const int strip_id = BmnHodoAddress::GetStripId(addr);
        auto* strip = fBmnHodoEvent->GetOrAddStrip(strip_id);
 
        if (vec.size() == 1) {
            auto* d = vec.front();
            if (d->fNpeaks > 1)
                strip->SetPileUpStatus(BmnHodoStrip::PileUp::Resolved);
        } else {
            strip->SetPileUpStatus(BmnHodoStrip::PileUp::Unresolved);
        }
    }
    return;
}
 
// Build map: strip_id -> vector of digis
void BmnHodoReconstructor::BuildFilteredStripMap()
{
    // 1. denoise
    auto denoised = GetDenoisedDigis();
 
    // 2. exp only: mark pileup
    if (fIsExp)
        MarkPileup(denoised);
 
    // 3. build map strip_id -> digis
    for (auto* digi : denoised)
        fBmnHodoEvent->GetOrAddStrip(digi->GetStripId())->AddDigit(digi);
}
 
void BmnHodoReconstructor::OnlineWrite(const std::unique_ptr<TTree>& resultTree)
{
    if (!IsActive())
        return;
 
    resultTree->Branch("HodoEvent", &fBmnHodoEvent);
    resultTree->Fill();
}
 
void BmnHodoReconstructor::Finish()
{
    printf("Work time of BmnHodoReconstructor: %4.2f sec.\n", fworkTime);
}