BmnAligner.h

Go to the documentation of this file.

/*
    BM@N alignment routine
    BM@N experiment at NICA complex, JINR, 2025
 
    Department: Math & Soft Group of HEP lab
    Author:     Igor Polev, polev@jinr.ru
 
    BmnAligner class declaration
    Alignment algorithm itself
*/
 
#ifndef BMNALIGNER_H
#define BMNALIGNER_H
 
#include "BmnAlignDefines.h"
#include "BmnAlignResult.h"
#include "BmnDataIterator.h"
#include "BmnDetectorModel.h"
#include "BmnHitRecord.h"
#include "BmnMeasureModel.h"
#include "BmnRamIterator.h"
#include "BmnRootIterator.h"
#include "BmnSimpleProgressBar.h"
#include "Math/BrentMinimizer1D.h"
#include "TPad.h"
#include "TStopwatch.h"
 
#include <array>
#include <deque>
#include <mutex>
#include <thread>   // IWYU pragma: export
#include <vector>
 
// different types of hits to be used in test tasks and real tasks
template<typename HitType = BmnHitRecord>
class BmnAligner
{
  public:
    BmnAligner(   // basic constructor
        const BmnMeasureModel* measureModel,
        const BmnDetectorModel* modelDetector,
        const char* constraintsPath = nullptr,
        Int_t nThreads = 1);
    BmnAligner(   // main constructor
        const BmnMeasureModel* measureModel,
        const BmnDetectorModel* modelDetector,
        const char* filesListPath,
        const char* constraintsPath = nullptr,
        Int_t nThreads = 1);
    virtual ~BmnAligner();
 
    // user interface methods
    Bool_t Align(Bool_t fixEvery = kFALSE);
    virtual Bool_t PrepareData(const char* tracksPath = nullptr);
    Bool_t SaveSolution(const char* jsonPath, Bool_t withAlpha = kFALSE);
    Bool_t LoadSolution(const char* jsonPath, Bool_t withAlpha = kTRUE, Bool_t withMSE = kTRUE);
    virtual void Draw();
    // getters
    Bool_t Initialized() const noexcept { return fInitialized; }
    Int_t TotalTracks() const noexcept { return fTotalTracks; }
    Int_t TotalHits() const noexcept { return fTotalHits; }
    Double_t OmegaScaleFactor() const noexcept { return fOmegaScaleFactor; }
    BmnAlignResult* GetCurrentResult() const noexcept { return fpResultCurrent; }
    BmnAlignResult* GetResult(Bool_t last = kTRUE) const noexcept
    {
        if (fIterationResults.empty())
            return nullptr;
        else if (last)
            return const_cast<BmnAlignResult*>(&fIterationResults.back());
        else
            return const_cast<BmnAlignResult*>(&fIterationResults.front());
    }
 
  protected:
    // main internal methods
    virtual Bool_t IterateAlignment();
    Double_t CalculateMSE(const SVectGLT& solution);       // does not store result
    virtual Double_t CalculateMSE(Bool_t withResiduals);   // stores result in fResultCurrent
    Bool_t SolveGlobal(const TMatrixD& mCL, const TVectorD& vBL);
    virtual void ReportResults() const;
    // internal service methods
    Bool_t LoadConstraints();
    void InitIterators();
    void PrepareDrawing();
    void DrawResiduals();
    inline void AddBlankSolution();
    inline void StoreZeroSolution(const std::vector<SVectLC>& AlphaZero);
    inline void MarkTime(TStopwatch& timer);
    void DoneInTime(TStopwatch& timer, const char* stage);
 
