BmnAlignConfig.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:     Claude (Opus 5)
    Supervisor: Igor Polev, polev@jinr.int
 
    BmnAlignConfig class implementation
*/
 
#include "BmnAlignConfig.h"
 
#include "fairlogger/Logger.h"
 
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <nlohmann/json.hpp>
#include <regex>
#include <string>
 
namespace
{
// A constraint RHS this small is meant as zero; see the comment at the use site.
constexpr Double_t SOFT_ZERO{1.0E-6};
 
// objects a config entry addresses; they point into the layout of the
// hierarchy being loaded, which does not change while a file is read
typedef std::vector<const NodeRef*> NodeList_t;
 
/*
 *  Locators shared by the constraints and the presigmas files. Each returns
 *  an empty string on success, or what was wrong with the value, which the
 *  caller logs together with the number of the entry.
 */
 
// "ParamName": one name or an array of names -> parameter indices
std::string ParamList(const nlohmann::json& js, VecInt_t& params)
{
    params.clear();
    // CAUTION: no braces here - "json list {x}" makes the array [x], not a copy of x
    const nlohmann::json list = js.is_array() ? js : nlohmann::json::array({js});
    for (const nlohmann::json& name : list) {
        const Int_t p{name.is_string() ? BmnParamLayout::ParamIndex(name.get<std::string>()) : -1};
        if (p < 0)
            return "unknown ParamName " + name.dump();
        params.push_back(p);
    }
    return params.empty() ? "empty ParamName array" : "";
}
 
// "Over": a set of detector elements, never composites
//   Elements - all of them, the target included
//   Tracker  - all but the target (all of them in a model without one)
//   Target   - the target only (nothing in a model without one)
std::string LeafScope(const nlohmann::json& js, const BmnAlignHierarchy& hierarchy, NodeList_t& nodes)
{
    nodes.clear();
    const std::string over{js.is_string() ? js.get<std::string>() : js.dump()};
    if (over != "Elements" && over != "Tracker" && over != "Target")
        return "unknown Over value '" + over + "', expected Elements / Tracker / Target";
    const BmnDetectorModel* model{hierarchy.DetModel()};
    const Int_t target{model->TargetID() < 0 ? -1 : model->Idx(model->TargetID())};
    for (const NodeRef& node : hierarchy.Layout().ByNode()) {
        if (node.IsComposite())
            continue;
        const Bool_t isTarget{node.leaf == target};
        if (over == "Elements" || (over == "Tracker" && !isTarget) || (over == "Target" && isTarget))
            nodes.push_back(&node);
    }
    return "";
}
 
// "NodeMatch": the composites whose whole name matches a regular expression
std::string MatchComposites(const nlohmann::json& js, const BmnAlignHierarchy& hierarchy, NodeList_t& nodes)
{
    nodes.clear();
    if (!js.is_string())
        return "NodeMatch " + js.dump() + " is not a string";
    std::regex re;
    try {
        re.assign(js.get<std::string>());
    } catch (const std::regex_error& e) {
        return "bad NodeMatch " + js.dump() + ": " + e.what();
    }
    for (const NodeRef& node : hierarchy.Layout().ByNode())
        if (node.IsComposite() && std::regex_match(node.name, re))
            nodes.push_back(&node);
    return "";
}
}   // anonymous namespace
 
