histo.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 (diagnostic macro).
 
    For every alignable element the classic residual-histogram set is filled
    and read to extract a first-order estimate of each of the 6 corrections.
    The reading rules follow directly from the linearized measurement model
    (see BmnMeasureModel::dXdA / dYdA), where the X-residual expands as
 
        rx = dX - dZ*beta_x - psi*y_bar - phi*(beta_x*y_bar) + theta*(beta_x*x_bar)
 
    so each correction is the coefficient of one covariate:
 
        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 is read from the mean of ry; the symmetric ry expansion gives an
     independent estimate of dZ, psi, phi, theta - reported as the alternative
     readings and drawn as the X/Y pairs next to the X-residual ones.)
 
    The single-covariate reading is deliberately approximate: the covariates
    are mildly correlated, so the method is meant to be iterated - align the
    reconstruction with the extracted corrections and re-run until every mean
    sits at 0 and every profile is flat.
 
    Run on a DST reconstructed with the current best alignment; the numbers
    printed are the 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 <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 <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
#endif
 
// POSIX terminal control for single-key (arrow) interactive detector switching.
// 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"};
// 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, 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
 
// 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 per detector: from X-residuals and from Y-residuals
struct DetAcc
{
    std::array<LinAcc, PARAM_CNT> x, y;
};
 
// per-hit drawn quantities kept for every detector so the picture can be
// redrawn instantly when the user browses to another detector (no re-read)
struct DetDraw
{
    VecDouble_t rx, ry, betaX, betaXY, betaXX, yBar,   // X-residual covariates
        betaY, betaYY, betaYX, xBar;                   // Y-residual covariates
};
 