    typedef std::array<Bool_t, BMN_MODULE_COUNT> DetectorFlags_t;
    typedef std::array<Int_t, BMN_MODULE_COUNT> DetectorHitCount_t;
    typedef ROOT::Math::BrentMinimizer1D Minimizer1D_t;
    // service variables
    BmnSimpleProgressBar fProgressBar;
    Bool_t fInitialized{kFALSE}, fMSEclcSlow{kFALSE};
    UInt_t fMSEclcCounter{0u}, fMSEclcTick{1u};
    Double_t fTotalCpuTime{0.0}, fTotalRealTime{0.0},
        // calculation state and parameters
        fOmegaScaleFactor{1.0};
    Int_t fIteration{-1}, fRegularIters{0}, fActiveDetCount{BMN_MODULE_COUNT};
    DetectorFlags_t fDetectorActive;
    std::vector<Bool_t>* fpTrackUsable{nullptr};   // TODO: replace with TBits
    Int_t fTracksUsableCnt{0}, fTotalTracks{0}, fTotalHits{0};
    DetectorHitCount_t fHitsPerDetector{};
    Minimizer1D_t fBrent;
    // task definition
    const BmnMeasureModel* fpMeasureModel;
    const BmnDetectorModel* fpDetModel;
    const char* fpConstraintsPath;
    std::vector<TVectorD> fConstraints{};   // array of constraint vectors c
    Int_t fConstraintsCnt{0};
    // data storage implementation
    BmnRootIterator* fpRootIterator{nullptr};          // iterator over ROOT files
    BmnRamIterator<HitType>* fpRamIterator{nullptr};   // iterator over RAM storage
    BmnDataIterator* fpIteratorMain{nullptr};          // currently used iterator
    HitType* fpHits{nullptr};                          // typed link to raw hits data storage
    Byte_t* fpHitsData{nullptr};                       // raw data storage for hits
    Int_t* fpFirstHits{nullptr};                       // indexes of first hit of each track
    // iteration results storage
    BmnAlignResult *fpResultLast{nullptr},          // link to result of previous iteration
        *fpResultCurrent{nullptr};                  // link to result of latest iteration
    std::deque<BmnAlignResult> fIterationResults;   // all results of calculations
 
    // drawing stuff - works only in ROOT macro mode
    // TODO: move drawing functionality to some other class (? BmnAlignResult)
    TPad *fpPad3D{nullptr}, *fpPadG1{nullptr}, *fpPadG2{nullptr}, *fpPadG3{nullptr}, *fpPadG4{nullptr};
    Double_t fMaxZvalue{0.0};
 
    // multithreading implementation
  private:
    std::vector<BmnDataIterator*> fpIterators;   // one iterator per thread
    std::vector<Int_t> fFirstTracksIdx;          // first track to process for each thread
    std::mutex fMutexGlobal;
    Int_t fThreadsCnt{1};
    // thread functions: must have last Int_t argument to support _BMN_ALIGNER_RUN_THREADS macro
    void ProcessTracks(TVectorD* pvB_dash, TMatrixD* pmC_dash, Int_t threadNum);
    void CalculateMSEpart(Bool_t withResiduals, Int_t threadNum);
 
  protected:
    Int_t GetThreadCount() const noexcept { return fThreadsCnt; }
    BmnSimpleProgressBar& GetProgressBar() noexcept { return fProgressBar; }
};
 
// macro is used instead of method for simplicity of variable function arguments implementation
// GetProgressBar() and GetThreadCount() methods are used instead of direct fields access to make
// the macro available for descendant classes while keeping multithread implementation private
#define _BMN_ALIGNER_RUN_THREADS(...)                              \
    {                                                              \
        std::thread** _THRDS = new std::thread*[GetThreadCount()]; \
        for (auto i = 0; i < GetThreadCount(); i++)                \
            _THRDS[i] = new std::thread(__VA_ARGS__, i);           \
        for (auto i = 0; i < GetThreadCount(); i++) {              \
            _THRDS[i]->join();                                     \
            delete _THRDS[i];                                      \
        }                                                          \
        delete[] _THRDS;                                           \
    }
 
// inline implementation
 
template<typename HitType>
inline BmnAligner<HitType>::BmnAligner(const BmnMeasureModel* measureModel,
                                       const BmnDetectorModel* modelDetector,
                                       const char* constraintsPath,
                                       Int_t nThreads)
    : fpMeasureModel(measureModel)
    , fpDetModel(modelDetector)
    , fpConstraintsPath(constraintsPath)
    , fThreadsCnt(nThreads)
{
    TH1D::AddDirectory(kFALSE);           // unbind all historgrams from default file
    TH1D::StatOverflows(kTRUE);           // do not ignore overflows in mean and stddev
    fBrent.SetNpx(BMN_BRENT_GRID_SIZE);   // lower bracketing points number to speed up Brent algorithm
    fDetectorActive.fill(kTRUE);          // parameters of all detectors are active initially
};
 
template<typename HitType>
inline void BmnAligner<HitType>::AddBlankSolution()
{
    fIterationResults.emplace_back(fTotalTracks);
    fpResultLast = fpResultCurrent;
    fpResultCurrent = &fIterationResults.back();
}
 
template<typename HitType>
inline void BmnAligner<HitType>::StoreZeroSolution(const std::vector<SVectLC>& AlphaZero)
{
    AddBlankSolution();
    fpResultCurrent->Alpha() = AlphaZero;
    CalculateMSE(kTRUE);
}
 
template<typename HitType>
inline void BmnAligner<HitType>::MarkTime(TStopwatch& timer)
{
    timer.Stop();
    fTotalCpuTime += timer.CpuTime();
    fTotalRealTime += timer.RealTime();
}
 
#endif   // BMNALIGNER_H