compareS.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 4.7 (Opus)
    Supervisor: Igor Polev, polev@jinr.int
 
    Visual comparison of two alignment solutions produced by
    BmnAligner::SaveSolution. Useful to compare an obtained solution
    against a known precise solution (e.g. the truth saved by
    BmnAligner::MakeTestTask).
*/
 
// headers below are for IDE syntax to work and they will be
// turned off during execution in ROOT macro mode
#ifndef __CLING__
#include "TAxis.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TPad.h"
 
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <vector>
// #include "TLegend.h"
#include "TLine.h"
#include "TString.h"
#endif
 
// nlohmann/json must be visible to CLING too. Its default header search path
// in ROOT macro mode may not include the locations used by BmnRoot.
#ifdef __CLING__
#pragma cling add_include_path("/home/igor/bmnroot/services")
#pragma cling add_include_path("/usr/include")
#endif
#include <nlohmann/json.hpp>
 
// number of correction parameters per detector (matches BMN_GLOBAL_PARAMS_PD)
#define PARAM_CNT 6
constexpr static const char* PARAM_NAMES[PARAM_CNT]{"dX",  "dY",    "dZ",
                                                    "phi", "theta", "psi"};   // matches BmnParamLayout::ParamNames()
constexpr static const char* PARAM_TITLES[PARAM_CNT]{"X shift, cm",       "Y shift, cm",         "Z shift, cm",
                                                     "phi rotation, rad", "theta rotation, rad", "psi rotation, rad"};
constexpr static const Double_t PARAM_RANGES[PARAM_CNT]{0.25, 0.25, 0.25, 0.02, 0.02, 0.02};
 
// canvas and pads configuration
#define RES_X 1920
#define RES_Y 1080
#define PADS_X 3
#define PADS_Y 2
 
struct CmpSolution
{
    std::vector<Int_t> ids;
    std::vector<std::vector<Double_t>> values;   // [detector][param]
    std::vector<std::vector<Double_t>> errors;   // [detector][param] - zeros if file omits "errors"
};
 
static Bool_t LoadCmpSolution(const char* path, CmpSolution& sol)
{
    std::ifstream ifs(path);
    if (!ifs.good()) {
        std::cerr << " <!> compareSolutions: failed to open " << path << std::endl;
        return kFALSE;
    }
    nlohmann::json js;
    try {
        ifs >> js;
    } catch (const std::exception& e) {
        std::cerr << " <!> compareSolutions: failed to parse JSON " << path << ":\n     " << e.what() << std::endl;
        return kFALSE;
    }
    if (!js.contains("CorrectionsPerDetector") || !js["CorrectionsPerDetector"].is_array()) {
        std::cerr << " <!> compareSolutions: 'CorrectionsPerDetector' missing or"
                     " not an array in "
                  << path << std::endl;
        return kFALSE;
    }
    const auto& arr = js["CorrectionsPerDetector"];
    sol.ids.clear();   // just in case
    sol.values.clear();
    sol.errors.clear();
    sol.ids.reserve(arr.size());
    sol.values.reserve(arr.size());
    sol.errors.reserve(arr.size());
    for (const auto& el : arr) {
        if (!el.contains("detector_id") || !el.contains("values") || !el["values"].is_array()
            || el["values"].size() != PARAM_CNT)
        {
            std::cerr << " <!> compareSolutions: malformed entry in " << path << std::endl;
            return kFALSE;
        }
        std::vector<Double_t> v;
        v.reserve(PARAM_CNT);
        for (const auto& x : el["values"])
            v.push_back(x.get<Double_t>());
        // "errors" is optional: fall back to zeros
        std::vector<Double_t> e(PARAM_CNT, 0.0);
        if (el.contains("errors") && el["errors"].is_array() && el["errors"].size() == PARAM_CNT) {
            for (Int_t k = 0; k < PARAM_CNT; k++)
                e[k] = el["errors"][k].get<Double_t>();
        }
        sol.ids.push_back(el["detector_id"].get<Int_t>());
        sol.values.push_back(std::move(v));
        sol.errors.push_back(std::move(e));
    }
    return kTRUE;
}
 
/**
 * @brief Visual comparison of two alignment solutions.
 *
 * Both files are read in the format BmnAligner::SaveSolution writes. Typical
 * use is checking an obtained solution against a known truth, such as the one
 * makett.C saves next to a synthetic task.
 *
 * @param obtainedPath .json solution under test
 * @param precisePath  .json solution to compare it against
 * @param pdfPath      .pdf file to save the canvas into, null to only display
 *
 * @return 0 on success, 1 if either solution could not be loaded
 */
