histoGL.C

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:     Claude (Opus 4.8)
    Supervisor: Igor Polev, polev@jinr.int
 
    Histogramming method of planar-detector alignment - WHOLE-DETECTOR view.
 
    Companion to histo.C. Where histo.C fills and reads the residual-histogram
    set separately for every alignable element, histoGL.C pools ALL hits of the
    whole detector into a single set of histograms - it reports the ensemble
    first-order corrections rather than per-element ones.
 
    The reading rules are identical to histo.C (see there for the derivation):
 
        histogram / profile of rx     read as         gives
        --------------------------    ----------      -----------------
        rx                            mean            dX
        rx  vs  beta_x                slope           dZ    = -slope
        rx  vs  y_bar                 slope           psi   = -slope
        rx  vs  beta_x * y_bar        slope           phi   = -slope
        rx  vs  beta_x * x_bar        slope           theta = +slope
 
    (dY and the alternative dZ/phi/theta/psi come from the symmetric ry
     expansion, exactly as in histo.C.)
 
    Two extra panels fill the top row: the running evolution of ONE selected
    correction as more and more tracks are accounted for - the X-solution on
    the left, the Y-solution on the right. Watching the two independent
    readings settle (and, for the shared parameters, converge toward each
    other) is a direct check of how much data the ensemble reading needs and
    whether the single-covariate approximation is self-consistent.
 
    Interactive: the Left / Right arrows no longer browse detector elements
    (there is only one, pooled histogram set now) - they cycle which
    correction the two evolution panels display. q or Esc quits.
 
    Run on a DST reconstructed with the current best alignment; the numbers
    printed are the ensemble corrections still to be applied.
*/
 
// headers below are for IDE syntax to work and they will be
// turned off during execution in ROOT macro mode
#ifndef __CLING__
#include "BmnAlignModels.h"   // IWYU pragma: keep
#include "BmnRootIterator.h"
#include "BmnSimpleProgressBar.h"
#include "BmnWeightedLine.h"
 
#include <TCanvas.h>
#include <TF1.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TLine.h>
#include <TPad.h>
#include <TPaveText.h>
#include <TProfile.h>
#include <TROOT.h>
#include <TString.h>
#include <TStyle.h>
#include <TSystem.h>
#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
#endif
 
// POSIX terminal control for single-key (arrow) interactive parameter cycling.
// Kept out of the CLING-skipped block: these are not autoloaded by ROOT.
#include <termios.h>
#include <unistd.h>
 
// correction parameters, ordered as BmnParamLayout::ParamNames().
// Local enum so the macro stays self-contained in interpreted ROOT mode
#define PARAM_CNT 6
enum eParam
{
    P_dX,
    P_dY,
    P_dZ,
    P_phi,
    P_theta,
    P_psi
};
constexpr static const char* PARAM_NAMES[PARAM_CNT]{"dX", "dY", "dZ", "phi", "theta", "psi"};
constexpr static const char* PARAM_UNIT[PARAM_CNT]{"cm", "cm", "cm", "rad", "rad", "rad"};
// how each parameter is read from its diagnostic: mean of residual (shift)
// or slope of residual vs covariate, and the sign that turns the reading
// into the correction (matches the rx expansion in the header)
constexpr static const Bool_t PARAM_IS_SHIFT[PARAM_CNT]{kTRUE, kTRUE, kFALSE, kFALSE, kFALSE, kFALSE};
// two independent readings per rotation/dZ: from the X-residual and from the
// Y-residual. Only the sign of psi differs between them (see dXdA / dYdA).
constexpr static const Double_t PARAM_SIGN_X[PARAM_CNT]{+1.0, +1.0, -1.0, -1.0, +1.0, -1.0};
constexpr static const Double_t PARAM_SIGN_Y[PARAM_CNT]{+1.0, +1.0, -1.0, -1.0, +1.0, +1.0};
constexpr static const char* PARAM_DIAG[PARAM_CNT]{// X-residual reading
                                                   "mean of X-residual",
                                                   "mean of Y-residual",
                                                   "slope of X-res vs track x-slope",
                                                   "slope of X-res vs (x-slope * y)",
                                                   "slope of X-res vs (x-slope * x)",
                                                   "slope of X-res vs y"};
