BmnAlignResult.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
 
    BmnAlignResult class header-only implementation
    Stores results of alignment iteration
*/
 
#ifndef BMNALIGNRESULT_H
#define BMNALIGNRESULT_H
 
#include "BmnAlignDefines.h"
#include "TH1D.h"
#include "TVectorD.h"
 
#include <algorithm>
#include <array>
#include <stdexcept>
#include <vector>
 
class BmnAlignResult
{
  public:
    typedef std::array<std::pair<Int_t, Double_t>, BMN_MODULE_COUNT> IdxValuePair_t;
    typedef std::array<SVectGL, BMN_MODULE_COUNT> ArraySVectGL_t;
    typedef std::array<SVectLC, BMN_MODULE_COUNT> ArraySVectLC_t;
    typedef std::array<Double_t, BMN_MODULE_COUNT> ArrayDouble_t;
    typedef std::vector<SVectLC> VectorSVectLC_t;
 
  private:
    VectorSVectLC_t fAlpha;    // array of local parameters for each track
    SVectGLT fAfull;           // global parameters for all detectors in one vector
    ArraySVectGL_t fA;         // array of global parameters for each detector
    TVectorD fL,               // Lagrange multipliers for global parameters constrains
        fSVDSig;               // singular values vector from SVD decomposition
    TH1D fResidualsX,          // histogram of residuals by X coordinate
        fResidualsY;           // histogram of residuals by Y coordinate
    ArrayDouble_t fMSE;        // MSE values per detector plane
    Double_t fMSEtotal;        // total value of target function at current solution
    const Int_t fAlphaCount;   // number of tracks (set in constructor)
 
  public:
    explicit BmnAlignResult(Int_t TrackCount);
    BmnAlignResult(const BmnAlignResult& other);
    BmnAlignResult& operator=(const BmnAlignResult& other);
    ~BmnAlignResult() = default;
    // default constructor is required for ROOT serialization and not intended to be used
    BmnAlignResult()
        : fAlphaCount(0) {};
 
    VectorSVectLC_t& Alpha() noexcept { return fAlpha; };
    SVectLC& Alpha(Int_t index) noexcept { return fAlpha[index]; };
    SVectGL& A(Int_t index) noexcept { return fA[index]; };
    SVectGLT& A() noexcept { return fAfull; };
    TVectorD& L() noexcept { return fL; };
    TVectorD& SVDSig() noexcept { return fSVDSig; };
    TH1D& ResidualsX() noexcept { return fResidualsX; };
    TH1D& ResidualsY() noexcept { return fResidualsY; };
 
    Double_t GetValueMSE(Int_t index) const noexcept { return fMSE[index]; };
    Double_t GetValueMSE() const noexcept { return fMSEtotal; };
    Int_t GetAlphaCount() const noexcept { return fAlphaCount; };
    ULong_t GetHistCounter() const noexcept;
    IdxValuePair_t GetSortedMSE() const;
 
    void AddValueMSE(Int_t detectorIndex, Double_t val);
    void ResetValueMSE();
 
    void SeparateA();
    void ConcatenateA();
};
 
// Inline implementations
 
inline BmnAlignResult::BmnAlignResult(Int_t TrackCount)
    : fAlpha(TrackCount)
    , fAfull()
    , fA()
    , fL()
    , fSVDSig()
    , fResidualsX(Form("residualsX_%lu", GetHistCounter()),
                  "Residuals by X",
                  BMN_RESIDUALS_BINS,
                  -BMN_RESIDUALS_RANGE_X,
                  BMN_RESIDUALS_RANGE_X)
    , fResidualsY(Form("residualsY_%lu", GetHistCounter()),
                  "Residuals by Y",
                  BMN_RESIDUALS_BINS,
                  -BMN_RESIDUALS_RANGE_Y,
                  BMN_RESIDUALS_RANGE_Y)
    , fMSE{}
    , fMSEtotal(0.0)
    , fAlphaCount(TrackCount)
{
    ResetValueMSE();
}
 
inline BmnAlignResult::BmnAlignResult(const BmnAlignResult& other)
    : fAlpha(other.fAlpha)
    , fAfull(other.fAfull)
    , fA(other.fA)
    , fL(other.fL)
    , fSVDSig(other.fSVDSig)
    , fResidualsX(other.fResidualsX)
    , fResidualsY(other.fResidualsY)
    , fMSE(other.fMSE)
    , fMSEtotal(other.fMSEtotal)
    , fAlphaCount(other.fAlphaCount)
{
    // Update histogram names to be unique using common counter
    fResidualsX.SetName(Form("residualsX_%lu", GetHistCounter()));
    fResidualsY.SetName(Form("residualsY_%lu", GetHistCounter()));
}
 
inline BmnAlignResult& BmnAlignResult::operator=(const BmnAlignResult& other)
{
    if (this == &other)
        return *this;
    if (fAlphaCount != other.fAlphaCount)
        throw std::runtime_error("BmnAlignResult: assignment between objects with different alphaCount");
 
    fL.ResizeTo(other.fL);
    fSVDSig.ResizeTo(other.fSVDSig);
    fAlpha = other.fAlpha;
    fAfull = other.fAfull;
    fA = other.fA;
    fL = other.fL;
    fSVDSig = other.fSVDSig;
    fResidualsX = other.fResidualsX;
    fResidualsY = other.fResidualsY;
    fMSE = other.fMSE;
    fMSEtotal = other.fMSEtotal;
 
    return *this;
}
 
inline BmnAlignResult::IdxValuePair_t BmnAlignResult::GetSortedMSE() const
{
    IdxValuePair_t result;
    for (Int_t i = 0; i < BMN_MODULE_COUNT; i++) {
        result[i] = std::make_pair(i, fMSE[i]);
    }
 
    // Sort by MSE value in ascending order (best detectors first)
    std::sort(result.begin(), result.end(),
              [](const std::pair<Int_t, Double_t>& a, const std::pair<Int_t, Double_t>& b) noexcept {
                  return a.second < b.second;
              }   // ascending order (best detectors first)
    );
 
    return result;
}
 
inline void BmnAlignResult::AddValueMSE(Int_t detectorIndex, Double_t val)
{
    if (detectorIndex < 0 || detectorIndex >= BMN_MODULE_COUNT)
        throw std::runtime_error("BmnAlignResult: detectorIndex out of range");
    fMSE[detectorIndex] += val;
    fMSEtotal += val;
}
 
inline void BmnAlignResult::ResetValueMSE()
{
    fMSE.fill(0.0);
    fMSEtotal = 0.0;
}
 
inline void BmnAlignResult::SeparateA()
{
    auto itr{fAfull.begin()};
    for (Int_t i = 0; i < BMN_MODULE_COUNT; i++, itr += BMN_GLOBAL_PARAMS_PD)
        fA[i].SetElements(itr, BMN_GLOBAL_PARAMS_PD);
}
 
inline void BmnAlignResult::ConcatenateA()
{
    for (Int_t i = 0, j = 0; i < BMN_MODULE_COUNT; i++, j += BMN_GLOBAL_PARAMS_PD)
        fAfull.Place_at(fA[i], j);
}
 
inline ULong_t BmnAlignResult::GetHistCounter() const noexcept
{
    // Common counter for all histograms in all BmnAlignResult instances
    static ULong_t sHistCounter{0};
    return ++sHistCounter;
};
 
#endif   // BMNALIGNRESULT_H