BmnAlignReport.cxx

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
 
    BmnAlignReport class implementation
*/
 
#include "BmnAlignReport.h"
 
#include "BmnAlignDefines.h"
#include "BmnRamIterator.h"
 
#include <TCanvas.h>
#include <TDatime.h>
#include <TGraph.h>
#include <TH1.h>
#include <TLatex.h>
#include <TPad.h>
#include <TPolyLine3D.h>
#include <TRandom.h>
#include <TString.h>
#include <TSystem.h>
#include <algorithm>
#include <cmath>
#include <fairlogger/Logger.h>
#include <fstream>
#include <vector>
 
namespace
{
constexpr Int_t DRAW_TRACKS_CNT{20};   // number of random tracks to draw
// A ratio or a norm this small is treated as zero rather than divided by.
// Far above machine epsilon on purpose: it is the scale below which a
// solved parameter carries no information, not a rounding guard.
constexpr Double_t SOFT_ZERO{1.0E-6};
}   // anonymous namespace
 
char* BmnAlignReport::Change(Double_t before, Double_t after)
{
    auto factor{std::abs(after / before)};
    if (factor < 1.0)
        return Form("%+.1f%%", 100.0 * (after / before - 1.0));
    else
        return Form("%+.1fom", std::log10(factor));
}
 
void BmnAlignReport::Results(const BmnAligner& aligner)
{
    LOG(state) << " > BmnAlignReport : final report ===============================================";
    if (aligner.Results().empty()) {
        LOG(error) << "<!> No results to report";
        LOG(state) << " < BmnAlignReport : report canceled ============================================";
        return;
    }
 
    const BmnAlignResult &before{aligner.ResultZero()}, &after{aligner.ResultLast()};
    Int_t idxLast{after.SVDSig().GetNrows() - 1};
    LOG(info) << Form("    %.2e -> %.2e (%s) : MSE value change", before.GetValueMSE(), after.GetValueMSE(),
                      Change(before.GetValueMSE(), after.GetValueMSE()));
    LOG(info) << Form("    %.2e -> %.2e (%s) : mean error in X change (abs)", std::abs(before.ResidualsX().GetMean()),
                      std::abs(after.ResidualsX().GetMean()),
                      Change(before.ResidualsX().GetMean(), after.ResidualsX().GetMean()));
    LOG(info) << Form("    %.2e -> %.2e (%s) : mean error in Y change (abs)", std::abs(before.ResidualsY().GetMean()),
                      std::abs(after.ResidualsY().GetMean()),
                      Change(before.ResidualsY().GetMean(), after.ResidualsY().GetMean()));
    LOG(info) << Form("    %.2e -> %.2e (%s) : error stddev in X change", before.ResidualsX().GetStdDev(),
                      after.ResidualsX().GetStdDev(),
                      Change(before.ResidualsX().GetStdDev(), after.ResidualsX().GetStdDev()));
    LOG(info) << Form("    %.2e -> %.2e (%s) : error stddev in Y change", before.ResidualsY().GetStdDev(),
                      after.ResidualsY().GetStdDev(),
                      Change(before.ResidualsY().GetStdDev(), after.ResidualsY().GetStdDev()));
    if (after.SVDSig().GetNrows() > 0) {
        LOG(info) << Form("    %.2e : latest SVD minimum singular value", after.SVDSig()(idxLast));
        LOG(info) << Form("    %.2e : latest SVD condition number", after.SVDSig()(0) / after.SVDSig()(idxLast));
    }
    const TVectorD& params{after.A()};
    const BmnAlignPriors& priors{aligner.Priors()};
    const Int_t constCount{priors.ConstCount()};
    if (constCount) {
        // A() is the effective correction and already carries what the leaf
        // coefficients imply, but a NodeName coefficient addresses the
        // composite's own share - both live in A(), in one layout
        const TVectorD violation{priors.rhs - priors.constraints * params};
        Int_t goodCount{0};
        for (Int_t i = 0; i < constCount; i++) {
            if (std::abs(violation[i]) <= SOFT_ZERO) {
                goodCount++;
                continue;
            }
            LOG(warning) << Form("    %.2e : constraint row #%d violation", std::abs(violation[i]), i + 1);
        }
        LOG(info) << Form("    %d of %d : constraints satisfied within %.0e", goodCount, constCount, SOFT_ZERO);
    }
    if (aligner.Hierarchy().CompositeCount()) {
        LOG(info) << "    Composite node corrections (own share, cm and rad):";
        TString header{Form("      %3s %-14s", "lvl", "node")};
        for (Int_t j = 0; j < BMN_GLOBAL_PARAMS_PD; j++)
            header += Form(" %11s", BmnParamLayout::ParamName(j));
        LOG(info) << header;
        for (const NodeRef& node : aligner.Layout().ByNode()) {
            if (!node.IsComposite())
                continue;
            TString line{Form("      %3d %-14s", node.level, node.name)};
            for (const ParamRef& p : node.params)
                line += Form(" %+11.4e", params(p.full));
            LOG(info) << line;
        }
    }
    // TODO: make errors report optional
    /*
    LOG(info) << "    Relative errors of parameters (sigma / |value|, '--' = undef):";
    const TVectorD& aVal {after.A()};
    const TVectorD& aErr {after.AErr()};
    for (Int_t i = 0; i < Elements(); i++) {
        const Int_t ID {aligner.DetectorModel()->EncodedID(i)};
        TString     line {Form("      det %2d (id %9d, %s) :", i, ID, aligner.DetectorModel()->DecodeID(ID).c_str())};
        for (Int_t j = 0; j < BMN_GLOBAL_PARAMS_PD; j++) {
            const Int_t g {i * BMN_GLOBAL_PARAMS_PD + j};
            if (aErr[g] <= 0.0 || std::fabs(aVal[g]) <= SOFT_ZERO)
                line += Form("   %-5s    --", BmnParamLayout::ParamName(j));
            else
                line += Form("   %-5s %5.1f%%", BmnParamLayout::ParamName(j), 100.0 * aErr[g] / std::fabs(aVal[g]));
        }
        LOG(info) << line;
    }
    */
    LOG(info) << "    Execution timings:";
    Int_t cpuTime{static_cast<Int_t>(aligner.TotalCpuTime())}, realTime{static_cast<Int_t>(aligner.TotalRealTime())};
    LOG(info) << Form("      %dh %02dm %02ds : total CPU time spent", cpuTime / 3600, (cpuTime / 60) % 60,
                      cpuTime % 60);
    LOG(info) << Form("      %dh %02dm %02ds : total real time spent", realTime / 3600, (realTime / 60) % 60,
                      realTime % 60);
    LOG(info) << Form("      x%.1f : multithreading boost factor", aligner.TotalCpuTime() / aligner.TotalRealTime());
    LOG(state) << " < BmnAlignReport : final report end ===========================================";
    LOG(state) << "   ---";
}
 
