Jungfraujoch/image_analysis/geom_refinement/PowderCalibration.h

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only

#pragma once

#include <optional>
#include <string>
#include <vector>

#include "../../common/AzimuthalIntegrationMapping.h"
#include "../../common/DiffractionExperiment.h"
#include "../../common/DiffractionGeometry.h"
#include "../../common/SpotToSave.h"
#include "PowderAutoSeed.h"
#include "RingOptimizer.h"   // RingFitUncertainty

// How the powder rings the detector geometry is fitted to are measured (rugnux --calibration).
enum class CalibrationMethod {
    Rings,   // the run-summed (q x azimuth) azimuthal profile: the ring measured at every azimuth
    Spots    // the pooled per-image spot lists: the ring sampled wherever the spot finder bit
};

struct CalibrationResult {
    DiffractionGeometry geometry;    // the fitted geometry
    size_t ring_points = 0;          // ring measurements the fit used
    // Scatter of those measurements about the fitted rings, and the standard error it implies on the
    // beam centre. Both in pixels - a calibration that has gone wrong (textured ice, one visible ring)
    // says so here, and that is the only warning a user gets.
    double rms_radial_pxl = 0.0;
    double beam_sigma_pxl = 0.0;
    // What the fit itself says about how well each parameter is determined, and how badly the tilt is
    // correlated with the beam centre. rms/beam_sigma above describe the SCATTER of the measurements;
    // this describes the FIT, and the two part company exactly where it matters - a two-ring tilt can
    // leave a small rms while being free to move tens of pixels of beam centre with it.
    RingFitUncertainty uncertainty;
    // The distance the rings themselves asked for before the fit ran, and what the header said. They
    // are reported rather than only used because a large gap between them is the one thing that says
    // the header was wrong - which is usually why the calibration was run at all. Zero when the profile
    // showed too few rings to measure a scale.
    float seed_distance_mm = 0.0f;
    float header_distance_mm = 0.0f;
    // Whether the tilt survived. It is refined by default, but a fit that cannot separate it from the
    // beam centre will still return one - spending tens of pixels of PONI on a tilt the rings do not
    // support - so a tilt smaller than a few times its own uncertainty is declined and the fit repeated
    // with rot1/rot2 pinned. significance is max(|rot| / sigma(rot)) over the two components, 0 where
    // the tilt was never a free parameter.
    bool tilt_refined = false;
    float tilt_significance = 0.0f;
    // How far the beam centre the SPOTS vote for is from the one the rings were fitted to, in pixels,
    // and where the spots put it. An independent cross-check that costs nothing, because the spots have
    // already been found: the circle through them reads nothing from the header, so it holds in exactly
    // the regime the profile does not - a profile binned about a badly wrong centre shows each ring
    // smeared across its sectors, and no seeding recovers that. A large disagreement is therefore not a
    // close call to arbitrate but a statement that this profile could not have been fitted, whatever
    // the residual says. Zero when no spots were available.
    std::optional<DiffractionGeometry> spots_geometry;
    float spots_disagreement_pxl = 0.0f;
};

// How many of its own sigmas the fitted tilt is away from zero - the number the gate above reads.
// Zero when the tilt was held fixed, or when the covariance could not be computed at all (which is
// itself a statement that the tilt is not separable).
float TiltSignificance(const DiffractionGeometry &geom, const RingFitUncertainty &unc);

// Below this the tilt is not a measurement. A tilt the rings genuinely resolve clears it by orders of
// magnitude - measured over a LaB6 distance series, 50, 33, 15 and 8 sigma at 110 to 300 mm - while the
// two-ring fit at 500 mm, whose tilt flips sign and drags the PONI 28 px for no gain in residual,
// reaches 0.1. Any value between 2 and 5 gives the same verdict on all five, so this is a statement
// about which regime a fit is in and not a tuned number.
//
// It is a REJECTION criterion only. Passing it does not certify a tilt: the tilt estimator is
// systematics-limited rather than variance-limited, and a coherent half-pixel error in the ring
// positions fakes a tilt of the size usually seen while leaving sigma small.
constexpr float TILT_MIN_SIGNIFICANCE = 3.0f;

// Both fits take the detector tilt as a free parameter unless refine_tilt is false, which holds
// rot1/rot2 at the value `geom` came in with and fits only the beam centre and the distance. The
// tilt is real and a PONI carries it, but a program that has nowhere to put it - XDS - is better
// given a geometry measured with the tilt pinned than one measured tilted and then flattened.

// Fit the geometry to the rings found in a run-summed azimuthal profile (CalibrationMethod::Rings).
//
// spots, where they are given, add one more beam-centre hypothesis: the centre of the circle through
// them, voted for by every triple (FindCircleCenter). It reads nothing from the header at all - only
// where the spots landed - so it is the one estimate that survives a header the profile itself cannot
// correct from. The rings' own estimate stops working once the beam centre is out by more than about
// half a ring spacing, because past that the azimuthally averaged profile shows each ring as the two
// horns of a sinusoid rather than as a ring; the circle through the spots has no such limit. Measured
// on a 110 mm LaB6 exposure it returns the same geometry from a header 400 px and 8x in distance wrong.
CalibrationResult CalibrateFromProfile(const std::vector<float> &profile,
                                       const AzimuthalIntegrationMapping &mapping,
                                       const DiffractionGeometry &geom,
                                       const std::vector<float> &calibrant_ring_q,
                                       bool refine_tilt = true,
                                       const std::vector<SpotToSave> &spots = {});

// Fit the geometry to a pooled spot list (CalibrationMethod::Spots): the beam centre from scratch off
// the Hough circle centre, then the same ring fit.
CalibrationResult CalibrateFromSpots(const std::vector<SpotToSave> &spots,
                                     const DiffractionGeometry &geom,
                                     const std::vector<float> &calibrant_ring_q,
                                     bool refine_tilt = true);

// Write the geometry as a pyFAI PONI file, the interchange format every azimuthal-integration tool
// reads. Throws if the file cannot be written.
void WritePoniFile(const std::string &path, const DiffractionExperiment &experiment,
                   const DiffractionGeometry &geom);

// The same result as JSON, for whatever reads it next rather than for a person.
//
// A PONI file says where the detector is and nothing else - not how well it is known, not whether the
// tilt in it was measured or merely carried over, and not in the spelling anything in this system
// takes. The object written here has two members. "dataset_settings" holds the geometry under the
// property names broker/jfjoch_api.yaml gives them, and holds nothing else, so it is a valid
// dataset_settings body as it stands and can be POSTed or merged without translating a single field;
// the four the schema requires are always present. "calibration" holds what the run knows about that
// geometry - the residual, the fit's own sigmas, whether the tilt survived its significance test, and
// where the spots independently put the beam - because a calibration that has gone wrong looks exactly
// like one that has not until those are read.
//
// beam_x_pxl/beam_y_pxl is the PONI, as everywhere in this system: on a tilted detector it is NOT where
// the direct beam lands, and that point is reported separately under "calibration". Throws if the file
// cannot be written.
void WriteCalibrationJson(const std::string &path, const DiffractionExperiment &experiment,
                          const CalibrationResult &result,
                          const std::string &calibrant, const std::string &method);