Bool_t BmnAlignConfig::ReadSolution(const char* jsonPath, const BmnAlignHierarchy& hierarchy, Solution& solution)
{
    solution = Solution{};
    solution.zeroComposites = hierarchy.CompositeCount() > 0;
 
    std::ifstream ifs(jsonPath);
    if (!ifs.good()) {
        LOG(error) << " <!> BmnAlignConfig::ReadSolution: failed to open JSON file";
        return kFALSE;
    }
    using namespace nlohmann;
    json js;
    try {
        ifs >> js;
 
        // check if solution is suitable
        json jsAlignmentType = BmnParamLayout::ParamNames();
        if (js.contains("CorrectionValues")) {
            if (js["CorrectionValues"] != jsAlignmentType) {
                LOG(error) << " <!> ReadSolution: wrong alignment type in JSON file\n"
                           << "     expected: " << jsAlignmentType << ", in file: " << js["CorrectionValues"];
                return kFALSE;
            } else
                LOG(info) << "    " << js["CorrectionValues"] << " : alignment type checked";
        } else
            LOG(warning) << "  ! CorrectionValues key is missing: assuming alignment type is correct";
        const Int_t elements{hierarchy.LeafCount()};
        if (js.contains("DetectorElements")) {
            if (static_cast<Int_t>(js["DetectorElements"]) != elements) {
                LOG(error) << " <!> ReadSolution: wrong detector configuration in JSON file\n"
                           << "     modules expected: " << elements << ", in file: " << js["DetectorElements"];
                return kFALSE;
            } else
                LOG(info) << "    " << js["DetectorElements"] << " : detector elements count checked";
        } else
            LOG(warning) << "  ! DetectorElements key is missing: assuming detector configuration is correct";
        if (!js.contains("CorrectionsPerDetector") || !js["CorrectionsPerDetector"].is_array()) {
            LOG(error) << "   ReadSolution: CorrectionsPerDetector key is missing or not an array";
            return kFALSE;
        }
 
        // read solution for each stored detector element
        Int_t elementsCount{static_cast<Int_t>(js["CorrectionsPerDetector"].size())};
        LOG_IF(warning, elementsCount != elements)
            << Form(" ! %d detector elements found in file, %d was expected", elementsCount, elements);
        LOG_IF(warning, elementsCount > elements) << " ! additional data validation is advised";
        LOG_IF(warning, elementsCount < elements) << " ! assuming zero values for missing detector elements";
        for (Int_t i = 0; i < elementsCount; i++) {
            json& element = js["CorrectionsPerDetector"][i];
            if (!element.contains("detector_id")) {
                LOG(warning) << "  ! DetectorID key is missing for some element, skipping.";
                continue;
            }
            Int_t detID = element["detector_id"].get<Int_t>();
            if (hierarchy.DetModel()->UnknownID(detID)) {
                LOG(warning) << "  ! Detector element ID " << detID << " (" << hierarchy.DetModel()->DecodeID(detID)
                             << ")"
                             << " is unknown or not used in current configuration, skipping";
                continue;
            }
            if (!element.contains("values") || !element["values"].is_array()
                || element["values"].size() != BMN_GLOBAL_PARAMS_PD)
            {
                LOG(warning) << "  ! Correction values for the element " << detID << " ("
                             << hierarchy.DetModel()->DecodeID(detID) << ")"
                             << " are invalid, skipping.";
                continue;
            }
            const Int_t detIdx{hierarchy.DetModel()->Idx(detID)};
            NodeValues leaf;
            leaf.node = detIdx;
            json& values = element["values"];
            for (Int_t j = 0; j < BMN_GLOBAL_PARAMS_PD; j++)
                leaf.values(j) = values[j].get<Double_t>();
            if (detID == hierarchy.DetModel()->TargetID())   // see the target note in the header
                leaf.values(idZ) += hierarchy.DetModel()->TargetShift();
            solution.leaves.push_back(leaf);
            // "errors" key (if present in the file) is intentionally ignored
        }
        LOG(info) << "    " << solution.leaves.size() << " detector elements loaded from JSON file";
 
        // Composite node parameters, see WriteSolution. Matched by name, so a
        // file written for a different hierarchy is reported node by node
        // instead of being misapplied by position. An absent or partial block
        // leaves the corresponding parameters at zero: the geometry loaded above
        // is unaffected, the whole of it just counts as the elements' own - which
        // is exactly what a file written before the levels existed means.
        const Int_t nCompNodes{hierarchy.CompositeCount()};
        if (nCompNodes) {
            if (js.contains("CompositeNodes") && js["CompositeNodes"].is_array())
                for (const json& entry : js["CompositeNodes"]) {
                    if (!entry.contains("name") || !entry.contains("values") || !entry["values"].is_array()
                        || entry["values"].size() != BMN_GLOBAL_PARAMS_PD)
                    {
                        LOG(warning) << "  ! malformed composite node entry, skipping";
                        continue;
                    }
                    const std::string name{entry["name"].get<std::string>()};
                    const NodeRef* node{hierarchy.Layout().Find(name)};
                    if (!node || !node->IsComposite()) {
                        LOG(warning) << "  ! '" << name
                                     << "' is not a composite node of"
                                        " the current hierarchy, skipping";
                        continue;
                    }
                    NodeValues comp;
                    comp.node = node->params[0].node;
                    for (const ParamRef& p : node->params)
                        comp.values(p.prm) = entry["values"][p.prm].get<Double_t>();
                    solution.composites.push_back(comp);
                }
            LOG(info) << "    " << solution.composites.size() << " of " << nCompNodes
                      << " composite node(s) loaded from JSON file";
            LOG_IF(warning, (Int_t)solution.composites.size() < nCompNodes)
                << "  ! the rest start at zero, so their share of the corrections"
                   " counts as the elements' own until the next iteration";
        }
 
        // TODO: load tracks approximation form JSON file if available
    } catch (std::exception& e) {
        LOG(error) << " <!> ReadSolution: failed to parse JSON file " << jsonPath;
        LOG(error) << e.what();
        return kFALSE;
    }
    return kTRUE;
}
 