void BmnAlignReport::ModesBelow(const BmnAligner& aligner, Double_t absThreshold, Bool_t dumpNew)
{
    LOG(state) << " > BmnAlignReport : weak modes analysis";
 
    // sanity checks
    if (aligner.Results().empty()) {
        LOG(error) << " <!> no iteration results available, aborting";
        return;
    }
    const BmnAlignResult& res{aligner.ResultLast()};
    const TVectorD& sig{res.SVDSig()};
    const TMatrixD& V{res.SVDv()};
    if (sig.GetNrows() == 0 || V.GetNrows() == 0) {
        LOG(error) << " <!> SVD data not found, aborting";
        return;
    }
    const Int_t nTotal{V.GetNrows()};
    // V comes from the SVD of the solved system, which is hierarchical
    const Int_t nParams{aligner.Hierarchy().AllParams()};
    if (nTotal < nParams) {
        LOG(error) << " <!> stored V matrix smaller than param count, aborting";
        return;
    }
 
    // Rows of V are full-layout parameters: leaves first, then the composites.
    // aligner.Layout() names both - indexing a detector table by
    // row / BMN_GLOBAL_PARAMS_PD would run off its end as soon as the mode
    // reaches into the composite block.
 
    // Lazy-opened sink for auto-generated new constraints (one JSON array per call,
    // only created if at least one unconstrained mode is encountered).
    std::ofstream newCnstrFile;
    TString newCnstrPath;
    Int_t newCnstrCount{0};
    Bool_t newCnstrFailed{kFALSE};
    auto openNewCnstrFile = [&]() {
        if (newCnstrFile.is_open() || newCnstrFailed)
            return;
        TDatime dt;
        const TString stamp{Form("%04d%02d%02d_%02d%02d%02d", dt.GetYear(), dt.GetMonth(), dt.GetDay(), dt.GetHour(),
                                 dt.GetMinute(), dt.GetSecond())};
        TString srcBase;
        if (aligner.ConstraintsPath()) {
            srcBase = gSystem->BaseName(aligner.ConstraintsPath());
            const Int_t dot{srcBase.Last('.')};
            if (dot >= 0)
                srcBase.Remove(dot);
        }
        newCnstrPath = "cnstr_new";
        if (!srcBase.IsNull())
            newCnstrPath += "_from-" + srcBase;
        newCnstrPath += "_" + stamp + ".json";
        newCnstrFile.open(newCnstrPath.Data(), std::ios::out);
        if (!newCnstrFile.is_open()) {
            LOG(error) << " <!> ReportModesBelow: failed to open " << newCnstrPath
                       << " for writing; further dumps suppressed";
            newCnstrFailed = kTRUE;
            return;
        }
        newCnstrFile << "[\n";
    };
 
    // The constraints have to be compared in the basis V lives in, that is as
    // the very rows IterateAlignment() appends to the solved system.
    const TMatrixD cnstrRows{aligner.Hierarchy().PromoteConstraints(aligner.Priors().constraints)};
    VecDouble_t cnstrNorms;
    for (Int_t j = 0; j < cnstrRows.GetNrows(); j++)
        cnstrNorms.push_back(std::sqrt(TVectorD(TMatrixDRow_const(cnstrRows, j)).Norm2Sqr()));
 
    // walk weak modes from the tail of the descending singular values
    LOG(info) << Form("    %.2e : eigenvalue threshold to consider", absThreshold);
    Int_t weakFound{0};
    for (Int_t k = sig.GetNrows() - 1; k >= 0; k--) {
        if (sig(k) > absThreshold)
            break;
        weakFound++;
 
        // split the mode into (param, multiplier) subspaces
        constexpr Double_t PARAM_ENERGY_MIN{0.10};
        const TVectorD vk{TMatrixDColumn_const(V, k)};
        const TVectorD vParam{vk.GetSub(0, nParams - 1)};
        const Double_t paramNorm2{vParam.Norm2Sqr()}, multNorm2{vk.Norm2Sqr() - paramNorm2};
        const Double_t paramFrac{paramNorm2 / std::max(paramNorm2 + multNorm2, SOFT_ZERO)};
        LOG(info) << "    --------------------------------------------------";
        LOG(info) << Form("  = weak mode #%d : eigenvalue = %.3e, param-subspace energy = %.1f%%", weakFound, sig(k),
                          100.0 * paramFrac);
        if (paramFrac < PARAM_ENERGY_MIN) {
            LOG(info) << "    --> multiplier-dominated, likely linearization artifact - skipped";
            continue;
        }
 
        const Double_t pNorm{std::sqrt(paramNorm2)};
 
        // Cosine similarity vs. user-defined constraints (in param subspace).
        // Evaluated before contributors so the print can be widened for unconstrained modes.
        constexpr Double_t OVERLAP_MIN{0.2};   // TODO: make this threshold configurable / meaningful
        Bool_t anyOverlap{kFALSE};
        if (aligner.Priors().ConstCount()) {
            std::vector<std::pair<Int_t, Double_t>> overlaps;   // (constraint #, cosine)
            overlaps.reserve(aligner.Priors().ConstCount());
            const TVectorD dots{cnstrRows * vParam};
            for (Int_t j = 0; j < aligner.Priors().ConstCount(); j++) {
                const Double_t cosine{std::fabs(dots[j]) / std::max(pNorm * cnstrNorms[j], SOFT_ZERO)};
                overlaps.emplace_back(j, cosine);
                if (cosine > OVERLAP_MIN)
                    anyOverlap = kTRUE;
            }
            LOG(info) << "    overlap with existing constraints (|cos angle|):";
            for (const auto& ov : overlaps)
                if (ov.second > OVERLAP_MIN)
                    LOG(info) << Form("      constraint row #%d : %.1f%% overlap", ov.first + 1, 100.0 * ov.second);
            LOG_IF(warning, !anyOverlap) << "  ! no overlap with any constraint detected"
                                            " -- likely new DoF to be fixed";
        }
 
        if (!anyOverlap) {
            // unconstrained mode: dump the full param eigenvector so the user
            // can read out the missing constraint directly. With dumpNew = kTRUE emit
            // a constraints-JSON snippet ready to paste into a constraints file.
            if (dumpNew) {
                openNewCnstrFile();
                if (newCnstrFile.is_open()) {
                    if (newCnstrCount > 0)
                        newCnstrFile << ",\n";
                    newCnstrFile << "    {\n";
                    newCnstrFile << Form("        \"Comment\": \"weak mode #%d, eigenvalue = %.3e, auto-dumped\",\n",
                                         weakFound, sig(k));
                    newCnstrFile << "        \"Coefficients\": [\n";
                    constexpr Double_t VALUE_SKIP_BELOW{0.5e-10};   // empirical threshold for export precision
                    for (Int_t r = 0; r < nParams; r++) {
                        const Double_t v{V(r, k) / pNorm};
                        if (std::fabs(v) < VALUE_SKIP_BELOW)
                            continue;   // loader will iterpret as zero
                        const ParamRef& p{aligner.Layout()[r]};
                        // a composite row has no detector to point at, so it is
                        // emitted the way the constraints loader addresses one
                        const TString quotedName{Form("\"%s\"", p.nodeName)};
                        newCnstrFile << (p.IsComposite()
                                             ? Form("            { \"NodeName\": %11s, \"ParamName\": \"%s\", "
                                                    "\"Value\": %.10f },\n",
                                                    quotedName.Data(), p.prmName, v)
                                             : Form("            { \"DetectorID\": %9d, \"ParamName\": \"%s\", "
                                                    "\"Value\": %.10f },\n",
                                                    aligner.DetectorModel()->EncodedID(p.leaf), p.prmName, v));
                    }
                    newCnstrFile << "            { \"RHS\": 0.0 }\n";
                    newCnstrFile << "        ]\n";
                    newCnstrFile << "    }";
                    newCnstrCount++;
                    LOG(info) << "    new constraint snippet written to " << newCnstrPath;
                }
            } else {
                LOG(info) << "    full eigenvector (param-subspace, normalized):";
                for (const ParamRef& p : aligner.Layout().All())
                    LOG(info) << Form("      %-30s : %+.6f", p.Label().Data(), V(p.full, k) / pNorm);
            }
        }
        if (anyOverlap || dumpNew) {
            // constrained mode: just the top entries by magnitude.
            std::vector<std::pair<Double_t, Int_t>> ranked;   // |entry|, row index
            ranked.reserve(nParams);
            for (Int_t r = 0; r < nParams; r++)
                ranked.emplace_back(std::fabs(V(r, k) / pNorm), r);
            std::partial_sort(ranked.begin(), ranked.begin() + std::min<Int_t>(BMN_GLOBAL_PARAMS_PD, ranked.size()),
                              ranked.end(),
                              [](const std::pair<Double_t, Int_t>& a, const std::pair<Double_t, Int_t>& b) {
                                  return a.first > b.first;
                              });
            LOG(info) << "    top contributors (param-subspace, normalized):";
            for (Int_t t = 0; t < BMN_GLOBAL_PARAMS_PD && t < (Int_t)ranked.size(); t++) {
                const Int_t row{ranked[t].second};
                const ParamRef& p{aligner.Layout()[row]};
                LOG(info) << Form("      %-30s : %+.6f", p.Label().Data(), V(row, k) / pNorm);
            }
        }
    }
 
    if (newCnstrFile.is_open()) {
        newCnstrFile << "\n]\n";
        newCnstrFile.close();
        LOG(info) << "    " << newCnstrCount << " new constraint(s) saved into " << newCnstrPath;
    }
    LOG(info) << "    ========================";
    LOG(info) << "    " << weakFound << " weak mode(s) detected below threshold";
    LOG(state) << " < BmnAlignReport : weak modes analysis done";
    LOG(state) << "   ---";
}
 
