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.int
 
    BmnAlignResult class header-only implementation
    Stores results of alignment iteration
*/
 
#ifndef BMNALIGNRESULT_H
#define BMNALIGNRESULT_H
 
#include "BmnAlignDefines.h"
 
#include <TH1D.h>
#include <TH2D.h>
#include <TMatrixD.h>
#include <TMatrixDSym.h>
#include <TVectorD.h>
#include <algorithm>
#include <atomic>
#include <stdexcept>
 
/*
 *  All per-parameter quantities below live in the *full* parameter layout the
 *  detector model's hierarchy defines: every element followed by every
 *  composite node, six corrections each (see BmnParamLayout). A leaf entry
 *  holds the *effective* correction - the element's own share plus everything
 *  its ancestors contribute - while a composite entry holds that node's own
 *  share; SolveGlobal() maintains the identity between them. Keeping the two
 *  levels in one vector is what lets an error or a covariance be quoted for a
 *  composite node at all.
 */
class BmnAlignResult
{
  private:
    // Binning of the histograms this class owns. Ranges are in cm and chosen for
    // the residual scale the alignment works at, so they belong to these
    // histograms and to nothing else.
    static constexpr Int_t RESIDUALS_BINS{201};
    static constexpr Double_t RESIDUALS_RANGE_X{0.2};
    static constexpr Double_t RESIDUALS_RANGE_Y{0.4};
    static constexpr Double_t TARGET_Z{0.0};            // measured Z position of the target
    static constexpr Double_t TARGET_THICKNESS{1.75};   // thickness of the target (along Z axis)
 
    VecSVectLC_t fAlpha;   // local parameters for each track
    TVectorD fParams,      // corrections, full layout (see above)
        fErrors,           // their 1-sigma errors (diagonal of V = C^-1); populated only
                           // by the non-regularized SolveGlobal call (mirrors SVDSig/SVDv)
        fGlCorr;           // global correlation coefficients, same layout
    VecSVectGL_t fA;       // per-element copy of the leaf block, for the per-hit
                           // code that wants one element's six corrections together
    TVectorD fL,           // Lagrange multipliers for global parameters constrains
        fSVDSig;           // singular values vector from SVD decomposition of MNE
    TMatrixD fSVDv;        // right singular vectors (columns) of the non-regularized MNE
    TMatrixDSym fCov;      // covariance in the full layout; 0x0 if not available.
                           // Stored only on the latest iteration result.
    TH1D fResidualsX,      // histogram of residuals by X coordinate
        fResidualsY,       // histogram of residuals by Y coordinate
        fVertexZ;          // histogram of primary vertex distribution by Z coordinate
    TH2D fResidualsSX,     // "syntetic" histogram of X residuals for visualization
        fResidualsSY;      // "syntetic" histogram of Y residuals for visualization
    Double_t fMSEtotal;    // total value of target function at current solution
    Int_t fAlphaCount;     // number of tracks (set in constructor)
    Int_t fElements;       // alignable detector elements, i.e. leaf nodes
 