constexpr static const char* PARAM_DIAG_Y[PARAM_CNT]{// Y-residual reading
                                                     "mean of X-residual",
                                                     "mean of Y-residual",
                                                     "slope of Y-res vs track y-slope",
                                                     "slope of Y-res vs (y-slope * y)",
                                                     "slope of Y-res vs (y-slope * x)",
                                                     "slope of Y-res vs x"};
constexpr static const char* COV_TITLE_X[PARAM_CNT]{"X residual, cm",
                                                    "Y residual, cm",
                                                    "track x-slope #beta_{x}",
                                                    "#beta_{x} #upoint y_{bar}, cm",
                                                    "#beta_{x} #upoint x_{bar}, cm",
                                                    "y_{bar}, cm"};
constexpr static const char* COV_TITLE_Y[PARAM_CNT]{"X residual, cm",
                                                    "Y residual, cm",
                                                    "track y-slope #beta_{y}",
                                                    "#beta_{y} #upoint y_{bar}, cm",
                                                    "#beta_{y} #upoint x_{bar}, cm",
                                                    "x_{bar}, cm"};
 
// slope-panel drawing mode (drawMode argument):
//   DM_PROFILE - both slope rows as 1D TProfile  (X in middle, Y in bottom)
//   DM_SCATTER - both slope rows as 2D TH2D       (X in middle, Y in bottom)
//   DM_MIX     - X-residual slopes only: 1D TProfile in middle, 2D TH2D in
//                bottom (Y-residual readings are computed but not drawn)
enum eDrawMode
{
    DM_PROFILE,
    DM_SCATTER,
    DM_MIX
};
 
// canvas and pads configuration: top row = shifts (dX, dY) + two evolution
// panels, middle row = X-residual slopes, bottom row = Y-residual (alternative)
// slopes; one param per column
#define RES_X 3440
#define RES_Y 1440
#define PADS_X 4
#define PADS_Y 3
 
// residual histogram ranges for the drawn shift panels
#define RESIDUAL_RANGE_X 0.02
#define RESIDUAL_RANGE_Y 0.4
#define RESIDUAL_BINS 201
#define PROFILE_BINS 201
#define PAVE_TEXT_SIZE 0.04   // ~half of the TPaveText auto-fit default
 
// evolution sampling: record every track for the first EVO_DENSE_HEAD tracks
// (to resolve the fast early convergence), then thin out to about
// EVO_TAIL_POINTS samples across the remaining range
#define EVO_DENSE_HEAD 200
#define EVO_TAIL_POINTS 600
 
// simple accumulator: mean of y and slope of a straight line y = a + b*x, unweighted
struct LinAcc
{
    BmnWeightedLine line;   // unit weights
    void Add(Double_t x, Double_t y) { line.Add(x, y); }
    Double_t MeanY() const { return line.N() > 0 ? line.Mean() : 0.0; }
    Double_t RmsY() const { return line.N() < 2 ? 0.0 : line.Rms(); }
    Double_t MeanYErr() const { return line.N() > 0 ? RmsY() / std::sqrt(line.N()) : 0.0; }
    Double_t Slope() const { return std::abs(line.Det()) > 0.0 ? line.Slope() : 0.0; }
    Double_t SlopeErr() const
    {
        if (!(std::abs(line.Det()) > 0.0) || line.N() < 3)
            return 0.0;
        // unit weights: the variance scale comes from the scatter itself
        Double_t s2{line.Chi2() / (line.N() - 2)};
        return s2 > 0.0 ? std::sqrt(s2 * line.SlopeVar()) : 0.0;
    }
    // reading turned into a correction: sign array selects X- or Y-solution
    Double_t Value(Int_t p, const Double_t* sign) const { return sign[p] * (PARAM_IS_SHIFT[p] ? MeanY() : Slope()); }
    Double_t Error(Int_t p) const { return PARAM_IS_SHIFT[p] ? MeanYErr() : SlopeErr(); }
};
 