void BmnAlignReport::Correlations(const BmnAligner& aligner, Double_t absThreshold, Int_t topK, Bool_t skipIntraDet)
{
    LOG(state) << " > BmnAlignReport : parameter correlation analysis";
 
    if (aligner.Results().empty()) {
        LOG(error) << " <!> no iteration results available, aborting";
        return;
    }
    const BmnAlignResult& res{aligner.ResultLast()};
    const TMatrixDSym& C{res.Cov()};
    const TVectorD& err{res.AErr()};
    if (C.GetNrows() != aligner.Hierarchy().AllParams()) {
        LOG(error) << " <!> covariance matrix not available, aborting";
        return;
    }
 
    // sigma source: fAerr (full-length, inherited across iterations); 0 means "no info"
    auto sigmaOf = [&](Int_t g) -> Double_t {
        const Double_t s{err(g)};
        return (s > 0.0 && s != TMath::Infinity()) ? s : 0.0;
    };
    // both levels: a leaf-versus-parent-composite correlation is exactly the
    // pair the priors arbitrate, so it is the most informative one available
    const BmnParamLayout::Range_t& params{aligner.Layout().All()};
    Int_t validCnt{0};
    for (const ParamRef& p : params)
        if (sigmaOf(p.full) > 0.0)
            validCnt++;
 
    // collect upper-triangle pairs with |rho| >= absThreshold
    struct Pair
    {
        Double_t rho;
        const ParamRef *r, *s;
    };
    std::vector<Pair> pairs;
    Int_t totalPairs{0};
    for (std::size_t i = 0; i < params.size(); i++) {
        const ParamRef& r{params[i]};
        const Double_t sr{sigmaOf(r.full)};
        if (sr == 0.0)
            continue;
        for (std::size_t j = i + 1; j < params.size(); j++) {
            const ParamRef& s{params[j]};
            const Double_t ss{sigmaOf(s.full)};
            if (ss == 0.0)
                continue;
            if (skipIntraDet && r.node == s.node)
                continue;
            totalPairs++;
            const Double_t rho{C(r.full, s.full) / (sr * ss)};
            if (std::fabs(rho) >= absThreshold)
                pairs.push_back({rho, &r, &s});
        }
    }
 
    std::sort(pairs.begin(), pairs.end(),
              [](const Pair& a, const Pair& b) { return std::fabs(a.rho) > std::fabs(b.rho); });
 
    LOG(info) << Form("    %d param(s) with valid sigma", validCnt);
 
    // per-parameter global correlation
    const TVectorD& glCorr{res.AglCorr()};
    struct GlCorr
    {
        Double_t corr;
        const ParamRef* p;
    };
    std::vector<GlCorr> glAbove;
    for (const ParamRef& p : params) {
        if (sigmaOf(p.full) == 0.0)
            continue;
        if (glCorr(p.full) >= absThreshold)
            glAbove.push_back({glCorr(p.full), &p});
    }
    std::sort(glAbove.begin(), glAbove.end(), [](const GlCorr& a, const GlCorr& b) { return a.corr > b.corr; });
    LOG(info) << "    --------------------------------------------------";
    LOG(info) << Form("    %d param(s) with rho_gl >= %.2f out of %d valid", (Int_t)glAbove.size(), absThreshold,
                      validCnt);
    const Int_t glShown{std::min(topK, (Int_t)glAbove.size())};
    if (glShown > 0) {
        LOG(info) << "    top params by global correlation:";
        LOG(info) << Form("      %-51s   %-7s   %-10s", "param", "rho_gl", "sigma");
        for (Int_t k = 0; k < glShown; k++) {
            const ParamRef& p{*glAbove[k].p};
            LOG(info) << Form("      %s %-6s    %+.4f   %.2e", p.nodeName, p.prmName, glAbove[k].corr, sigmaOf(p.full));
        }
    }
 
    LOG(info) << "    --------------------------------------------------";
    LOG(info) << Form("    %d pair(s) with |rho| >= %.2f out of %d pair(s) considered", (Int_t)pairs.size(),
                      absThreshold, totalPairs);
    const Int_t shown{std::min(topK, (Int_t)pairs.size())};
    if (shown == 0) {
        LOG(state) << " < BmnAlignReport : correlation analysis done";
        LOG(state) << "   ---";
        return;
    }
    LOG(info) << "    top correlated pairs:";
    LOG(info) << Form("      %-51s  %-51s   %-7s   %-10s %-10s", "param i", "param j", "rho", "sigma_i", "sigma_j");
    for (Int_t k = 0; k < shown; k++) {
        const ParamRef &r{*pairs[k].r}, &s{*pairs[k].s};
        LOG(info) << Form("      %s %-6s   %s  %-6s   %+.4f   %.2e   %.2e", r.nodeName, r.prmName, s.nodeName,
                          s.prmName, pairs[k].rho, sigmaOf(r.full), sigmaOf(s.full));
    }
 
    LOG(state) << " < BmnAlignReport : correlation analysis done";
    LOG(state) << "   ---";
}
 