Bool_t BmnAlignConfig::LoadConstraints(const char* jsonPath, const BmnAlignHierarchy& hierarchy, BmnAlignPriors& priors)
{
    /*
     * JSON schema for constraints file
     * --------------------------------
     * Top-level value: array of entries. An entry is either one constraint
     * written out coefficient by coefficient (constraint object), or a typed
     * object the loader turns into a group of constraints.
     *
     * Detector elements are named by DetectorID, never by their index: an
     * index means a different element in every model (full, half, without
     * the target), and nothing could tell a file written for another model.
     *
     * Constraint object:
     *   - "Coefficients" (array, required): list of coefficient objects.
     *   - any other string-valued attribute (e.g. "Comment"): ignored. Used for
     *     comments. Order of attributes is irrelevant.
     *   - "IDComment" (optional, any value): ignored. service/comment_ids.C
     *     writes the text form of the DetectorIDs there, in any object type.
     *
     * Coefficient object (one of two mutually exclusive schemes per constraint):
     *
     *   Scheme A (locator) — must use *exactly one* of {ParamName, ParamIndex}
     *     and *exactly one* of {DetectorID, NodeName}:
     *       { "ParamName": "dX",   "DetectorID": 201392146, "Value": 1.0 }
     *       { "ParamIndex": 0,     "DetectorID": 201392146, "Value": 1.0 }
     *       { "ParamName": "theta","NodeName":   "GEM",     "Value": 1.0 }
     *
     *   NodeName addresses a composite node of the hierarchy (a leaf may also be
     *   named this way). Note the difference in meaning: a coefficient on a leaf
     *   parameter constrains that element's *effective* correction, the sum of
     *   its own parameters and everything its ancestors contribute, while a
     *   coefficient on a composite constrains that composite's own parameters.
     *   Leaf coefficients therefore keep working unchanged when a hierarchy is
     *   introduced - they still fix the same physical degree of freedom.
     *
     *   Scheme B (direct):
     *       { "ConstrIndex": 17, "Value": 1.0 }   // 0..AllParams()-1
     *
     *   RHS:                                          // at most one per constraint
     *       { "RHS": 0.0 }                            // no other attributes allowed
     *
     *   Mixing schemes within one constraint is an error.
     *   Missing coefficients are zero (no warning).
     *
     * Typed object: the loader builds the constraints from the hierarchy and
     * the detector model, so one entry stays right for every model it is used
     * with - full or half, flat or hierarchical, with or without the target.
     *   - "Type" (required): Fix | Sum | Shear | NodeEqual, see below.
     *   - "ParamName" (required): one name or an array of names. The entry is
     *     applied to each of them in turn.
     *   - exactly one of:
     *       "Over":      "Elements" | "Tracker" | "Target" - detector elements:
     *                    all, all but the target, the target only
     *       "NodeMatch": regular expression; selects the composites whose
     *                    whole name matches it, e.g. "STA0[6-9][UL]"
     *       "ID":        array of DetectorIDs, NodeEqual only
     *   - "RHS" (optional, not for NodeEqual): the right-hand side of every row
     *     made, 0.0 if absent.
     *   - "Comment", "IDComment" (optional). Nothing else is allowed: in a
     *     typed object a string value is not a comment, so a misspelt key
     *     cannot vanish.
     *
     *   Type       rows made                     each row says
     *   Fix        one per object                P = RHS
     *   Sum        one                           sum of P over the objects = RHS
     *   Shear      one, Over only                sum of w_i P_i = RHS
     *   NodeEqual  n-1 per group of n elements   P of an element = P of the first one
     *
     *   Shear weights are w_i = (z_i - z_mean) / max |z_i - z_mean|, with z of
     *   each element from the detector model geometry. Straight tracks do not
     *   see a dX or dY growing linearly with z: it goes into their slopes. Lin3
     *   sees dZ only as -dZ * slope, so a dZ growing linearly with z (a scale)
     *   goes the same way.
     *   A NodeEqual group is the elements anywhere below one matched composite,
     *   or the ID list. A group of a single element makes no rows and is
     *   reported.
     *
     *   As for coefficients, rows on elements constrain their effective
     *   corrections, rows on composites their own parameters.
     *   An entry that addresses nothing in the current model makes no rows,
     *   only a warning, so one file can serve several models. Each typed entry
     *   logs the numbers of the rows it made: later messages (BmnAlignPriors,
     *   BmnAlignReport) count rows, not entries.
     */
    if (!jsonPath) {
        LOG(warning) << "  ! constraints file in not provided, proceeding without constraints";
        return kTRUE;
    }
    LOG(info) << "    " << jsonPath << " : constraints file provided";
    std::ifstream ifs(jsonPath);
    if (!ifs.good()) {
        LOG(error) << " <!> BmnAlignConfig::LoadConstraints: failed to open JSON file " << jsonPath;
        return kFALSE;
    }
    nlohmann::json jsonData;
    try {
        ifs >> jsonData;
    } catch (const std::exception& e) {
        LOG(error) << " <!> BmnAlignConfig::LoadConstraints: failed to parse JSON file " << jsonPath;
        LOG(error) << e.what();
        return kFALSE;
    }
    if (!jsonData.is_array()) {
        LOG(error) << " <!> BmnAlignConfig::LoadConstraints: top-level JSON is not an array in " << jsonPath;
        return kFALSE;
    }
    const Int_t entriesCnt{static_cast<Int_t>(jsonData.size())};
    if (entriesCnt <= 0) {
        LOG(error) << " <!> BmnAlignConfig::LoadConstraints: empty constraints array in " << jsonPath;
        return kFALSE;
    }
 
    auto detIdxByID = [&hierarchy](Int_t encID) -> Int_t {
        try {
            return hierarchy.DetModel()->Idx(encID);
        } catch (const std::out_of_range&) {
            return -1;
        }
    };
    auto isCoefKey = [](const std::string& k) -> Bool_t {
        return k == "ParamName" || k == "ParamIndex" || k == "DetectorID" || k == "NodeName" || k == "ConstrIndex"
               || k == "Value" || k == "RHS";
    };
    auto isTypedKey = [](const std::string& k) -> Bool_t {
        return k == "Type" || k == "ParamName" || k == "Over" || k == "NodeMatch" || k == "ID" || k == "RHS"
               || k == "Comment" || k == "IDComment";
    };
 
    // a typed entry makes any number of rows, so they are collected first
    const Int_t nFull{hierarchy.AllParams()};
    std::vector<TVectorD> rows;
    VecDouble_t rowsRhs;
    Int_t errors{0};
    auto addRow = [&rows, &rowsRhs](const TVectorD& row, Double_t rhs) {
        rows.push_back(row);
        rowsRhs.push_back(rhs);
    };
    priors.inhomogeneous = kFALSE;
    auto noteRhs = [&priors](Double_t rhs, Int_t entry) {
        if (std::fabs(rhs) <= SOFT_ZERO)
            return;
        priors.inhomogeneous = kTRUE;
        LOG(warning) << "  ! constraint #" << entry
                     << " has non-zero RHS -- (unlikely) violation is possible"
                        " (see comment in BmnAligner::SolveGlobal)";
    };
 
    for (Int_t i = 0; i < entriesCnt; i++) {
        const auto& jsConstr = jsonData[i];
        if (!jsConstr.is_object()) {
            LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << " is not an object";
            errors++;
            continue;
        }
 
        // ==== typed object ====
        if (jsConstr.contains("Type")) {
            Int_t entryErrors{0};
            auto fail = [&](const std::string& what) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << " : " << what;
                entryErrors++;
            };
            for (auto it = jsConstr.begin(); it != jsConstr.end(); ++it)
                if (!isTypedKey(it.key()))
                    fail("unexpected attribute '" + it.key() + "'");
            const nlohmann::json& jsType{jsConstr.at("Type")};
            const std::string type{jsType.is_string() ? jsType.get<std::string>() : ""};
            const Bool_t fix{type == "Fix"}, sum{type == "Sum"}, shear{type == "Shear"}, equal{type == "NodeEqual"};
            if (!fix && !sum && !shear && !equal)
                fail("unknown Type " + jsType.dump() + ", expected Fix / Sum / Shear / NodeEqual");
            VecInt_t params;
            if (!jsConstr.contains("ParamName"))
                fail("ParamName is required");
            else if (const std::string what{ParamList(jsConstr.at("ParamName"), params)}; !what.empty())
                fail(what);
            const Bool_t hasOver{jsConstr.contains("Over")}, hasMatch{jsConstr.contains("NodeMatch")},
                hasID{jsConstr.contains("ID")}, hasRHS{jsConstr.contains("RHS")};
            if (hasOver + hasMatch + hasID != 1)
                fail("exactly one of Over / NodeMatch / ID required");
            if (hasID && !equal)
                fail("ID is for NodeEqual only");
            if (hasOver && equal)
                fail("NodeEqual takes NodeMatch or ID, not Over");
            if (hasMatch && shear)
                fail("Shear weighs elements by their z, so it takes Over, not NodeMatch");
            if (hasRHS && equal)
                fail("NodeEqual takes no RHS");
            if (hasRHS && !jsConstr.at("RHS").is_number())
                fail("RHS " + jsConstr.at("RHS").dump() + " is not a number");
            if (entryErrors) {
                errors += entryErrors;
                continue;
            }
            const Double_t rhs{hasRHS ? jsConstr.at("RHS").get<Double_t>() : 0.0};
 
            // the addressed objects; only NodeEqual has more than one group of them
            std::vector<NodeList_t> groups(1);
            std::string what;
            if (hasOver)
                what = LeafScope(jsConstr.at("Over"), hierarchy, groups[0]);
            else if (hasMatch)
                what = MatchComposites(jsConstr.at("NodeMatch"), hierarchy, groups[0]);
            else if (const nlohmann::json & ids{jsConstr.at("ID")}; !ids.is_array())
                what = "ID " + ids.dump() + " is not an array";
            else
                for (const nlohmann::json& id : ids) {
                    const Int_t detIdx{id.is_number_integer() ? detIdxByID(id.get<Int_t>()) : -1};
                    if (detIdx < 0) {
                        what = "unknown DetectorID " + id.dump();
                        if (id.is_number_integer())
                            what += " (" + hierarchy.DetModel()->DecodeID(id.get<Int_t>()) + ")";
                        break;
                    }
                    if (std::find(groups[0].begin(), groups[0].end(), &hierarchy.Layout().ByNode()[detIdx])
                        != groups[0].end())
                    {
                        what = "ID lists " + id.dump() + " twice";
                        break;
                    }
                    groups[0].push_back(&hierarchy.Layout().ByNode()[detIdx]);
                }
            if (what.empty() && shear
                && hierarchy.DetModel()->Geometry().size() != static_cast<std::size_t>(hierarchy.LeafCount()))
                what = "Shear needs the element z, but the detector model has no geometry";
            if (!what.empty()) {
                fail(what);
                errors += entryErrors;
                continue;
            }
 
            if (equal && hasMatch) {   // one group per matched composite: the elements anywhere below it
                const NodeList_t composites{groups[0]};
                groups.clear();
                for (const NodeRef* comp : composites) {
                    NodeList_t group;
                    for (const NodeRef& node : hierarchy.Layout().ByNode()) {
                        if (node.IsComposite())
                            continue;
                        for (Int_t n = hierarchy.Nodes()[node.leaf].parent; n >= 0; n = hierarchy.Nodes()[n].parent)
                            if (hierarchy.Nodes()[n].name == comp->name) {
                                group.push_back(&node);
                                break;
                            }
                    }
                    LOG_IF(warning, group.size() < 2) << "  ! constraint #" << i + 1 << " : composite '" << comp->name
                                                      << "' holds a single element, skipped";
                    if (group.size() >= 2)
                        groups.push_back(group);
                }
            } else if (equal && groups[0].size() < 2) {
                LOG(warning) << "  ! constraint #" << i + 1 << " : the ID list names a single element, skipped";
                groups.clear();
            }
 
            VecDouble_t weight;   // Shear only, in the order of groups[0]
            if (shear && !groups[0].empty()) {
                const auto& geometry{hierarchy.DetModel()->Geometry()};
                Double_t zMean{0.0}, spread{0.0};
                for (const NodeRef* node : groups[0])
                    zMean += geometry[node->leaf].z;
                zMean /= groups[0].size();
                for (const NodeRef* node : groups[0])
                    spread = std::max(spread, std::fabs(geometry[node->leaf].z - zMean));
                if (!(spread > 0.0)) {
                    fail("Shear over elements which all lie at one z");
                    errors += entryErrors;
                    continue;
                }
                for (const NodeRef* node : groups[0])
                    weight.push_back((geometry[node->leaf].z - zMean) / spread);
            }
 
            const std::size_t firstRow{rows.size()};
            for (const Int_t p : params) {
                if (fix)
                    for (const NodeRef* node : groups[0]) {
                        TVectorD row(nFull);
                        row(node->params[p].full) = 1.0;
                        addRow(row, rhs);
                    }
                if ((sum || shear) && !groups[0].empty()) {
                    TVectorD row(nFull);
                    for (std::size_t k = 0; k < groups[0].size(); k++)
                        row(groups[0][k]->params[p].full) = shear ? weight[k] : 1.0;
                    addRow(row, rhs);
                }
                if (equal)
                    for (const NodeList_t& group : groups)
                        for (std::size_t k = 1; k < group.size(); k++) {
                            TVectorD row(nFull);
                            row(group[0]->params[p].full) = -1.0;
                            row(group[k]->params[p].full) = 1.0;
                            addRow(row, 0.0);
                        }
            }
            if (rows.size() == firstRow) {
                LOG(warning) << "  ! constraint #" << i + 1 << " (" << type
                             << ") addresses nothing in this model, no rows made";
                continue;
            }
            noteRhs(rhs, i + 1);
            LOG(info) << Form("    constraint #%d (%s) : %zu row(s), rows #%zu..#%zu", i + 1, type.c_str(),
                              rows.size() - firstRow, firstRow + 1, rows.size());
            continue;
        }
 
        // ==== constraint object ====
        const nlohmann::json* jsCoeffs{nullptr};
        for (auto it = jsConstr.begin(); it != jsConstr.end(); ++it) {
            if (it.value().is_string() || it.key() == "IDComment")
                continue;   // comment-like attribute, ignored
            if (it.key() == "Coefficients" && it.value().is_array())
                jsCoeffs = &it.value();
            else {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1
                           << " has unexpected attribute '" << it.key() << "'";
                errors++;
            }
        }
        if (!jsCoeffs) {
            LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1
                       << " has no 'Coefficients' array";
            errors++;
            continue;
        }
 
        enum Scheme
        {
            Unknown,
            Locator,
            Direct
        };
        Scheme scheme{Unknown};
        Bool_t hasRHS{kFALSE};
        Bool_t anyLhs{kFALSE};
        TVectorD row(nFull);
        Double_t rowRhs{0.0};
 
        for (Int_t cIdx = 0; cIdx < (Int_t)jsCoeffs->size(); cIdx++) {
            const auto& jsCoef = (*jsCoeffs)[cIdx];
            if (!jsCoef.is_object()) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                           << cIdx + 1 << " is not an object";
                errors++;
                continue;
            }
            // Reject unknown keys inside a coefficient (typo guard).
            Bool_t badKey{kFALSE};
            for (auto it = jsCoef.begin(); it != jsCoef.end(); ++it)
                if (!isCoefKey(it.key())) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " has unknown attribute '" << it.key() << "'";
                    badKey = kTRUE;
                }
            if (badKey) {
                errors++;
                continue;
            }
 
            const Bool_t hasParamName{jsCoef.contains("ParamName")};
            const Bool_t hasParamIndex{jsCoef.contains("ParamIndex")};
            const Bool_t hasDetectorID{jsCoef.contains("DetectorID")};
            const Bool_t hasNodeName{jsCoef.contains("NodeName")};
            const Bool_t hasConstrIndex{jsCoef.contains("ConstrIndex")};
            const Bool_t hasValue{jsCoef.contains("Value")};
            const Bool_t hasRHSKey{jsCoef.contains("RHS")};
 
            // ---- RHS coefficient ----
            if (hasRHSKey) {
                if (hasParamName || hasParamIndex || hasDetectorID || hasNodeName || hasConstrIndex || hasValue) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " has RHS together with other attributes";
                    errors++;
                    continue;
                }
                if (hasRHS) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1
                               << " has multiple RHS entries";
                    errors++;
                    continue;
                }
                rowRhs = jsCoef["RHS"].get<Double_t>();
                noteRhs(rowRhs, i + 1);
                hasRHS = kTRUE;
                continue;
            }
 
            // ---- LHS coefficient ----
            if (!hasValue) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                           << cIdx + 1 << " has neither Value nor RHS";
                errors++;
                continue;
            }
            const Double_t value{jsCoef["Value"].get<Double_t>()};
 
            if (hasConstrIndex) {
                if (hasParamName || hasParamIndex || hasDetectorID || hasNodeName) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " mixes ConstrIndex with locator attributes";
                    errors++;
                    continue;
                }
                if (scheme == Locator) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1
                               << " mixes locator and direct schemes";
                    errors++;
                    continue;
                }
                scheme = Direct;
                const Int_t idx{jsCoef["ConstrIndex"].get<Int_t>()};
                if (idx < 0 || idx >= hierarchy.AllParams()) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " : ConstrIndex " << idx << " out of range [0, " << hierarchy.AllParams()
                               << ")";
                    errors++;
                    continue;
                }
                row(idx) = value;
                anyLhs = kTRUE;
                continue;
            }
 
            // locator scheme
            if (hasParamName == hasParamIndex) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                           << cIdx + 1 << " : exactly one of ParamName / ParamIndex required";
                errors++;
                continue;
            }
            if (hasDetectorID + hasNodeName != 1) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                           << cIdx + 1 << " : exactly one of DetectorID / NodeName required";
                errors++;
                continue;
            }
            if (scheme == Direct) {
                LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1
                           << " mixes locator and direct schemes";
                errors++;
                continue;
            }
            scheme = Locator;
 
            Int_t paramIdx{-1};
            if (hasParamName) {
                const std::string name{jsCoef["ParamName"].get<std::string>()};
                paramIdx = BmnParamLayout::ParamIndex(name);
                if (paramIdx < 0) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " : unknown ParamName '" << name << "'";
                    errors++;
                    continue;
                }
            } else {
                paramIdx = jsCoef["ParamIndex"].get<Int_t>();
                if (paramIdx < 0 || paramIdx >= BMN_GLOBAL_PARAMS_PD) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " : ParamIndex " << paramIdx << " out of range [0, "
                               << BMN_GLOBAL_PARAMS_PD << ")";
                    errors++;
                    continue;
                }
            }
 
            // the addressed object, however the entry chose to name it
            const NodeRef* target{nullptr};
            if (hasNodeName) {
                const std::string name{jsCoef["NodeName"].get<std::string>()};
                target = hierarchy.Layout().Find(name);
                if (!target) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " : unknown NodeName '" << name << "'";
                    errors++;
                    continue;
                }
            } else {
                const Int_t encID{jsCoef["DetectorID"].get<Int_t>()};
                const Int_t detIdx{detIdxByID(encID)};
                if (detIdx < 0) {
                    LOG(error) << " <!> BmnAlignConfig::LoadConstraints: constraint #" << i + 1 << ", coefficient #"
                               << cIdx + 1 << " : unknown DetectorID " << encID << " ("
                               << hierarchy.DetModel()->DecodeID(encID) << ")";
                    errors++;
                    continue;
                }
                target = &hierarchy.Layout().ByNode()[detIdx];
            }
            row(target->params[paramIdx].full) = value;
            anyLhs = kTRUE;
        }
        if (!anyLhs) {
            LOG(warning) << "  ! constraint #" << i + 1 << " has all-zero LHS - it is either trivial or unsatisfiable";
        }
        addRow(row, rowRhs);
    }
 
    if (errors) {
        LOG(error) << " <!> BmnAlignConfig::LoadConstraints: " << errors << " error(s) detected, aborting";
        return kFALSE;
    }
    const Int_t rowsCnt{static_cast<Int_t>(rows.size())};
    priors.constraints.ResizeTo(rowsCnt, nFull);
    priors.rhs.ResizeTo(rowsCnt);
    for (Int_t r = 0; r < rowsCnt; r++) {
        TMatrixDRow(priors.constraints, r) = rows[r];
        priors.rhs(r) = rowsRhs[r];
    }
    LOG(info) << "    " << rowsCnt << " constraint row(s) loaded from " << entriesCnt << " entries";
    return kTRUE;
}
 