// two independent readings for the whole detector: from X- and from Y-residuals
struct DetAcc
{
    std::array<LinAcc, PARAM_CNT> x, y;
};
 
// per-hit drawn quantities for the whole detector, pooled over all elements
struct DetDraw
{
    VecDouble_t rx, ry, betaX, betaXY, betaXX, yBar,   // X-residual covariates
        betaY, betaYY, betaYX, xBar;                   // Y-residual covariates
};
 
// running-parameter history: one value/error series per parameter and solution,
// sampled at the track counts stored in n
struct EvoHist
{
    VecDouble_t n;   // cumulative tracks at each sample
    std::array<VecDouble_t, PARAM_CNT> valX, errX, valY, errY;
    void Sample(Long64_t tracks, const DetAcc& a)
    {
        n.push_back(static_cast<Double_t>(tracks));
        for (Int_t p = 0; p < PARAM_CNT; p++) {
            valX[p].push_back(a.x[p].Value(p, PARAM_SIGN_X));
            errX[p].push_back(a.x[p].Error(p));
            valY[p].push_back(a.y[p].Value(p, PARAM_SIGN_Y));
            errY[p].push_back(a.y[p].Error(p));
        }
    }
};
 
// put the controlling terminal into non-canonical, no-echo mode (on == true)
// so single keystrokes - including arrow escape sequences - are read at once
static void histoGLSetRawMode(bool on, struct termios* saved)
{
    if (on) {
        tcgetattr(STDIN_FILENO, saved);
        struct termios raw{*saved};
        raw.c_lflag &= ~(ICANON | ECHO);
        raw.c_cc[VMIN] = 0;   // non-blocking read
        raw.c_cc[VTIME] = 0;
        tcsetattr(STDIN_FILENO, TCSANOW, &raw);
    } else {
        tcsetattr(STDIN_FILENO, TCSANOW, saved);
    }
}
 
// one drawn diagnostic panel (a shift histogram or a slope profile / 2D scatter)
struct Panel
{
    Int_t row, col, p;
    Bool_t isShift, isY, as2D;   // as2D: draw as TH2D instead of TProfile
    const VecDouble_t* res;      // residual values
    const VecDouble_t* cov;      // covariate values (nullptr for shift)
    Double_t resRange;
    const char* diag;
    const char* covTitle;
    const LinAcc* a;   // accumulator for the value/error box
    const Double_t* sign;
};
 
/**
 * @brief Whole-detector histogramming estimate of the alignment corrections.
 *
 * Companion to histo.C: the same reading rules, but every hit of the detector
 * is pooled into one histogram set, so the numbers are ensemble corrections
 * rather than per-element ones. Useful for judging a global shift or rotation
 * of the whole assembly.
 *
 * @param fileListPath text file listing the .root files to read
 * @param selParam     correction shown in the two evolution panels, one of the
 *                     P_dX .. P_psi constants
 * @param drawMode     slope panel style: DM_PROFILE, DM_SCATTER or DM_MIX
 * @param pdfPath      .pdf file to save the canvas into, null to only display
 *
 * @return 0 on success, non-zero if the input data could not be read
 */