TH1* BmnAlignReport::ClonePlot(const TH1& source)
{
    TH1* copy{static_cast<TH1*>(source.Clone())};
    copy->SetDirectory(nullptr);
    return copy;
}
 
BmnAlignReport::Pads BmnAlignReport::PrepareCanvas()
{
    Pads pads;
    // prepare canvas and pads
    static const UInt_t ww{1600}, wh{1200};
    TCanvas* pCnv = new TCanvas(   // RAM is not released in macro mode
        "canvas", "BM@N alignment", ww, wh);
    pCnv->SetCanvasSize(ww - 20, wh - 70);
 
    // divide canvas into pads
    // pay attention: dynamic RAM is not released in ROOT macro mode
 
    // palce for synthetic histograms by X
    pads.resSXbefore = new TPad("padG1", "G1", 0.0, 0.0, 0.5, 0.34);
    pads.resSXafter = new TPad("padG2", "G2", 0.5, 0.0, 1.0, 0.34);
    // palce for synthetic histograms by Y
    pads.resSYbefore = new TPad("padG3", "G3", 0.0, 0.34, 0.5, 0.67);
    pads.resSYafter = new TPad("padG4", "G4", 0.5, 0.34, 1.0, 0.67);
    // palce for total residuals
    pads.resX = new TPad("padG5", "G5", 0.0, 0.67, 0.25, 1.0);
    pads.resY = new TPad("padG6", "G6", 0.25, 0.67, 0.5, 1.0);
    // palce for SVD singular values
    pads.svd = new TPad("padG7", "G7", 0.5, 0.67, 0.75, 1.0);
    // palce for 3D tracks visualization
    pads.tracks = new TPad("padG8", "G8", 0.75, 0.67, 1.0, 1.0);
 
    pads.resSXbefore->Draw();
    pads.resSXafter->Draw();
    pads.resSYbefore->Draw();
    pads.resSYafter->Draw();
    pads.resX->Draw();
    pads.resY->Draw();
    pads.svd->Draw();
    pads.tracks->Draw();
 
    return pads;
}
 