Bool_t BmnAlignConfig::LoadPresigmas(const char* jsonPath, const BmnAlignHierarchy& hierarchy, BmnAlignPriors& priors)
{
    /*
     * JSON schema for presigmas file
     * ------------------------------
     * Top-level value: array of presigma objects. Each object attaches a prior
     * uncertainty to one global parameter (of one detector element, or of all
     * of them at once), turning it into a pseudo-measurement of that parameter
     * - see BmnAligner::ApplyPresigmas for the mathematics.
     *
     * Presigma object:
     *   - exactly one of {ParamName, ParamIndex} (required): which correction.
     *       "ParamName":  "dX" | "dY" | "dZ" | "phi" | "theta" | "psi", or an
     *                     array of them for the same prior on each
     *       "ParamIndex": 0..BMN_GLOBAL_PARAMS_PD-1
     *   - at most one of {DetectorID, NodeName, Over, NodeMatch} (optional):
     *       which objects. Elements are named by DetectorID, never by index
     *       (see LoadConstraints).
     *       "NodeName"  addresses one hierarchy node (a composite, or a leaf).
     *       "Over"      "Elements" | "Tracker" | "Target": detector elements -
     *                   all, all but the target, the target only.
     *       "NodeMatch" regular expression: the composites whose whole name
     *                   matches it.
     *       If all are absent, the entry applies to every detector element,
     *       like "Over": "Elements". Composite nodes are reached by NodeName
     *       or NodeMatch only, so that adding a level cannot silently shift
     *       the balance between levels. An entry that addresses nothing in
     *       the current model is reported and has no effect.
     *   - "Sigma" (required, > 0): prior uncertainty, in the units of the
     *       parameter (cm for dX/dY/dZ, rad for phi/theta/psi).
     *   - "Ref" (optional, default 0.0): value the prior pulls towards, in the
     *       same units. Zero means the nominal (starting) geometry.
     *   - any other string-valued attribute (e.g. "Comment"): ignored. Used for
     *     comments. Order of attributes is irrelevant.
     *   - "IDComment" (optional, any value): ignored, see LoadConstraints.
     *
     * Later entries override earlier ones for the same parameter, so a
     * detector-wide default may be stated first and then refined per element.
     * Parameters never mentioned get no prior at all.
     */
    if (!jsonPath) {
        LOG(info) << "    presigmas file is not provided, proceeding without priors";
        return kTRUE;
    }
    LOG(info) << "    " << jsonPath << " : presigmas file provided";
    std::ifstream ifs(jsonPath);
    if (!ifs.good()) {
        LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: failed to open JSON file " << jsonPath;
        return kFALSE;
    }
    nlohmann::json jsonData;
    try {
        ifs >> jsonData;
    } catch (const std::exception& e) {
        LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: failed to parse JSON file " << jsonPath;
        LOG(error) << e.what();
        return kFALSE;
    }
    if (!jsonData.is_array()) {
        LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: top-level JSON is not an array in " << jsonPath;
        return kFALSE;
    }
    const Int_t entriesCnt = jsonData.size();
    if (entriesCnt <= 0) {
        LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: empty presigmas array in " << jsonPath;
        return kFALSE;
    }
 
    auto detIdxByID = [&hierarchy](Int_t encID) -> Int_t {
        try {
            return hierarchy.DetModel()->Idx(encID);
        } catch (const std::out_of_range&) {
            return -1;
        }
    };
    auto isPresigmaKey = [](const std::string& k) -> Bool_t {
        return k == "ParamName" || k == "ParamIndex" || k == "DetectorID" || k == "NodeName" || k == "Over"
               || k == "NodeMatch" || k == "Sigma" || k == "Ref" || k == "IDComment";
    };
 
    priors.sigma.assign(hierarchy.AllParams(), 0.0);
    priors.ref.assign(hierarchy.AllParams(), 0.0);
    Int_t errors{0};
 
    for (Int_t i = 0; i < entriesCnt; i++) {
        const auto& jsEntry = jsonData[i];
        if (!jsEntry.is_object()) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " is not an object";
            errors++;
            continue;
        }
        // Reject unknown keys (typo guard), string values are comment-like
        Bool_t badKey{kFALSE};
        for (auto it = jsEntry.begin(); it != jsEntry.end(); ++it) {
            if (it.value().is_string())
                continue;   // comment-like attribute, ignored
            if (!isPresigmaKey(it.key())) {
                LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " has unknown attribute '"
                           << it.key() << "'";
                badKey = kTRUE;
            }
        }
        if (badKey) {
            errors++;
            continue;
        }
 
        const Bool_t hasParamName{jsEntry.contains("ParamName")};
        const Bool_t hasParamIndex{jsEntry.contains("ParamIndex")};
        const Bool_t hasDetectorID{jsEntry.contains("DetectorID")};
        const Bool_t hasNodeName{jsEntry.contains("NodeName")};
        const Bool_t hasOver{jsEntry.contains("Over")};
        const Bool_t hasMatch{jsEntry.contains("NodeMatch")};
 
        if (hasParamName == hasParamIndex) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1
                       << " : exactly one of ParamName / ParamIndex required";
            errors++;
            continue;
        }
        if (hasDetectorID + hasNodeName + hasOver + hasMatch > 1) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1
                       << " : at most one of DetectorID / NodeName / Over / NodeMatch allowed";
            errors++;
            continue;
        }
        if (!jsEntry.contains("Sigma")) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " has no Sigma";
            errors++;
            continue;
        }
        const Double_t sigma{jsEntry["Sigma"].get<Double_t>()};
        if (!(sigma > 0.0)) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " : Sigma must be positive, got "
                       << sigma;
            errors++;
            continue;
        }
        const Double_t ref{jsEntry.contains("Ref") ? jsEntry["Ref"].get<Double_t>() : 0.0};
 
        VecInt_t params;
        if (hasParamName) {
            const std::string what{ParamList(jsEntry["ParamName"], params)};
            if (!what.empty()) {
                LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " : " << what;
                errors++;
                continue;
            }
        } else {
            const Int_t paramIdx{jsEntry["ParamIndex"].get<Int_t>()};
            if (paramIdx < 0 || paramIdx >= BMN_GLOBAL_PARAMS_PD) {
                LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " : ParamIndex " << paramIdx
                           << " out of range [0, " << BMN_GLOBAL_PARAMS_PD << ")";
                errors++;
                continue;
            }
            params.push_back(paramIdx);
        }
 
        NodeList_t nodes;
        std::string what;
        if (hasNodeName) {
            const std::string name{jsEntry["NodeName"].get<std::string>()};
            if (const NodeRef * node{hierarchy.Layout().Find(name)})
                nodes.push_back(node);
            else
                what = "unknown NodeName '" + name + "'";
        } else if (hasDetectorID) {
            const Int_t encID{jsEntry["DetectorID"].get<Int_t>()};
            const Int_t detIdx{detIdxByID(encID)};
            if (detIdx >= 0)
                nodes.push_back(&hierarchy.Layout().ByNode()[detIdx]);
            else
                what = Form("unknown DetectorID %d (%s)", encID, hierarchy.DetModel()->DecodeID(encID).c_str());
        } else if (hasMatch)
            what = MatchComposites(jsEntry["NodeMatch"], hierarchy, nodes);
        else   // no locator at all means every element, as "Over": "Elements" does
            what = LeafScope(hasOver ? jsEntry["Over"] : nlohmann::json("Elements"), hierarchy, nodes);
        if (!what.empty()) {
            LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: presigma #" << i + 1 << " : " << what;
            errors++;
            continue;
        }
        LOG_IF(warning, nodes.empty()) << "  ! presigma #" << i + 1 << " addresses nothing in this model";
 
        for (const Int_t p : params)
            for (const NodeRef* node : nodes) {
                priors.sigma[node->params[p].full] = sigma;
                priors.ref[node->params[p].full] = ref;
            }
    }
 
    if (errors) {
        LOG(error) << " <!> BmnAlignConfig::LoadPresigmas: " << errors << " error(s) detected, aborting";
        return kFALSE;
    }
 
    priors.presigmaCount = 0;
    for (Double_t sigma : priors.sigma)
        if (sigma > 0.0)
            priors.presigmaCount++;
    if (!priors.UsePresigmas()) {
        LOG(warning) << "  ! no parameter got a prior, proceeding without presigmas";
        return kTRUE;
    }
 
    LOG(info) << "    " << priors.presigmaCount << " of " << hierarchy.AllParams() << " parameters got a prior ("
              << entriesCnt << " entries)";
    std::array<Int_t, BMN_GLOBAL_PARAMS_PD> perCorrection{};
    for (const ParamRef& p : hierarchy.Layout().Leaves())
        if (priors.sigma[p.full] > 0.0)
            perCorrection[p.prm]++;
    for (Int_t p = 0; p < BMN_GLOBAL_PARAMS_PD; p++)
        LOG_IF(info, perCorrection[p]) << "      " << BmnParamLayout::ParamName(p) << " : " << perCorrection[p]
                                       << " element(s) constrained";
    for (const NodeRef& node : hierarchy.Layout().ByNode()) {
        if (!node.IsComposite())
            continue;
        TString named;
        for (const ParamRef& p : node.params)
            if (priors.sigma[p.full] > 0.0)
                named += Form("%s=%.3g ", p.prmName, priors.sigma[p.full]);
        LOG(info) << "      composite '" << node.name << "' : " << (named.IsNull() ? "no prior" : named.Data());
    }
    return kTRUE;
}
 