int histoGL(const char* fileListPath = "/home/igor/bmnroot/align_new/cfg/run8/7466.noa_rlerr.list",
            Int_t selParam = P_dZ,     // correction shown in the two evolution panels
            Int_t drawMode = DM_MIX,   // slope panels: DM_PROFILE / DM_SCATTER / DM_MIX
            const char* pdfPath = nullptr)
{
    // init section
    BmnRun8cgTopDM detModel;
    detModel.Init();
    if (selParam < 0 || selParam >= PARAM_CNT) {
        std::cerr << " <!> histoGL: selParam " << selParam << " out of range [0, " << PARAM_CNT << ")" << std::endl;
        return 1;
    }
    if (drawMode < DM_PROFILE || drawMode > DM_MIX) {
        std::cerr << " <!> histoGL: drawMode " << drawMode << " out of range [" << DM_PROFILE << ", " << DM_MIX << "]"
                  << std::endl;
        return 1;
    }
    BmnSimpleProgressBar progressBar;
    BmnRootIterator itr(fileListPath, detModel, BmnDefaultTrackFilter(), BmnDefaultHitFilter(detModel), &progressBar);
    if (!itr.Initialized()) {
        std::cerr << " <!> histoGL: failed to initialize iterator on " << fileListPath << std::endl;
        return 1;
    }
 
    std::cout << "\nInitializing ROOT reader for " << fileListPath << "\n";
    itr.CountElements();   // count suitable tracks/hits, calibrate per-track progress
    const Long64_t totalTracks{itr.GetTotalTracks()};
    // thin the evolution tail to about EVO_TAIL_POINTS samples
    const Long64_t evoStride{std::max<Long64_t>(1, totalTracks / EVO_TAIL_POINTS)};
 
    DetAcc acc;     // pooled X- and Y-residual accumulators (whole detector)
    DetDraw draw;   // pooled per-hit records for drawing (whole detector)
    EvoHist evo;    // running-correction history vs tracks accounted
 
    // walk tracks; fit a straight line per track, then read residuals
    std::cout << "Reading data...\n";
    Long64_t usedTracks{0}, usedHits{0};
    progressBar.Init();
    itr.ResetAll();
    while (!itr.EndOfTracks()) {
        // weighted straight-line fit  x = alpha_x + beta_x * z,  y = alpha_y + beta_y * z
        BmnWeightedLine lineX, lineY;
        itr.ResetHits();
        while (!itr.EndOfHits()) {
            lineX.Add(itr.HitZ(), itr.HitX(), itr.HitWx());
            lineY.Add(itr.HitZ(), itr.HitY(), itr.HitWy());
            itr.NextHit();
        }
        if (!(lineX.Det() > 0.0) || !(lineY.Det() > 0.0)) {   // degenerate track (no z spread)
            itr.NextTrack();
            continue;
        }
        Double_t betaX{lineX.Slope()}, alphaX{lineX.At0()}, betaY{lineY.Slope()}, alphaY{lineY.At0()};
 
        // residuals and covariates (bar = straight-line prediction),
        // pooled over the whole detector (no per-element split)
        itr.ResetHits();
        while (!itr.EndOfHits()) {
            Double_t z{itr.HitZ()}, xBar{alphaX + betaX * z}, yBar{alphaY + betaY * z}, rx{itr.HitX() - xBar},
                ry{itr.HitY() - yBar};
 
            // X-residual readings
            acc.x[P_dX].Add(0.0, rx);               // mean  -> dX
            acc.x[P_dZ].Add(betaX, rx);             // slope -> -dZ
            acc.x[P_phi].Add(betaX * yBar, rx);     // slope -> -phi
            acc.x[P_theta].Add(betaX * xBar, rx);   // slope -> +theta
            acc.x[P_psi].Add(yBar, rx);             // slope -> -psi
            // Y-residual readings (alternative solution)
            acc.y[P_dY].Add(0.0, ry);               // mean  -> dY
            acc.y[P_dZ].Add(betaY, ry);             // slope -> -dZ
            acc.y[P_phi].Add(betaY * yBar, ry);     // slope -> -phi
            acc.y[P_theta].Add(betaY * xBar, ry);   // slope -> +theta
            acc.y[P_psi].Add(xBar, ry);             // slope -> +psi
 
            draw.rx.push_back(rx);
            draw.ry.push_back(ry);
            draw.betaX.push_back(betaX);
            draw.betaXY.push_back(betaX * yBar);
            draw.betaXX.push_back(betaX * xBar);
            draw.yBar.push_back(yBar);
            draw.betaY.push_back(betaY);
            draw.betaYY.push_back(betaY * yBar);
            draw.betaYX.push_back(betaY * xBar);
            draw.xBar.push_back(xBar);
            usedHits++;
            itr.NextHit();
        }
        usedTracks++;
        // sample the running corrections: dense at the start, thinned later
        if (usedTracks <= EVO_DENSE_HEAD || usedTracks % evoStride == 0)
            evo.Sample(usedTracks, acc);
        itr.NextTrack();
    }
    progressBar.Clear();
 
    if (usedTracks == 0) {
        std::cerr << " <!> histoGL: no usable tracks found" << std::endl;
        return 1;
    }
    // guarantee the final state is the last evolution sample
    if (evo.n.empty() || evo.n.back() != static_cast<Double_t>(usedTracks))
        evo.Sample(usedTracks, acc);
 
    // numeric report: single ensemble reading for the whole detector
 
    using namespace std;
    cout << "\nWhole-detector histogramming: ensemble first-order corrections\n";
    cout << usedTracks << " tracks, " << usedHits << " hits\n\n";
 
    cout << fixed << setprecision(6);
    cout << "  " << left << setw(7) << "dX" << " = " << right << setw(11) << acc.x[P_dX].Value(P_dX, PARAM_SIGN_X)
         << " +/- " << setw(10) << acc.x[P_dX].Error(P_dX) << "   (" << PARAM_DIAG[P_dX] << ")\n";
    cout << "  " << left << setw(7) << "dY" << " = " << right << setw(11) << acc.y[P_dY].Value(P_dY, PARAM_SIGN_Y)
         << " +/- " << setw(10) << acc.y[P_dY].Error(P_dY) << "   (" << PARAM_DIAG[P_dY] << ")\n";
    for (Int_t p = P_dZ; p < PARAM_CNT; p++)
        cout << "  " << left << setw(7) << PARAM_NAMES[p] << right << "  X: " << setw(11)
             << acc.x[p].Value(p, PARAM_SIGN_X) << " +/- " << setw(9) << acc.x[p].Error(p) << "   Y: " << setw(11)
             << acc.y[p].Value(p, PARAM_SIGN_Y) << " +/- " << setw(9) << acc.y[p].Error(p) << "\n";
    cout << "  RMS(rx) = " << acc.x[P_dX].RmsY() << "   RMS(ry) = " << acc.y[P_dY].RmsY() << "\n" << flush;
 
    // drawn objects are detached from any TDirectory and marked kCanDelete so
    // that a pad Clear() on each redraw disposes of the previous picture
    TH1::AddDirectory(kFALSE);
    gStyle->SetOptFit(1);
    gStyle->SetOptStat("emr");
 
    TCanvas* pCnv = new TCanvas(   // RAM not released - acceptable in macro mode
        "canvas", "BM@N histogramming (whole detector)", RES_X, RES_Y);
    pCnv->SetCanvasSize(RES_X - 20, RES_Y - 96);
 
    // build the pad grid once; the static diagnostics never change while the
    // user cycles the selected parameter, so only the two evolution pads are
    // redrawn on a key press
    pCnv->cd();
    const Double_t padW{1.0 / PADS_X}, padH{1.0 / PADS_Y};
    TPad* pPads[PADS_Y][PADS_X];
    for (UInt_t i = 0; i < PADS_Y; i++)
        for (UInt_t j = 0; j < PADS_X; j++) {
            pPads[i][j] = new TPad(Form("pad%u%u", i, j), Form("pad%u%u", i, j), padW * j, padH * i, padW * (j + 1),
                                   padH * (i + 1));
            pPads[i][j]->Draw();
        }
    // top row is physical index PADS_Y-1; its two rightmost pads host the
    // evolution panels (X-solution then Y-solution)
    TPad* pEvoX{pPads[PADS_Y - 1][2]};
    TPad* pEvoY{pPads[PADS_Y - 1][3]};
 
    // static panels: top row = shifts (dX, dY), middle = X-residual slopes,
    // bottom = Y-residual (alternative) slopes; one param per column - exactly
    // the histo.C set, but pooled over the whole detector
    const Int_t slopeP[4]{P_dZ, P_phi, P_theta, P_psi};
    const VecDouble_t* covX[4]{&draw.betaX, &draw.betaXY, &draw.betaXX, &draw.yBar};
    const VecDouble_t* covY[4]{&draw.betaY, &draw.betaYY, &draw.betaYX, &draw.xBar};
    std::vector<Panel> panels;
    panels.push_back({0, 0, P_dX, kTRUE, kFALSE, kFALSE, &draw.rx, nullptr, RESIDUAL_RANGE_X, PARAM_DIAG[P_dX],
                      COV_TITLE_X[P_dX], &acc.x[P_dX], PARAM_SIGN_X});
    panels.push_back({0, 1, P_dY, kTRUE, kTRUE, kFALSE, &draw.ry, nullptr, RESIDUAL_RANGE_Y, PARAM_DIAG[P_dY],
                      COV_TITLE_Y[P_dY], &acc.y[P_dY], PARAM_SIGN_Y});
    for (Int_t c = 0; c < 4; c++) {
        Int_t p{slopeP[c]};
        // middle row: X-residual slope, 2D only in pure-scatter mode
        panels.push_back({1, c, p, kFALSE, kFALSE, drawMode == DM_SCATTER, &draw.rx, covX[c], RESIDUAL_RANGE_X,
                          PARAM_DIAG[p], COV_TITLE_X[p], &acc.x[p], PARAM_SIGN_X});
        // bottom row: Y-residual slope normally; in mix mode the same X-residual
        // slope drawn as 2D (Y readings stay computed, not drawn)
        if (drawMode == DM_MIX)
            panels.push_back({2, c, p, kFALSE, kFALSE, kTRUE, &draw.rx, covX[c], RESIDUAL_RANGE_X, PARAM_DIAG[p],
                              COV_TITLE_X[p], &acc.x[p], PARAM_SIGN_X});
        else
            panels.push_back({2, c, p, kFALSE, kTRUE, drawMode == DM_SCATTER, &draw.ry, covY[c], RESIDUAL_RANGE_Y,
                              PARAM_DIAG_Y[p], COV_TITLE_Y[p], &acc.y[p], PARAM_SIGN_Y});
    }
 
    for (const Panel& pn : panels) {
        pPads[PADS_Y - 1 - pn.row][pn.col]->cd();
        const char* tag{pn.isShift ? "" : (pn.isY ? " (Y)" : " (X)")};
        const char* suff{pn.isY ? "y" : "x"};
        const char* axis{pn.isY ? "Y" : "X"};
 
        // per-pad report box with the extracted value
        TPaveText* pave = new TPaveText(0.70, 0.86, 0.99, 0.92, "BR NDC");
        pave->SetBit(kCanDelete);
        pave->SetFillColor(0);
        pave->SetTextSize(PAVE_TEXT_SIZE);
        pave->AddText(
            Form("%s%s = %+.5f #pm %.5f", PARAM_NAMES[pn.p], tag, pn.a->Value(pn.p, pn.sign), pn.a->Error(pn.p)));
 
        if (pn.isShift) {
            // shift: 1D residual histogram, correction = mean (dashed line)
            TH1D* h = new TH1D(Form("h_%s", PARAM_NAMES[pn.p]), Form("%s;%s;hits", pn.diag, pn.covTitle), RESIDUAL_BINS,
                               -pn.resRange, pn.resRange);
            h->SetBit(kCanDelete);
            for (Double_t r : *pn.res)
                h->Fill(r);
            h->SetLineColor(kAzure + 2);
            h->SetStats(0);
            h->Draw();
            TLine* mean = new TLine(h->GetMean(), 0.0, h->GetMean(), h->GetMaximum());
            mean->SetBit(kCanDelete);
            mean->SetLineColor(kRed);
            mean->SetLineStyle(2);
            mean->SetLineWidth(2);
            mean->Draw();
        } else {
            // rotation / dZ: residual against its covariate, slope fit
            const VecDouble_t& cov{*pn.cov};
            const VecDouble_t& res{*pn.res};
            Double_t lo{+1e30}, hi{-1e30};
            for (Double_t c : cov) {
                if (c < lo)
                    lo = c;
                if (c > hi)
                    hi = c;
            }
            if (hi <= lo) {
                lo -= 1e-3;
                hi += 1e-3;
            }
 
            Double_t margin{0.05 * (hi - lo)};
            if (pn.as2D) {
                // 2D scatter of residual vs covariate; pol1 fit is the slope
                // (TH2::Fit weights each column by content, like the profile)
                TH2D* h2 = new TH2D(Form("h2_%s_%s", PARAM_NAMES[pn.p], suff),
                                    Form("%s;%s;%s residual, cm", pn.diag, pn.covTitle, axis), PROFILE_BINS,
                                    lo - margin, hi + margin, RESIDUAL_BINS, -pn.resRange, pn.resRange);
                h2->SetBit(kCanDelete);
                for (std::size_t i = 0; i < cov.size(); i++)
                    h2->Fill(cov[i], res[i]);
                h2->Fit("pol1", "Q");
                if (TF1* fit = h2->GetFunction("pol1")) {
                    fit->SetLineColor(kRed);
                    fit->SetLineWidth(2);
                }
                h2->SetStats(0);
                h2->Draw("COLZ");
            } else {
                // residual profiled against its covariate, slope fit
                TProfile* pr = new TProfile(Form("p_%s_%s", PARAM_NAMES[pn.p], suff),
                                            Form("%s;%s;#LT %s residual #GT, cm", pn.diag, pn.covTitle, axis),
                                            PROFILE_BINS, lo - margin, hi + margin);
                pr->SetBit(kCanDelete);
                for (std::size_t i = 0; i < cov.size(); i++)
                    pr->Fill(cov[i], res[i]);
                pr->SetMarkerStyle(kFullCircle);
                pr->SetMarkerColor(kAzure + 2);
                pr->SetLineColor(kAzure + 2);
                pr->Fit("pol1", "Q");
                pr->SetStats(0);
                pr->Draw();
            }
        }
        pave->Draw();
    }
 
    // (re)draw one evolution panel: the running value of parameter p (X- or
    // Y-solution) with a statistical-error band, against the number of tracks
    auto drawEvo = [&](TPad* pad, Int_t p, Bool_t useY) {
        pad->cd();
        pad->Clear();
        pad->SetLogx();
        pad->SetGridx();
        pad->SetGridy();
 
        const LinAcc& a{useY ? acc.y[p] : acc.x[p]};
        const char* sol{useY ? "Y-solution" : "X-solution"};
        const VecDouble_t& vv{useY ? evo.valY[p] : evo.valX[p]};
        const VecDouble_t& ve{useY ? evo.errY[p] : evo.errX[p]};
 
        // dX has no independent Y reading and dY none in X - say so plainly
        if (a.line.N() <= 0 || vv.empty()) {
            TPaveText* msg = new TPaveText(0.05, 0.42, 0.95, 0.58, "BR NDC");
            msg->SetBit(kCanDelete);
            msg->SetFillColor(0);
            msg->SetTextSize(PAVE_TEXT_SIZE);
            msg->AddText(Form("%s : no independent %s", PARAM_NAMES[p], sol));
            msg->Draw();
            pad->Modified();
            pad->Update();
            return;
        }
 
        const Int_t n{static_cast<Int_t>(vv.size())};
        TGraphErrors* band = new TGraphErrors(n);   // value with +/- error band
        band->SetBit(kCanDelete);
        TGraph* line = new TGraph(n);   // central value line
        line->SetBit(kCanDelete);
        for (Int_t i = 0; i < n; i++) {
            band->SetPoint(i, evo.n[i], vv[i]);
            band->SetPointError(i, 0.0, ve[i]);
            line->SetPoint(i, evo.n[i], vv[i]);
        }
        const Double_t vFinal{vv.back()}, eFinal{ve.back()};
 
        // scale to the converged tail so the plateau - not the noisy first few
        // tracks - fills the panel
        Double_t ylo{+1e30}, yhi{-1e30};
        for (Int_t i = n / 2; i < n; i++) {
            ylo = std::min(ylo, vv[i] - ve[i]);
            yhi = std::max(yhi, vv[i] + ve[i]);
        }
        if (!(yhi > ylo)) {
            ylo = vFinal - 1e-6;
            yhi = vFinal + 1e-6;
        }
        const Double_t span{0.25 * (yhi - ylo)};
        ylo -= span;
        yhi += span;
 
        band->SetTitle(Form("%s  (%s);tracks accounted;%s, %s", PARAM_NAMES[p], sol, PARAM_NAMES[p], PARAM_UNIT[p]));
        band->SetFillColorAlpha(kAzure - 9, 0.7);
        band->SetLineColor(kAzure - 9);
        band->SetMinimum(ylo);
        band->SetMaximum(yhi);
        band->Draw("A3");
        band->GetXaxis()->SetLimits(evo.n.front(), evo.n.back());
 
        line->SetLineColor(kAzure + 2);
        line->SetLineWidth(2);
        line->Draw("L SAME");
 
        // final value as a red dashed reference
        TLine* ref = new TLine(evo.n.front(), vFinal, evo.n.back(), vFinal);
        ref->SetBit(kCanDelete);
        ref->SetLineColor(kRed);
        ref->SetLineStyle(2);
        ref->SetLineWidth(2);
        ref->Draw();
 
        TPaveText* box = new TPaveText(0.50, 0.82, 0.99, 0.92, "BR NDC");
        box->SetBit(kCanDelete);
        box->SetFillColor(0);
        box->SetTextSize(PAVE_TEXT_SIZE);
        box->AddText(Form("%s = %+.5f #pm %.5f", PARAM_NAMES[p], vFinal, eFinal));
        box->Draw();
 
        pad->Modified();
        pad->Update();
    };
 
    // redraw both evolution panels for the currently selected parameter
    auto showEvo = [&](Int_t p) {
        drawEvo(pEvoX, p, kFALSE);
        drawEvo(pEvoY, p, kTRUE);
        pCnv->Update();
    };
 
    showEvo(selParam);
 
    if (pdfPath && pdfPath[0])
        pCnv->SaveAs(pdfPath);
 
    // interactive: Left/Right arrows cycle the selected correction shown in the
    // two evolution panels (wrapping around), q or Esc quits. Skipped in batch
    // mode or when stdin is not a terminal (e.g. piped input) - there the
    // default parameter is drawn once.
    if (!gROOT->IsBatch() && isatty(STDIN_FILENO)) {
        cout << "\n[ Right / Left arrows: next / previous correction"
                " (evolution panels) | q: quit ]\n"
             << flush;
        struct termios saved;
        histoGLSetRawMode(true, &saved);
        Bool_t running{kTRUE};
        while (running) {
            gSystem->ProcessEvents();
            // leave if the user closed the canvas window
            if (!gROOT->GetListOfCanvases()->FindObject(pCnv))
                break;
            unsigned char c{0};
            if (read(STDIN_FILENO, &c, 1) <= 0) {   // nothing pending
                gSystem->Sleep(20);
                continue;
            }
            if (c == 'q' || c == 'Q') {
                running = kFALSE;
            } else if (c == 27) {   // ESC: arrow escape sequence or a bare Esc
                unsigned char a{0}, b{0};
                gSystem->Sleep(1);   // let the rest of the sequence arrive
                if (read(STDIN_FILENO, &a, 1) <= 0) {
                    running = kFALSE;
                    continue;
                }
                read(STDIN_FILENO, &b, 1);
                if (a == '[' && b == 'C')
                    selParam = (selParam + 1) % PARAM_CNT;
                else if (a == '[' && b == 'D')
                    selParam = (selParam - 1 + PARAM_CNT) % PARAM_CNT;
                else
                    continue;
                cout << "  selected correction: " << PARAM_NAMES[selParam] << "\n" << flush;
                showEvo(selParam);
            }
        }
        histoGLSetRawMode(false, &saved);
    }
 
    return 0;
}