  public:
    BmnAlignResult(Int_t TrackCount, Double_t MaxZvalue, Int_t Elements, Int_t AllParams);
    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) {};
 
    // non-const accessors
 
    VecSVectLC_t& Alpha() noexcept { return fAlpha; };
    SVectLC& Alpha(Int_t index) noexcept { return fAlpha[index]; };
    SVectGL& A(Int_t index) noexcept { return fA[index]; };
    TVectorD& A() noexcept { return fParams; };
    TVectorD& AErr() noexcept { return fErrors; };
    Double_t& AErr(Int_t index) noexcept { return fErrors(index); };
    TVectorD& AglCorr() noexcept { return fGlCorr; };
    Double_t& AglCorr(Int_t index) noexcept { return fGlCorr(index); };
    TVectorD& L() noexcept { return fL; };
    TVectorD& SVDSig() noexcept { return fSVDSig; };
    TMatrixD& SVDv() noexcept { return fSVDv; };
    TMatrixDSym& Cov() noexcept { return fCov; };
 
    TH1D& ResidualsX() noexcept { return fResidualsX; };
    TH1D& ResidualsY() noexcept { return fResidualsY; };
    TH1D& VertexZ() noexcept { return fVertexZ; };
    TH2D& ResidualsSX() noexcept { return fResidualsSX; };
    TH2D& ResidualsSY() noexcept { return fResidualsSY; };
 
    // const getters
 
    const VecSVectLC_t& Alpha() const noexcept { return fAlpha; };
    const SVectLC& Alpha(Int_t index) const noexcept { return fAlpha[index]; };
    const SVectGL& A(Int_t index) const noexcept { return fA[index]; };
    const TVectorD& A() const noexcept { return fParams; };
    const TVectorD& AErr() const noexcept { return fErrors; };
    const Double_t& AErr(Int_t index) const noexcept { return fErrors(index); };
    const TVectorD& AglCorr() const noexcept { return fGlCorr; };
    const Double_t& AglCorr(Int_t index) const noexcept { return fGlCorr(index); };
    const TVectorD& L() const noexcept { return fL; };
    const TVectorD& SVDSig() const noexcept { return fSVDSig; };
    const TMatrixD& SVDv() const noexcept { return fSVDv; };
    const TMatrixDSym& Cov() const noexcept { return fCov; };
 
    const TH1D& ResidualsX() const noexcept { return fResidualsX; };
    const TH1D& ResidualsY() const noexcept { return fResidualsY; };
    const TH1D& VertexZ() const noexcept { return fVertexZ; };
    const TH2D& ResidualsSX() const noexcept { return fResidualsSX; };
    const TH2D& ResidualsSY() const noexcept { return fResidualsSY; };
 
    Double_t GetValueMSE() const noexcept { return fMSEtotal; };
    Int_t GetAlphaCount() const noexcept { return fAlphaCount; };
    ULong_t GetHistCounter() const noexcept;
 
    void AddValueMSE(Double_t val) noexcept { fMSEtotal += val; };
    void ResetValueMSE() noexcept { fMSEtotal = 0.0; };
 
    void SeparateA();
    void ConcatenateA();
};
 
// Inline implementations
 
inline BmnAlignResult::BmnAlignResult(Int_t TrackCount, Double_t MaxZvalue, Int_t Elements, Int_t AllParams)
    : fAlpha(TrackCount)
    , fParams(AllParams)
    , fErrors(AllParams)
    , fGlCorr(AllParams)
    , fA(Elements)
    , fL()
    , fSVDSig()
    , fSVDv()
    , fCov(AllParams)
    , fResidualsX(Form("residualsX_%lu", GetHistCounter()),
                  "Total residuals by X (red: aligned)",
                  RESIDUALS_BINS,
                  -RESIDUALS_RANGE_X,
                  RESIDUALS_RANGE_X)
    , fResidualsY(Form("residualsY_%lu", GetHistCounter()),
                  "Total residuals by Y (red: aligned)",
                  RESIDUALS_BINS,
                  -RESIDUALS_RANGE_Y,
                  RESIDUALS_RANGE_Y)
    , fVertexZ(Form("vertexZ_%lu", GetHistCounter()),
               "Primary vertex distribution;Z, cm;tracks",
               RESIDUALS_BINS,
               TARGET_Z - TARGET_THICKNESS * 2.0,
               TARGET_Z + TARGET_THICKNESS * 2.0)
    , fResidualsSX(Form("residualsSX_%lu", GetHistCounter()),
                   "Review of X residuals along beam axis in FSD region;Z, cm;X residual, cm",
                   RESIDUALS_BINS * 10,
                   14.0,
                   49.0,   // Z range to cover all FSD stations
                   RESIDUALS_BINS,
                   -RESIDUALS_RANGE_X,
                   RESIDUALS_RANGE_X)
    , fResidualsSY(Form("residualsSY_%lu", GetHistCounter()),
                   "Review of Y residuals along beam axis in FSD region;Z, cm;Y residual, cm",
                   RESIDUALS_BINS * 10,
                   14.0,
                   49.0,   // Z range to cover all FSD stations
                   RESIDUALS_BINS,
                   -RESIDUALS_RANGE_Y,
                   RESIDUALS_RANGE_Y)
    , fMSEtotal(0.0)
    , fAlphaCount(TrackCount)
    , fElements(Elements)
{}
 