// 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 histoSetRawMode(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 Per-element histogramming estimate of the alignment corrections.
 *
 * Diagnostic alternative to the global fit: for one alignable element at a time
 * it fills the classic residual-histogram set and reads a first-order estimate
 * of each of the 6 corrections off it, following the linearized measurement
 * model derived in the file header. Navigate between elements interactively
 * with the keyboard.
 *
 * @param fileListPath text file listing the .root files to read
 * @param drawDetIdx   index of the detector element to start on,
 *                     0 .. BmnDetectorModel::Alignables() - 1
 * @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 histo(const char* fileListPath = "/home/igor/bmnroot/align_new/cfg/run8/7466.noa_rlerr.list",
          Int_t drawDetIdx = 0,      // detector index to visualize, 0 .. Alignables() - 1
          Int_t drawMode = DM_MIX,   // slope panels: DM_PROFILE / DM_SCATTER / DM_MIX
          const char* pdfPath = nullptr)
{
    // init section
    BmnRun8cgTopDM detModel;
    detModel.Init();
    const Int_t nDet{detModel.Alignables()};
    if (drawDetIdx < 0 || drawDetIdx >= nDet) {
        std::cerr << " <!> histo: drawDetIdx " << drawDetIdx << " out of range [0, " << nDet << ")" << std::endl;
        return 1;
    }
    if (drawMode < DM_PROFILE || drawMode > DM_MIX) {
        std::cerr << " <!> histo: 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 << " <!> histo: 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
 
    std::vector<DetAcc> acc(nDet);     // per-detector X- and Y-residual accumulators
    std::vector<DetDraw> draw(nDet);   // per-detector per-hit records for drawing
 
    // walk tracks; fit a straight line per track, then read residuals
    std::cout << "Reading data...\n";
    Long64_t usedTracks{0}, usedHits{0};
    progressBar.Init();
    for (itr.ResetAll(); !itr.EndOfTracks(); itr.NextTrack()) {
        // weighted straight-line fit  x = alpha_x + beta_x * z,  y = alpha_y + beta_y * z
        BmnWeightedLine lineX, lineY;
        for (itr.ResetHits(); !itr.EndOfHits(); itr.NextHit()) {
            lineX.Add(itr.HitZ(), itr.HitX(), itr.HitWx());
            lineY.Add(itr.HitZ(), itr.HitY(), itr.HitWy());
        }
        if (!(lineX.Det() > 0.0) || !(lineY.Det() > 0.0))   // degenerate track (no z spread)
            continue;
 
        usedTracks++;
        Double_t betaX{lineX.Slope()}, alphaX{lineX.At0()}, betaY{lineY.Slope()}, alphaY{lineY.At0()};
        // residuals and covariates (bar = straight-line prediction)
        for (itr.ResetHits(); !itr.EndOfHits(); itr.NextHit(), usedHits++) {
            Int_t d{detModel.Idx(itr.HitDetectorID())};
            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[d].x[P_dX].Add(0.0, rx);               // mean  -> dX
            acc[d].x[P_dZ].Add(betaX, rx);             // slope -> -dZ
            acc[d].x[P_phi].Add(betaX * yBar, rx);     // slope -> -phi
            acc[d].x[P_theta].Add(betaX * xBar, rx);   // slope -> +theta
            acc[d].x[P_psi].Add(yBar, rx);             // slope -> -psi
            // Y-residual readings (alternative solution)
            acc[d].y[P_dY].Add(0.0, ry);               // mean  -> dY
            acc[d].y[P_dZ].Add(betaY, ry);             // slope -> -dZ
            acc[d].y[P_phi].Add(betaY * yBar, ry);     // slope -> -phi
            acc[d].y[P_theta].Add(betaY * xBar, ry);   // slope -> +theta
            acc[d].y[P_psi].Add(xBar, ry);             // slope -> +psi
 
            DetDraw& dd{draw[d]};
            dd.rx.push_back(rx);
            dd.ry.push_back(ry);
            dd.betaX.push_back(betaX);
            dd.betaXY.push_back(betaX * yBar);
            dd.betaXX.push_back(betaX * xBar);
            dd.yBar.push_back(yBar);
            dd.betaY.push_back(betaY);
            dd.betaYY.push_back(betaY * yBar);
            dd.betaYX.push_back(betaY * xBar);
            dd.xBar.push_back(xBar);
        }   // loop over hits
    }   // loop over tracks
    progressBar.Clear();
 
    if (usedTracks == 0) {
        std::cerr << " <!> histo: no usable tracks found" << std::endl;
        return 1;
    }
 
    // numeric report: extracted correction per detector
 
    using namespace std;
    cout << "\nHistogramming method: first-order corrections read from residuals\n";
    cout << usedTracks << " tracks, " << usedHits << " hits\n\n";
 
    // primary table: shifts from residual means, dZ/phi/theta/psi from X-residual
    cout << "Primary readings (shifts from means, dZ/phi/theta/psi from X-residual):\n";
    cout << "  " << left << setw(8) << "det_idx" << setw(10) << "det_id" << setw(22) << "element";
    cout << right << setw(11) << "samples";
    for (Int_t p = 0; p < PARAM_CNT; p++)
        cout << setw(13) << PARAM_NAMES[p];
    cout << setw(13) << "RMS(rx)" << setw(13) << "RMS(ry)" << "\n";
    for (Int_t d = 0; d < nDet; d++) {
        cout << "  " << left << setw(8) << d << setw(10) << detModel.EncodedID(d) << setw(22)
             << detModel.DecodeID(detModel.EncodedID(d));
        cout << right << setw(11) << acc[d].x[P_dX].line.N();
        cout << fixed << setprecision(6);
        cout << setw(13) << acc[d].x[P_dX].Value(P_dX, PARAM_SIGN_X) << setw(13)
             << acc[d].y[P_dY].Value(P_dY, PARAM_SIGN_Y);
        for (Int_t p = P_dZ; p < PARAM_CNT; p++)
            cout << setw(13) << acc[d].x[p].Value(p, PARAM_SIGN_X);
        cout << setw(13) << acc[d].x[P_dX].RmsY() << setw(13) << acc[d].y[P_dY].RmsY() << "\n";
    }
 
    // alternative table: dZ/phi/theta/psi from the Y-residual slopes
    cout << "\nAlternative readings (dZ/phi/theta/psi from Y-residual):\n";
    cout << "  " << left << setw(8) << "det_idx" << setw(10) << "det_id" << setw(22) << "element" << right;
    for (Int_t p = P_dZ; p < PARAM_CNT; p++)
        cout << setw(13) << PARAM_NAMES[p];
    cout << "\n";
    for (Int_t d = 0; d < nDet; d++) {
        cout << "  " << left << setw(8) << d << setw(10) << detModel.EncodedID(d) << setw(22)
             << detModel.DecodeID(detModel.EncodedID(d));
        cout << right << fixed << setprecision(6);
        for (Int_t p = P_dZ; p < PARAM_CNT; p++)
            cout << setw(13) << acc[d].y[p].Value(p, PARAM_SIGN_Y);
        cout << "\n";
    }
 
    // drawn objects are detached from any TDirectory and marked kCanDelete so
    // that a canvas 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", RES_X, RES_Y);
    pCnv->SetCanvasSize(RES_X - 20, RES_Y - 96);
 
    // (re)draw all diagnostics for detector idx and print its detailed readout
    auto show = [&](Int_t idx) {
        const DetAcc& da{acc[idx]};
        const DetDraw& dd{draw[idx]};
 
        // detailed readout: X- and Y-solutions side by side
        cout << "\nExtracted corrections for detector " << idx << " (id " << detModel.EncodedID(idx) << ", "
             << detModel.DecodeID(detModel.EncodedID(idx)) << ", " << da.x[P_dX].line.N() << " samples):\n";
        cout << fixed << setprecision(6);
        cout << "  " << left << setw(7) << "dX" << " = " << right << setw(11) << da.x[P_dX].Value(P_dX, PARAM_SIGN_X)
             << " +/- " << setw(10) << da.x[P_dX].Error(P_dX) << "   (" << PARAM_DIAG[P_dX] << ")\n";
        cout << "  " << left << setw(7) << "dY" << " = " << right << setw(11) << da.y[P_dY].Value(P_dY, PARAM_SIGN_Y)
             << " +/- " << setw(10) << da.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)
                 << da.x[p].Value(p, PARAM_SIGN_X) << " +/- " << setw(9) << da.x[p].Error(p) << "   Y: " << setw(11)
                 << da.y[p].Value(p, PARAM_SIGN_Y) << " +/- " << setw(9) << da.y[p].Error(p) << "\n";
        cout << flush;
 
        // fresh canvas: drop the previous picture and rebuild the pad grid
        pCnv->SetTitle(Form("BM@N histogramming : detector %d (id %d)", idx, detModel.EncodedID(idx)));
        pCnv->Clear();
        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();
            }
 
        // build the panel list: top row = shifts, middle row = X-residual
        // slopes, bottom row = Y-residual (alternative) slopes; one per column
        const Int_t slopeP[4]{P_dZ, P_phi, P_theta, P_psi};
        const VecDouble_t* covX[4]{&dd.betaX, &dd.betaXY, &dd.betaXX, &dd.yBar};
        const VecDouble_t* covY[4]{&dd.betaY, &dd.betaYY, &dd.betaYX, &dd.xBar};
        std::vector<Panel> panels;
        panels.push_back({0, 0, P_dX, kTRUE, kFALSE, kFALSE, &dd.rx, nullptr, RESIDUAL_RANGE_X, PARAM_DIAG[P_dX],
                          COV_TITLE_X[P_dX], &da.x[P_dX], PARAM_SIGN_X});
        panels.push_back({0, 1, P_dY, kTRUE, kTRUE, kFALSE, &dd.ry, nullptr, RESIDUAL_RANGE_Y, PARAM_DIAG[P_dY],
                          COV_TITLE_Y[P_dY], &da.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, &dd.rx, covX[c], RESIDUAL_RANGE_X,
                              PARAM_DIAG[p], COV_TITLE_X[p], &da.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, &dd.rx, covX[c], RESIDUAL_RANGE_X, PARAM_DIAG[p],
                                  COV_TITLE_X[p], &da.x[p], PARAM_SIGN_X});
            else
                panels.push_back({2, c, p, kFALSE, kTRUE, drawMode == DM_SCATTER, &dd.ry, covY[c], RESIDUAL_RANGE_Y,
                                  PARAM_DIAG_Y[p], COV_TITLE_Y[p], &da.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_%d", PARAM_NAMES[pn.p], idx), 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_%d", PARAM_NAMES[pn.p], suff, idx),
                                        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_%d", PARAM_NAMES[pn.p], suff, idx),
                                                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();
        }
        pCnv->Update();
    };
 
    show(drawDetIdx);
 
    if (pdfPath && pdfPath[0])
        pCnv->SaveAs(pdfPath);
 
    // interactive browser: Right/Left arrows step to the next/previous detector
    // (wrapping around), q or Esc quits. Skipped in batch mode or when stdin is
    // not a terminal (e.g. piped input) - there the single detector is drawn once.
    if (!gROOT->IsBatch() && isatty(STDIN_FILENO)) {
        cout << "\n[ Right / Left arrows: next / previous detector | q: quit ]\n" << flush;
        struct termios saved;
        histoSetRawMode(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')
                    drawDetIdx = (drawDetIdx + 1) % nDet;
                else if (a == '[' && b == 'D')
                    drawDetIdx = (drawDetIdx - 1 + nDet) % nDet;
                else
                    continue;
                show(drawDetIdx);
            }
        }
        histoSetRawMode(false, &saved);
    }
 
    return 0;
}