Bool_t BmnAlignConfig::WriteSolution(const char* jsonPath,
                                     const BmnAlignHierarchy& hierarchy,
                                     const BmnAlignResult& result,
                                     Int_t totalTracks,
                                     Bool_t withAlpha)
{
    LOG(state) << " > BmnAlignConfig : saving result to file";
    LOG(info) << "    " << jsonPath << " : file to place result";
 
    std::ofstream ofs(jsonPath);
    if (!ofs.good()) {
        LOG(error) << " <!> BmnAlignConfig::WriteSolution: failed to open output file";
        return kFALSE;
    }
 
    using namespace nlohmann;
    json js;   // all data to be saved
 
    // alignment result
    js["CorrectionValues"] = BmnParamLayout::ParamNames();
    js["DetectorElements"] = hierarchy.DetModel()->Alignables();
    js["CorrectionsPerDetector"] = json::array();
    // Composite nodes of the hierarchy. Not part of the geometry - the element
    // corrections are the *effective* ones and already contain J * Pcomp. What
    // the second block records is how that geometry splits between the levels,
    // which is lost otherwise and cannot be recovered from the elements alone.
    const TVectorD& solutionParams{result.A()};
    if (hierarchy.CompositeCount())
        js["CompositeNodes"] = json::array();
    for (const NodeRef& node : hierarchy.Layout().ByNode()) {
        json element;
        element["values"] = json::array();
        if (node.IsComposite()) {
            element["name"] = node.name;
            element["errors"] = json::array();
            for (const ParamRef& p : node.params) {
                element["values"].push_back(solutionParams(p.full));
                element["errors"].push_back(result.AErr(p.full));
            }
            js["CompositeNodes"].push_back(element);
        } else {
            const Int_t detID{hierarchy.DetModel()->EncodedID(node.leaf)};
            SVectGL values{result.A(node.leaf)};
            if (detID == hierarchy.DetModel()->TargetID())   // see the target note in the header
                values(idZ) -= hierarchy.DetModel()->TargetShift();
            element["detector_id"] = detID;
            element["errors"] = json::array();
            for (const ParamRef& p : node.params) {
                element["values"].push_back(values(p.prm));
                element["errors"].push_back(result.AErr(p.full));
            }
            js["CorrectionsPerDetector"].push_back(element);
        }
    }
 
    // tracks approximation
    if (withAlpha) {
        js["TracksCount"] = totalTracks;
        js["TracksParameters"] = json::array();
        for (const SVectLC& alpha : result.Alpha()) {
            json element = json::array();
            for (Int_t j = 0; j < BMN_LOCAL_PARAMS_PT; j++)
                element.push_back(alpha(j));
            js["TracksParameters"].push_back(element);
        }
    }
 
    // write data
    try {
        ofs << std::setw(4) << js << std::endl;
    } catch (std::exception& e) {
        LOG(error) << " <!> BmnAlignConfig::WriteSolution: failed to write JSON data to file\n" << e.what();
        return kFALSE;
    }
 
    LOG(state) << " < BmnAlignConfig : solution saved";
    LOG(state) << "   ---";
    return kTRUE;
}