inline BmnAlignResult::BmnAlignResult(const BmnAlignResult& other)
    : fAlpha(other.fAlpha)
    , fParams(other.fParams)
    , fErrors(other.fErrors)
    , fGlCorr(other.fGlCorr)
    , fA(other.fA)
    , fL(other.fL)
    , fSVDSig(other.fSVDSig)
    , fSVDv(other.fSVDv)
    , fCov(other.fCov)
    , fResidualsX(other.fResidualsX)
    , fResidualsY(other.fResidualsY)
    , fVertexZ(other.fVertexZ)
    , fResidualsSX(other.fResidualsSX)
    , fResidualsSY(other.fResidualsSY)
    , fMSEtotal(other.fMSEtotal)
    , fAlphaCount(other.fAlphaCount)
    , fElements(other.fElements)
{
    // Update histogram names to be unique using common counter
    fResidualsX.SetName(Form("residualsX_%lu", GetHistCounter()));
    fResidualsY.SetName(Form("residualsY_%lu", GetHistCounter()));
    fVertexZ.SetName(Form("vertexZ_%lu", GetHistCounter()));
    fResidualsSX.SetName(Form("residualsSX_%lu", GetHistCounter()));
    fResidualsSY.SetName(Form("residualsSY_%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);
    fParams.ResizeTo(other.fParams);
    fErrors.ResizeTo(other.fErrors);
    fGlCorr.ResizeTo(other.fGlCorr);
    fSVDSig.ResizeTo(other.fSVDSig);
    fSVDv.ResizeTo(other.fSVDv);
    fCov.ResizeTo(other.fCov);
    fAlpha = other.fAlpha;
    fParams = other.fParams;
    fErrors = other.fErrors;
    fGlCorr = other.fGlCorr;
    fA = other.fA;
    fL = other.fL;
    fSVDSig = other.fSVDSig;
    fSVDv = other.fSVDv;
    fCov = other.fCov;
    fResidualsX = other.fResidualsX;
    fResidualsY = other.fResidualsY;
    fVertexZ = other.fVertexZ;
    fResidualsSX = other.fResidualsSX;
    fResidualsSY = other.fResidualsSY;
    fMSEtotal = other.fMSEtotal;
    fElements = other.fElements;
 
    return *this;
}
 
// fA duplicates the leaf block of fParams so that the per-hit code can take a
// whole element at once. The two directions below are the sync points; the leaf
// block is the authoritative copy.
inline void BmnAlignResult::SeparateA()
{
    const Double_t* src{fParams.GetMatrixArray()};
    for (Int_t i = 0; i < fElements; i++, src += BMN_GLOBAL_PARAMS_PD)
        fA[i].SetElements(src, BMN_GLOBAL_PARAMS_PD);
}
 
inline void BmnAlignResult::ConcatenateA()
{
    Double_t* dst{fParams.GetMatrixArray()};
    for (Int_t i = 0; i < fElements; i++, dst += BMN_GLOBAL_PARAMS_PD)
        std::copy(fA[i].begin(), fA[i].end(), dst);
}
 
inline ULong_t BmnAlignResult::GetHistCounter() const noexcept
{
    // Common counter for all histograms in all BmnAlignResult instances.
    // Atomic because CalculateMSEpart() calls this concurrently from worker threads.
    static std::atomic<ULong_t> sHistCounter{0};
    return ++sHistCounter;
};
 
#endif   // BMNALIGNRESULT_H