int compareS(
    const char* obtainedPath = "/home/igor/bmnroot/align_new/cfg/solutions/sl9555.112e_lin3.04.json",
    const char* precisePath =
        "/home/igor/data.cache/tt9555.112e_lin3.04.solution.json",   //"/home/igor/bmnroot/align_new/cfg/solutions/Ruf_AZ_run9.json",
    const char* pdfPath = nullptr)
{
    // load both solutions
    CmpSolution sObt, sTru;
    if (!LoadCmpSolution(obtainedPath, sObt))
        return 1;
    if (!LoadCmpSolution(precisePath, sTru))
        return 1;
    if (sObt.ids.size() != sTru.ids.size()) {
        std::cerr << " <!> compareSolutions: detector counts differ (" << sObt.ids.size() << " obtained vs "
                  << sTru.ids.size() << " precise)" << std::endl;
        return 1;
    }
    const Int_t nDet = static_cast<Int_t>(sObt.ids.size());
 
    // align entries by detector_id (handles possible reordering)
    std::map<Int_t, Int_t> truIdxByID;
    for (Int_t i = 0; i < nDet; i++)
        truIdxByID[sTru.ids[i]] = i;
 
    std::vector<Int_t> ids(nDet);
    std::vector<std::vector<Double_t>> tVal(nDet, std::vector<Double_t>(PARAM_CNT, 0.0)),
        oVal(nDet, std::vector<Double_t>(PARAM_CNT, 0.0)), oErr(nDet, std::vector<Double_t>(PARAM_CNT, 0.0));
    for (Int_t i = 0; i < nDet; i++) {
        Int_t id = sObt.ids[i];
        auto it = truIdxByID.find(id);
        if (it == truIdxByID.end()) {
            std::cerr << " <!> compareSolutions: detector_id " << id << " not found in precise solution" << std::endl;
            return 1;
        }
        ids[i] = id;
        for (Int_t p = 0; p < PARAM_CNT; p++) {
            oVal[i][p] = sObt.values[i][p];
            oErr[i][p] = sObt.errors[i][p];
            tVal[i][p] = sTru.values[it->second][p];
        }
    }
 
    // numeric summary to stdout
    using namespace std;
    cout << "Per-parameter summary of (obtained - precise):\n";
    cout << "  " << left << setw(8) << "param"
         << right
         //<< setw(15) << "mean true"
         //<< setw(15) << "mean obt"
         //<< setw(15) << "mean diff"
         << setw(15) << "RMS diff" << setw(19) << "mean|diff/true|" << setw(18) << "max|diff/true|" << setw(10)
         << "det_idx" << setw(12) << "det_id"
         << "\n";
    constexpr Double_t TRUTH_EPS{1e-6};
    for (Int_t p = 0; p < PARAM_CNT; p++) {
        Double_t s{0.0}, s2{0.0}, t{0.0}, o{0.0}, sRel{0.0}, maxRel{-1.0};
        Int_t maxIdx{-1};
        for (Int_t i = 0; i < nDet; i++) {
            Double_t d{oVal[i][p] - tVal[i][p]};
            s += d;
            s2 += d * d;
            t += tVal[i][p];
            o += oVal[i][p];
            if (fabs(tVal[i][p]) < TRUTH_EPS)
                continue;   // skip near-zero truth
            Double_t r{fabs(d / tVal[i][p])};
            sRel += r;
            if (r > maxRel) {
                maxRel = r;
                maxIdx = i;
            }
        }
        cout << "  " << left << setw(8) << PARAM_NAMES[p] << right << fixed
             << setprecision(6)
             //<< setw(15) << t / nDet
             //<< setw(15) << o / nDet
             //<< setw(15) << s / nDet
             << setw(15) << sqrt(s2 / nDet) << setw(19) << sRel / nDet;
        if (maxIdx >= 0)
            cout << fixed << setprecision(1) << setw(18) << maxRel << setw(10) << maxIdx << setw(12) << ids[maxIdx];
        else
            cout << setw(18) << "n/a" << setw(10) << "-" << setw(12) << "-";
        cout << "\n";
    }
 
    // canvas + pads (RAM not released - acceptable in ROOT macro mode)
    TCanvas* pCnv = new TCanvas(
        "canvas", "BM@N alignment solutions comparison : BAR = true value, DOT = obtained (1-sigma conf), LINE = diff",
        RES_X, RES_Y);
    pCnv->SetCanvasSize(RES_X - 20, RES_Y - 96);
    constexpr 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();
        }
 
    // 6 panels: top row = shifts (dX, dY, dZ), bottom row = rotations (phi, theta, psi)
    for (Int_t p = 0; p < PARAM_CNT; p++) {
        Int_t col = p % PADS_X;
        Int_t row = p / PADS_X;
        // row 0 = shifts on top, row 1 = rotations below
        pPads[PADS_Y - 1 - row][col]->cd();
 
        std::vector<Double_t> x(nDet), xE(nDet, 0.0), yT(nDet), yO(nDet), yE(nDet), yD(nDet);
        Double_t yMin{+1e30}, yMax{-1e30};
        for (Int_t i = 0; i < nDet; i++) {
            x[i] = static_cast<Double_t>(i);
            yT[i] = tVal[i][p];
            yO[i] = oVal[i][p];
            yE[i] = oErr[i][p];
            yD[i] = oVal[i][p] - tVal[i][p];
            if (yT[i] < yMin)
                yMin = yT[i];
            if (yO[i] - yE[i] < yMin)
                yMin = yO[i] - yE[i];   // account for error bar
            if (yD[i] < yMin)
                yMin = yD[i];
            if (yT[i] > yMax)
                yMax = yT[i];
            if (yO[i] + yE[i] > yMax)
                yMax = yO[i] + yE[i];
            if (yD[i] > yMax)
                yMax = yD[i];
        }
        Double_t pad{0.1 * std::max(yMax - yMin, 1e-6)};
        if (yMax == yMin) {
            yMin -= 1.0e-3;
            yMax += 1.0e-3;
        }
 
        // per-series graphs
        // truth: teal bars
        TGraph* gT = new TGraph(nDet, x.data(), yT.data());
        gT->SetFillColor(kTeal);
        gT->SetLineColor(kTeal);
 
        // obtained: black markers + 1-sigma vertical error bars, no connecting line
        TGraphErrors* gO = new TGraphErrors(nDet, x.data(), yO.data(), xE.data(), yE.data());
        gO->SetMarkerStyle(kCircle);
        gO->SetMarkerColor(kBlack);
        gO->SetLineColor(kRed);
 
        // diff: red line, no markers
        TGraph* gD = new TGraph(nDet, x.data(), yD.data());
        gD->SetLineColor(kRed);
        gD->SetLineWidth(2);
 
        // bundle into a multigraph - acts as the axis frame
        TMultiGraph* mg = new TMultiGraph();
        mg->SetName(Form("mg_%s", PARAM_NAMES[p]));
        mg->SetTitle(Form("%s;detector index;", PARAM_TITLES[p]));
        mg->Add(gT, "B");   // bars first (drawn at back)
        mg->Add(gD, "L");   // diff line on top of bars
        mg->Add(gO, "P");   // obtained markers on top of everything
        mg->Draw("A");
        mg->GetXaxis()->SetLimits(-0.5, nDet - 0.5);
        mg->GetXaxis()->SetNdivisions(nDet);
        if (PARAM_RANGES[p] > 0.0)
            mg->GetYaxis()->SetRangeUser(-PARAM_RANGES[p], +PARAM_RANGES[p]);
        else
            mg->GetYaxis()->SetRangeUser(yMin - pad, yMax + pad);
 
        // zero reference line (drawn after multigraph so pad coords are set)
        TLine* zero = new TLine(-0.5, 0.0, nDet - 0.5, 0.0);
        zero->SetLineColor(kGray + 1);
        zero->SetLineStyle(3);
        zero->Draw();
 
        // legends removed to avoid overlaping with data
        // moved to window title in some way
        /*
        TLegend* leg = new TLegend(0.65, 0.74, 0.95, 0.92);
        leg->SetTextSize(0.035);
        leg->AddEntry(gT, "Precise (truth)",  "f");
        leg->AddEntry(gO, "Obtained",         "p");
        leg->AddEntry(gD, "Diff (obt - tru)", "l");
        leg->Draw();
        */
    }
 
    if (pdfPath && pdfPath[0])
        pCnv->SaveAs(pdfPath);
 
    return 0;
}