void BmnAlignReport::DrawResiduals(const BmnAligner& aligner, const Pads& pads, std::size_t idx)
{
    // BmnAlignResult hands its histograms out through non-const accessors
    const BmnAligner::Results_t& res{aligner.Results()};
    TH1* pHist;
 
    // top row: 1D total residuals by X (G20) and by Y (G21); reference is always [0]
 
    pads.resX->cd();
    pHist = ClonePlot(res[0].ResidualsX());
    pHist->SetStats(0);
    pHist->SetFillColor(kGray);
    pHist->SetLineColor(kGray);
    pHist->GetXaxis()->SetTitle("X residual, cm");
    pHist->GetYaxis()->SetTitle("measurements");
    pHist->GetYaxis()->SetRangeUser(0.0, res[idx].ResidualsX().GetMaximum() * 1.1);
    pHist->GetYaxis()->SetLabelSize(0);
    pHist->Draw("BAR");
    pHist = ClonePlot(res[idx].ResidualsX());
    pHist->SetStats(0);
    pHist->SetLineColor(kRed);
    pHist->Draw("CSAME");
 
    pads.resY->cd();
    pHist = ClonePlot(res[0].ResidualsY());
    pHist->SetStats(0);
    pHist->SetFillColor(kGray);
    pHist->SetLineColor(kGray);
    pHist->GetXaxis()->SetTitle("Y residual, cm");
    pHist->GetYaxis()->SetTitle("measurements");
    pHist->GetYaxis()->SetRangeUser(0.0, res[idx].ResidualsY().GetMaximum() * 1.1);
    pHist->GetYaxis()->SetLabelSize(0);
    pHist->Draw("BAR");
    pHist = ClonePlot(res[idx].ResidualsY());
    pHist->SetStats(0);
    pHist->SetLineColor(kRed);
    pHist->Draw("CSAME");
 
    // top-middle: SVD singular values for last iteration (G22), log Y
 
    pads.svd->cd();
    pads.svd->SetLogy();
    const TVectorD& vSig{res[idx].SVDSig()};
    const Int_t nSig{vSig.GetNrows()};
    if (nSig > 0) {
        // Identify the tail of smallest eigenvalues. SVD output is descending,
        // so we measure log-drops between consecutive pairs, take the two latest
        // "significant" drops, and start the plot BMN_SVD_TAIL_MARGIN points
        // before the earlier of them. Falls back to fewer drops (or to the last
        // BMN_SVD_TAIL_MIN points) if no clear jumps are found; tail length is
        // capped to no less than BMN_SVD_TAIL_MIN points either way.
        constexpr Int_t BMN_SVD_TAIL_MIN{75};
        constexpr Int_t BMN_SVD_TAIL_MARGIN{3};
        constexpr Int_t BMN_SVD_DROPS_KEEP{2};
        constexpr Double_t LOG2{0.6931471805599453};   // ln(2) - "factor of 2" threshold
        Int_t startIdx{0};
        if (nSig >= 2) {
            Double_t dropMax{0.0};
            VecDouble_t drop(nSig - 1, 0.0);
            for (Int_t i = 0; i < nSig - 1; i++) {
                Double_t a{std::max(vSig(i), SOFT_ZERO)}, b{std::max(vSig(i + 1), SOFT_ZERO)};
                drop[i] = std::log(a) - std::log(b);
                if (drop[i] > dropMax)
                    dropMax = drop[i];
            }
            const Double_t thr{std::max(0.5 * dropMax, LOG2)};
            Int_t earliestKeptJump{-1}, kept{0};
            for (Int_t i = nSig - 2; i >= 0 && kept < BMN_SVD_DROPS_KEEP; i--)
                if (drop[i] >= thr) {
                    earliestKeptJump = i;
                    kept++;
                }
            startIdx = (earliestKeptJump >= 0) ? std::max(0, earliestKeptJump - BMN_SVD_TAIL_MARGIN + 1)
                                               : std::max(0, nSig - BMN_SVD_TAIL_MIN);
        }
        if (nSig - startIdx < BMN_SVD_TAIL_MIN)
            startIdx = std::max(0, nSig - BMN_SVD_TAIL_MIN);
        const Int_t nTail{nSig - startIdx};
 
        // RAM not released - acceptable in ROOT macro mode (see PrepareCanvas)
        TGraph* pGsig = new TGraph(nTail);
        pGsig->SetName(Form("svdSig_%lu", res[idx].GetHistCounter()));
        pGsig->SetTitle("MNE SVD overview (tail);eigenvalue index;eigenvalue");
        for (Int_t i = 0; i < nTail; i++)
            pGsig->SetPoint(i, static_cast<Double_t>(startIdx + i), vSig(startIdx + i));
        pGsig->SetFillColor(kGray);
        pGsig->SetLineColor(kGray);
        pGsig->Draw("AB");
 
        // annotate the latest (minimal) eigenvalue above its bar
        const Double_t sigMin{vSig(nSig - 1)};
        TLatex* pSigMin = new TLatex();
        pSigMin->SetTextAlign(21);   // h-centered, v-bottom
        pSigMin->SetTextFont(62);    // bold Helvetica
        pSigMin->SetTextSize(0.05);
        pSigMin->SetTextColor(kRed);
        pSigMin->DrawLatex(static_cast<Double_t>(nSig - 1),
                           sigMin * 1.6,   // padded above the bar (pad is logy)
                           Form("#lambda_{min} = %.3g", sigMin));
    }
 
    // middle row: synthetic Y residuals over Z, before (G10) and after (G11)
 
    pads.resSYbefore->cd();
    pHist = ClonePlot(res[0].ResidualsSY());
    pHist->SetStats(0);
    pHist->SetTitle(Form("%s - before alignment", pHist->GetTitle()));
    pHist->Draw("COL");
 
    pads.resSYafter->cd();
    pHist = ClonePlot(res[idx].ResidualsSY());
    pHist->SetStats(0);
    pHist->SetTitle(Form("%s - after alignment", pHist->GetTitle()));
    pHist->Draw("COL");
 
    // bottom row: synthetic X residuals over Z, before (G00) and after (G01)
 
    pads.resSXbefore->cd();
    pHist = ClonePlot(res[0].ResidualsSX());
    pHist->SetStats(0);
    pHist->SetTitle(Form("%s - before alignment", pHist->GetTitle()));
    pHist->Draw("COL");
 
    pads.resSXafter->cd();
    pHist = ClonePlot(res[idx].ResidualsSX());
    pHist->SetStats(0);
    pHist->SetTitle(Form("%s - after alignment", pHist->GetTitle()));
    pHist->Draw("COL");
}
 
void BmnAlignReport::Draw(const BmnAligner& aligner, Int_t number)
{
    if (!aligner.Initialized())
        return;
    LOG(state) << " > BmnAlignReport : drawing started";
 
    const BmnAligner::Results_t& results{aligner.Results()};
    std::size_t resultsCnt{results.size()}, resultIdx{number < 0 ? number + resultsCnt : number - 1};
    LOG(info) << "    " << resultsCnt << " results available to draw";
    if (resultIdx < 0 || resultIdx >= resultsCnt) {
        LOG(error) << "<!> result with number " << number << " not found";
        LOG(state) << " < BmnAlignReport : drawing terminated";
        return;
    }
    LOG(info) << "    drawing " << std::abs(number) << " saved results counting from "
              << (number < 0 ? "the end" : "the beginning");
 
    const Pads pads{PrepareCanvas()};
    // TODO: 3D track samples must be optional or excluded
 
    // draw some tracks in 3D pad
    BmnRamIterator* pIterator{aligner.RamIterator()};
    pads.tracks->cd();
    const Double_t maxZdraw{aligner.MaxZvalue() * 1.1};
    const Double_t trackProbability{static_cast<Double_t>(DRAW_TRACKS_CNT)
                                    / static_cast<Double_t>(aligner.TracksUsable())};
    TRandom trueRnd(0);
    // RAM not released - acceptable in ROOT macro mode (see PrepareCanvas).
    // Lines must remain alive for the pad to keep drawing them on redraw.
    TPolyLine3D *pTrackLines = new TPolyLine3D[DRAW_TRACKS_CNT], *pApproxLines = new TPolyLine3D[DRAW_TRACKS_CNT];
 
    aligner.ProgressBar()->Init();
    Int_t iTrack{0}, iTrackDrawn{0};
 
    for (pIterator->ResetAll(); iTrackDrawn < DRAW_TRACKS_CNT && !pIterator->EndOfTracks();
         pIterator->NextTrack(), iTrack++)
    {
        // skip tracks according to probability
        if (!aligner.TrackUsable()->at(iTrack) || trueRnd.Uniform() > trackProbability)
            continue;
 
        // draw track itself
        Int_t iHit{0};
        for (; !pIterator->EndOfHits(); pIterator->NextHit(), iHit++) {
            // x and z are swapped for visual convenience
            pTrackLines[iTrackDrawn].SetPoint(iHit, pIterator->HitZ(), pIterator->HitY(), pIterator->HitX());
        }
        pTrackLines[iTrackDrawn].SetLineColor(kBlue);
        pTrackLines[iTrackDrawn].Draw();
 
        // draw linear approximation of the track
        const SVectLC& Alpha{aligner.ResultLast().Alpha(iTrack)};
        pApproxLines[iTrackDrawn].SetPoint(0, 0.0, Alpha(iY0), Alpha(iX0));
        pApproxLines[iTrackDrawn].SetPoint(1, maxZdraw, Alpha(iY0) + maxZdraw * Alpha(iY1),
                                           Alpha(iX0) + maxZdraw * Alpha(iX1));
        pApproxLines[iTrackDrawn].SetLineColor(kRed);
        pApproxLines[iTrackDrawn].SetLineStyle(kDotted);
        pApproxLines[iTrackDrawn].Draw();
        iTrackDrawn++;
    }   // loop over tracks
    aligner.ProgressBar()->Clear();
 
    DrawResiduals(aligner, pads, resultIdx);
 
    LOG(state) << " < BmnAlignReport : drawing finished";
    LOG(state) << "   ---";
}