leonarski_fandClaude Opus 5
ef5da29319
A powder calibration is run because the file's geometry is in doubt, so a fit that quietly returns part of that file has answered nothing - and it is indistinguishable from one that worked, down to the residual and the sigmas arranged around it. On one of four LaB6 exposures of one detector the tilt came out at 2.93x its own sigma, a hundredth under the significance gate, so it was declined and pinned - at the master's hardcoded rot1 -0.08, rot2 -0.22 deg. That is eight times the tilt just refused, on no evidence, and worth 10 px of PONI at 190 mm. rugnux printed it to four decimal places, wrote the .poni, and exited 0. Judge the result on provenance instead of on any residual: a geometry is a measurement only if every parameter in it came from this data. Two ways out of the fits do not qualify - a covariance that never conditioned, so the fit cannot say what it determined, and a declined tilt pinned at a non-zero value from the file. A declined tilt over a file stating no tilt still qualifies, because reporting no tilt is then exactly what was measured; so does --no-refine-tilt, because a hold that was asked for is a stated choice and not a silent substitution. No single number separates the four. rms is 2.465 px against 1.44-1.64; the significance of all four lies between 2.93 and 4.47, so the gate is nearly a coin flip at these distances and moving it would only recalibrate on one population; and the failed fit has the TIGHTEST parameter sigmas of the set, because pinning the tilt removes the tilt/centre correlation that inflates a good fit's. The spot cross-check reads 13.5 px against 0.98-2.66, but 10.4 px of that is the pinned tilt moving the PONI - the same defect one step downstream, not independent evidence. On a failure rugnux says so, writes no .poni - a PONI file states where the detector is and has no field in which to say it does not know - writes the JSON with converged false and the reason beside it, and exits non-zero. The re-binning pass now prefers a converged refit over a non-converged one whatever its residual, so a tilt an earlier pass measured is not what a later one gets pinned at. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
572 lines
34 KiB
C++
572 lines
34 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <fstream>
|
|
|
|
#include <nlohmann/json.hpp>
|
|
|
|
#include <spdlog/fmt/fmt.h>
|
|
|
|
#include "PowderCalibration.h"
|
|
#include "../../common/GitInfo.h"
|
|
#include "../../common/JFJochMath.h"
|
|
#include "AssignSpotsToRings.h" // FindCircleCenter
|
|
#include "RingOptimizer.h"
|
|
#include "RingsFromProfile.h"
|
|
|
|
namespace {
|
|
|
|
// How well the ring points sit on the fitted rings, in a unit a user can judge: the radial distance in
|
|
// pixels between where a point is and where the fitted geometry puts its ring. The fit's own residual
|
|
// is in q, so it is divided by the local dq/dr - measured by stepping one pixel outward along the radius
|
|
// rather than assumed, since dq/dr varies with two-theta and with the tilt.
|
|
//
|
|
// The beam centre enters a ring's radius as r(phi) = R + dx cos(phi) + dy sin(phi), so fitting it to n
|
|
// points of scatter s leaves the textbook var = 2 s^2 / n on each of dx and dy. That is the number that
|
|
// separates a beam centre that was measured from one that was merely reported.
|
|
CalibrationResult Summarize(const DiffractionGeometry &fitted,
|
|
const std::vector<RingOptimizerInput> &points,
|
|
const RingFitUncertainty &unc) {
|
|
CalibrationResult result;
|
|
result.geometry = fitted;
|
|
result.uncertainty = unc;
|
|
|
|
const float cx = fitted.GetBeamX_pxl();
|
|
const float cy = fitted.GetBeamY_pxl();
|
|
double sum_sq = 0.0;
|
|
for (const auto &p : points) {
|
|
const float r = std::hypot(p.x - cx, p.y - cy);
|
|
if (!(r > 1.0f))
|
|
continue;
|
|
const float q = fitted.PxlToQ(p.x, p.y);
|
|
const float dq_dr = fitted.PxlToQ(p.x + (p.x - cx) / r, p.y + (p.y - cy) / r) - q;
|
|
if (!(std::abs(dq_dr) > 0.0f))
|
|
continue;
|
|
const double dr = (q - p.q_expected) / dq_dr;
|
|
sum_sq += dr * dr;
|
|
++result.ring_points;
|
|
}
|
|
if (result.ring_points > 0) {
|
|
result.rms_radial_pxl = std::sqrt(sum_sq / static_cast<double>(result.ring_points));
|
|
result.beam_sigma_pxl = result.rms_radial_pxl
|
|
* std::sqrt(2.0 / static_cast<double>(result.ring_points));
|
|
}
|
|
return result;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
// Is the geometry about to be returned a measurement of this data, or does it carry a number out of the
|
|
// input file wearing the appearance of one? A calibration is run precisely because that file's geometry
|
|
// is in doubt, so handing part of it back - to four decimal places, with a residual and a set of sigmas
|
|
// arranged around it - is the one failure a caller cannot see. Nothing downstream can see it either: a
|
|
// .poni states where the detector is and has no field for how that was arrived at.
|
|
//
|
|
// Two ways out of the fits below produce such a geometry. Neither is a threshold; both are statements
|
|
// about where a number came from.
|
|
//
|
|
// The covariance never conditioned. RingOptimizer reports valid = false when the normal matrix is
|
|
// singular at the solution - some direction in parameter space costs the fit nothing at all - and a
|
|
// fit that cannot say what it determined has not determined it.
|
|
//
|
|
// The tilt was declined, and the value pinned in its place is not zero. Declining is the statement
|
|
// that these rings cannot tell a tilt from a shift of the beam centre. Pinning then writes the INPUT
|
|
// FILE's tilt into the answer, which those same rings support no better, and which is not even the
|
|
// tilt that was just refused: measured on a 190 mm powder exposure, the gate rejected a tilt of
|
|
// 0.026 deg at 2.9 sigma and returned the file's 0.22 deg - eight times larger, on no evidence, and
|
|
// worth 10 px of PONI at that distance. Where the file's tilt is zero the two agree and the output is
|
|
// honest: the fit found no tilt and reports none, which is what a long-distance fit should do. Where
|
|
// it is not zero, the answer carries an angle from the header at whatever size the header stated.
|
|
//
|
|
// --no-refine-tilt is not this. There the user asked for the file's tilt to be held, and a held tilt
|
|
// that was asked for is a stated choice rather than a silent substitution.
|
|
void JudgeCalibration(CalibrationResult &result, const DiffractionGeometry &header,
|
|
bool refine_tilt) {
|
|
if (!result.uncertainty.valid) {
|
|
result.converged = false;
|
|
result.reason = "the fit is degenerate at its solution - its covariance does not condition, so "
|
|
"it cannot say what it determined";
|
|
return;
|
|
}
|
|
if (refine_tilt && !result.tilt_refined
|
|
&& (header.GetPoniRot1_rad() != 0.0f || header.GetPoniRot2_rad() != 0.0f)) {
|
|
constexpr double RAD_TO_DEG = 180.0 / PI;
|
|
const double poni_pxl = std::hypot(header.GetPoniRot1_rad(), header.GetPoniRot2_rad())
|
|
* result.geometry.GetDetectorDistance_mm() / result.geometry.GetPixelSize_mm();
|
|
result.converged = false;
|
|
result.reason = fmt::format(
|
|
"the tilt was declined at {:.1f}x its own sigma and pinned at the input file's "
|
|
"rot1 {:+.4f} deg, rot2 {:+.4f} deg - a tilt these rings measured no better than the "
|
|
"one they refused, and worth {:.1f} px of PONI at this distance",
|
|
result.tilt_significance,
|
|
header.GetPoniRot1_rad() * RAD_TO_DEG, header.GetPoniRot2_rad() * RAD_TO_DEG,
|
|
poni_pxl);
|
|
}
|
|
}
|
|
|
|
float TiltSignificance(const DiffractionGeometry &geom, const RingFitUncertainty &unc) {
|
|
if (!unc.valid)
|
|
return 0.0f;
|
|
float significance = 0.0f;
|
|
if (unc.sigma_rot1_rad > 0.0)
|
|
significance = std::max(significance,
|
|
std::abs(geom.GetPoniRot1_rad()) / static_cast<float>(unc.sigma_rot1_rad));
|
|
if (unc.sigma_rot2_rad > 0.0)
|
|
significance = std::max(significance,
|
|
std::abs(geom.GetPoniRot2_rad()) / static_cast<float>(unc.sigma_rot2_rad));
|
|
return significance;
|
|
}
|
|
|
|
CalibrationResult CalibrateFromProfile(const std::vector<float> &profile,
|
|
const AzimuthalIntegrationMapping &mapping,
|
|
const DiffractionGeometry &geom,
|
|
const std::vector<float> &calibrant_ring_q,
|
|
bool refine_tilt,
|
|
const std::vector<SpotToSave> &spots) {
|
|
// Where to start. The ring search below is local - each ring is looked for inside a window a few
|
|
// pixels of radius wide - so a header distance more than a percent or two out puts every ring
|
|
// outside its own window, and what the fit then converges on is noise. Ask the rings what the
|
|
// distance is rather than believing the header (see PowderAutoSeed.h), and because a powder pattern
|
|
// has genuine distance aliases, ask for several answers and fit them all.
|
|
const auto observed = RingRadiiFromProfile(profile, mapping, geom);
|
|
const auto [radius_min, radius_max] = ProfileRadiusRange_pxl(mapping, geom);
|
|
const auto candidates = CandidateDistancesFromPowderRings(observed, calibrant_ring_q, geom,
|
|
radius_min, radius_max);
|
|
|
|
// ...and ask them where the beam is, for the same reason. The extraction looks for each ring within
|
|
// a window a few pixels of radius wide, so a header centre more than about ten pixels out puts the
|
|
// ring outside that window over much of the turn - and the fit then reads a cos(phi) signal off
|
|
// whatever sectors are left, which is how a 20 px error used to end 31 px wrong. The rings answer
|
|
// this without a calibrant and without a distance: a powder ring is a conic centred on the beam, so
|
|
// a wrong centre makes EVERY ring's radius oscillate once per turn by the same amount.
|
|
// Three places a beam centre can come from, and they fail in different regimes - which is the whole
|
|
// reason to carry all of them. The header is right on a well-configured instrument and is what a
|
|
// calibration is run to check. The rings' own wobble is exact while the error stays under about half
|
|
// a ring spacing, and stops meaning anything beyond that, where each ring reaches the azimuthally
|
|
// averaged profile as two horns rather than one peak. The circle through the spots reads nothing
|
|
// from the header whatsoever, so it holds where both of the others have given up - measured, it
|
|
// returns the same answer from a header 400 px out in the centre and eight times out in distance.
|
|
// None of them is trusted; each is fitted and the residual chooses.
|
|
std::vector<DiffractionGeometry> centres = {geom};
|
|
const auto add_centre = [&](float x, float y) {
|
|
// Within a pixel of one already on the list it is not another hypothesis, it is the same one,
|
|
// and fitting it again would only double the work for two answers nothing can tell apart.
|
|
for (const auto &c : centres)
|
|
if (std::hypot(x - c.GetBeamX_pxl(), y - c.GetBeamY_pxl()) < 1.0f)
|
|
return;
|
|
DiffractionGeometry candidate = geom;
|
|
centres.push_back(candidate.BeamX_pxl(x).BeamY_pxl(y));
|
|
};
|
|
|
|
if (const auto offset = BeamCentreOffsetFromProfile(profile, mapping, geom, observed))
|
|
add_centre(geom.GetBeamX_pxl() + offset->first, geom.GetBeamY_pxl() + offset->second);
|
|
|
|
// What the spots make of it, as a COMPLETE geometry rather than only a centre. Taking the circle
|
|
// centre alone is not enough: the distance candidates are read off the azimuthally averaged profile,
|
|
// which is averaged about the header's centre, and once that is badly wrong the profile no longer
|
|
// shows rings at all - a ring smeared over hundreds of pixels cannot be un-smeared by reading its
|
|
// bins differently. The spots have no such problem, because a spot's position is a fact about the
|
|
// image; GuessGeometry votes for the circle centre, clusters the radii into rings and takes the
|
|
// distance from the innermost one, which is exactly the path that makes --calibration spots immune
|
|
// to the header. Let it produce one whole starting geometry of its own.
|
|
std::optional<DiffractionGeometry> from_spots;
|
|
if (spots.size() >= 3) {
|
|
try {
|
|
DiffractionGeometry guess = geom;
|
|
GuessGeometry(guess, spots, calibrant_ring_q, refine_tilt);
|
|
from_spots = guess;
|
|
} catch (const std::exception &) {
|
|
// No circle, or no ring clusters - the spots simply have nothing to say here. The profile's
|
|
// own hypotheses stand on their own, so there is nothing to report and nothing to stop.
|
|
}
|
|
}
|
|
|
|
// One starting distance, fitted to convergence. A seed only has to land in the fit's basin, not on
|
|
// the answer: it is measured from blended peaks in the azimuthally-averaged profile and is good to
|
|
// about a per cent, which is close enough to converge from but far enough to sit every search window
|
|
// a few pixels off its ring - and an off-centre window takes its background off the ring's own
|
|
// flank, which costs both points and residual. So re-extract at the geometry the fit converged to
|
|
// and fit again. Nothing is re-read from disk, so the whole loop is free.
|
|
struct Attempt {
|
|
DiffractionGeometry geometry;
|
|
std::vector<RingOptimizerInput> points;
|
|
RingFitUncertainty uncertainty;
|
|
double rms_radial_pxl = 0.0;
|
|
};
|
|
const auto fit_from = [&](const DiffractionGeometry &start, bool seeded,
|
|
bool tilt) -> std::vector<Attempt> {
|
|
DiffractionGeometry current = start;
|
|
Attempt attempt;
|
|
constexpr int MAX_PASSES = 3;
|
|
for (int pass = 0; pass < MAX_PASSES; ++pass) {
|
|
// The profile's own geometry is the search on the first pass of an unseeded attempt, and
|
|
// the geometry converged to on every pass after - which is where the rings have been
|
|
// measured to be, rather than where the header guessed.
|
|
const DiffractionGeometry *search = (pass == 0 && !seeded) ? nullptr : ¤t;
|
|
auto pass_points = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q,
|
|
0.06f, 3.0f, search);
|
|
if (pass_points.empty())
|
|
break;
|
|
|
|
RingFitUncertainty pass_unc;
|
|
const auto pass_fitted = RingOptimizer(current, tilt).Run(pass_points, &pass_unc);
|
|
const float moved = std::abs(pass_fitted.GetDetectorDistance_mm()
|
|
- current.GetDetectorDistance_mm());
|
|
attempt.points = std::move(pass_points);
|
|
attempt.uncertainty = pass_unc;
|
|
attempt.geometry = pass_fitted;
|
|
current = pass_fitted;
|
|
// A tenth of a micron of distance moves the outermost ring by far less than a thousandth of
|
|
// a pixel, so there is nothing left for another pass to find.
|
|
if (moved < 1e-4f)
|
|
break;
|
|
}
|
|
if (attempt.points.empty())
|
|
return {};
|
|
|
|
// Two ways to take the final measurement, and they genuinely disagree about which is better.
|
|
//
|
|
// Following the rings sector by sector is what lets a badly placed beam centre be recovered at
|
|
// all - a centre wrong by (dx, dy) puts a ring at a different q in every sector, and one window
|
|
// centred on one q finds it only where that oscillation happens to be small. But it costs
|
|
// precision once they have been found, for a reason worth stating: a window that moves with phi
|
|
// makes every systematic of the peak finder - where the background line is taken, how the
|
|
// centroid sits in the window - vary with phi too, and phi is exactly the axis the beam centre
|
|
// is read off. Measured, rms 0.415 -> 0.525 px on a good 110 mm fit, and 0.831 when the window
|
|
// followed the fitted tilt as well.
|
|
//
|
|
// The alternative is a window that is the same in every sector, which the binned geometry with
|
|
// only its DISTANCE replaced gives by construction - a distance error moves every sector's
|
|
// window equally, where a centre or a tilt error does not. That is more precise where it works
|
|
// and finds nothing where the centre is far out. So do both and let the same rule that ranks
|
|
// everything else here decide.
|
|
std::vector<Attempt> out;
|
|
out.push_back(attempt);
|
|
|
|
DiffractionGeometry measure_geom = geom;
|
|
measure_geom.DetectorDistance_mm(current.GetDetectorDistance_mm());
|
|
auto measured = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q,
|
|
0.06f, 3.0f, &measure_geom);
|
|
if (!measured.empty()) {
|
|
Attempt fixed_window;
|
|
RingFitUncertainty measured_unc;
|
|
fixed_window.geometry = RingOptimizer(current, tilt).Run(measured, &measured_unc);
|
|
fixed_window.points = std::move(measured);
|
|
fixed_window.uncertainty = measured_unc;
|
|
out.push_back(std::move(fixed_window));
|
|
}
|
|
|
|
for (auto &a : out)
|
|
a.rms_radial_pxl = Summarize(a.geometry, a.points, a.uncertainty).rms_radial_pxl;
|
|
return out;
|
|
};
|
|
|
|
// Every combination of what the rings said and what the header said, fitted, with the residual left
|
|
// to choose. The header is a hypothesis like any other here, neither trusted nor discarded - and so
|
|
// is the seeded beam centre, which is NOT simply better than the header's.
|
|
//
|
|
// The centre seed reads the once-per-turn wobble of the ring radii, and a tilt puts a term of that
|
|
// same shape there too - one that grows as the ring's radius squared, where a centre error does not.
|
|
// Pooling the rings into one offset therefore absorbs part of the tilt into the centre, which is
|
|
// worth several pixels on a genuinely tilted detector and made a good 110 mm fit worse when it was
|
|
// simply believed. What it buys is capture range, and only where the header centre is far enough out
|
|
// that the extraction would otherwise find the rings over a fraction of the turn. Offering it as an
|
|
// alternative start costs one more fit each and needs no rule about when it applies.
|
|
struct Provenance { float seed_mm; double match; };
|
|
struct Start { DiffractionGeometry geometry; bool tracked; Provenance provenance; };
|
|
|
|
std::vector<Start> starts;
|
|
// The spots' geometry is a start in its own right, not one to cross with the profile's distances -
|
|
// its centre and its distance are measured together and belong together.
|
|
if (from_spots)
|
|
starts.push_back({*from_spots, true, {from_spots->GetDetectorDistance_mm(), 0.0}});
|
|
for (size_t centre = 0; centre < centres.size(); ++centre) {
|
|
for (size_t i = 0; i <= candidates.size(); ++i) {
|
|
const bool seeded_distance = i < candidates.size();
|
|
DiffractionGeometry start = centres[centre];
|
|
if (seeded_distance)
|
|
start.DetectorDistance_mm(candidates[i].distance_mm);
|
|
starts.push_back({start, seeded_distance || centre > 0,
|
|
{seeded_distance ? candidates[i].distance_mm : 0.0f,
|
|
seeded_distance ? candidates[i].score : 0.0}});
|
|
}
|
|
}
|
|
|
|
std::vector<std::pair<Attempt, Provenance>> attempts;
|
|
for (const auto &start : starts)
|
|
for (auto &attempt : fit_from(start.geometry, start.tracked, refine_tilt))
|
|
attempts.emplace_back(std::move(attempt), start.provenance);
|
|
|
|
// Residual alone cannot rank these: a starting distance so wrong that only one ring point survives
|
|
// leaves a residual of exactly zero, and would win every time. How many ring measurements an
|
|
// attempt explains is evidence in its own right, and the first thing to compare - an attempt that
|
|
// accounts for half as much of the pattern is not in the running whatever it does with what is
|
|
// left. Among those that explain a comparable amount, the residual decides, and decides clearly,
|
|
// because only the true distance makes every ring fit at once - on the aliased LaB6 case the two
|
|
// candidates differ by 0.4 px against 5.2 px, which is not a close call.
|
|
size_t most_points = 0;
|
|
for (const auto &[attempt, provenance] : attempts)
|
|
most_points = std::max(most_points, attempt.points.size());
|
|
|
|
std::optional<Attempt> best;
|
|
Provenance best_provenance{};
|
|
for (auto &[attempt, provenance] : attempts) {
|
|
if (attempt.points.size() * 2 < most_points)
|
|
continue;
|
|
if (!best || attempt.rms_radial_pxl < best->rms_radial_pxl) {
|
|
best = std::move(attempt);
|
|
best_provenance = provenance;
|
|
}
|
|
}
|
|
|
|
if (!best)
|
|
throw JFJochException(JFJochExceptionCategory::CalibrationError,
|
|
"No powder ring found in the summed azimuthal profile");
|
|
|
|
// Did the tilt earn its place? A single ring cannot separate a tilt from a beam-centre shift at all
|
|
// - both move a ring's radius as cos(phi), and only how that amplitude grows with the ring's radius
|
|
// tells them apart - and two barely-sampled rings cannot either. The fit still returns a tilt in
|
|
// that case, because nothing stopped it, and pays for it with the beam centre: at 500 mm on the LaB6
|
|
// series it swings to the opposite sign of every shorter distance, drags the PONI 28 px, and buys a
|
|
// residual of 0.960 px against 0.962 px pinned. The covariance says so plainly - 0.1 sigma, and a
|
|
// PONI quoted to +-180 px - so ask it, and where the answer is no, fit again with the tilt held.
|
|
// The rule is on the geometry that gets REPORTED, not on what the last fit happened to do: a tilt
|
|
// may only survive if the final fit had it free AND it cleared the test. Both halves are needed.
|
|
// RingOptimizer pins the tilt by itself when every point it is given lies on one ring, and on a
|
|
// barely-sampled pattern a later pass can land in exactly that state - which used to freeze the tilt
|
|
// at whatever the FIRST pass produced and hand it back with sigma zero, i.e. an unmeasured tilt
|
|
// wearing the appearance of a fixed one. Refitting from the header's tilt is what makes the
|
|
// reported geometry honest in both cases.
|
|
const bool tilt_was_free = best->uncertainty.valid
|
|
&& (best->uncertainty.sigma_rot1_rad > 0.0 || best->uncertainty.sigma_rot2_rad > 0.0);
|
|
const float significance = tilt_was_free ? TiltSignificance(best->geometry, best->uncertainty) : 0.0f;
|
|
bool tilt_refined = refine_tilt && tilt_was_free && significance >= TILT_MIN_SIGNIFICANCE;
|
|
// Not `refine_tilt && !tilt_refined`: with --no-refine-tilt the tilt is never free, so
|
|
// tilt_refined is false and this branch was skipped entirely - leaving the rotations at the ZERO
|
|
// that GuessInitialGeometry sets, while the CLI reported them as held at the header value. A
|
|
// declined tilt and a tilt nobody asked to refine both need the header put back and the rest
|
|
// refitted around it.
|
|
if (!tilt_refined) {
|
|
// From the header's tilt, not from the one being declined. fit_from now takes a whole geometry,
|
|
// so without this the "pinned" refit would pin rot1/rot2 at exactly the unvalidated values the
|
|
// gate just rejected - which is the same fault the gate exists to catch, reintroduced one level
|
|
// up.
|
|
DiffractionGeometry pinned_start = best->geometry;
|
|
pinned_start.PoniRot1_rad(geom.GetPoniRot1_rad()).PoniRot2_rad(geom.GetPoniRot2_rad());
|
|
std::optional<Attempt> pinned;
|
|
for (auto &a : fit_from(pinned_start, true, false))
|
|
if (!pinned || a.points.size() > pinned->points.size()
|
|
|| (a.points.size() == pinned->points.size()
|
|
&& a.rms_radial_pxl < pinned->rms_radial_pxl))
|
|
pinned = std::move(a);
|
|
if (pinned)
|
|
best = std::move(*pinned);
|
|
}
|
|
|
|
auto result = Summarize(best->geometry, best->points, best->uncertainty);
|
|
if (from_spots) {
|
|
result.spots_geometry = from_spots;
|
|
result.spots_disagreement_pxl =
|
|
std::hypot(from_spots->GetBeamX_pxl() - best->geometry.GetBeamX_pxl(),
|
|
from_spots->GetBeamY_pxl() - best->geometry.GetBeamY_pxl());
|
|
}
|
|
result.tilt_refined = tilt_refined;
|
|
result.tilt_significance = significance;
|
|
result.seed_distance_mm = best_provenance.seed_mm;
|
|
result.header_distance_mm = geom.GetDetectorDistance_mm();
|
|
JudgeCalibration(result, geom, refine_tilt);
|
|
return result;
|
|
}
|
|
|
|
CalibrationResult CalibrateFromSpots(const std::vector<SpotToSave> &spots,
|
|
const DiffractionGeometry &geom,
|
|
const std::vector<float> &calibrant_ring_q,
|
|
bool refine_tilt) {
|
|
DiffractionGeometry fitted = geom;
|
|
// From scratch (Hough circle centre + ring clustering), then refined: the guess pins the centre to a
|
|
// whole pixel and only sees the spots its clustering kept, so the refine re-matches every spot at
|
|
// that geometry.
|
|
RingFitUncertainty unc;
|
|
GuessGeometry(fitted, spots, calibrant_ring_q, refine_tilt);
|
|
OptimizeGeometry(fitted, spots, calibrant_ring_q, refine_tilt, &unc);
|
|
|
|
// The same question the profile path asks, for the same reason and with the same answer if no: a
|
|
// tilt may only be reported if this fit had it free and it stands clear of its own uncertainty.
|
|
// Pinning it means putting rot1/rot2 back where the geometry came in and refitting the beam centre
|
|
// and distance around that, not simply deleting an angle from the answer.
|
|
const bool tilt_was_free = unc.valid
|
|
&& (unc.sigma_rot1_rad > 0.0 || unc.sigma_rot2_rad > 0.0);
|
|
const float significance = tilt_was_free ? TiltSignificance(fitted, unc) : 0.0f;
|
|
const bool tilt_refined = refine_tilt && tilt_was_free && significance >= TILT_MIN_SIGNIFICANCE;
|
|
// Not `refine_tilt && !tilt_refined`: with --no-refine-tilt the tilt is never free, so
|
|
// tilt_refined is false and this branch was skipped entirely - leaving the rotations at the ZERO
|
|
// that GuessInitialGeometry sets, while the CLI reported them as held at the header value. A
|
|
// declined tilt and a tilt nobody asked to refine both need the header put back and the rest
|
|
// refitted around it.
|
|
if (!tilt_refined) {
|
|
fitted.PoniRot1_rad(geom.GetPoniRot1_rad()).PoniRot2_rad(geom.GetPoniRot2_rad());
|
|
OptimizeGeometry(fitted, spots, calibrant_ring_q, false, &unc);
|
|
}
|
|
|
|
auto result = Summarize(fitted, AssignSpotsToRings(fitted, spots, calibrant_ring_q), unc);
|
|
result.tilt_refined = tilt_refined;
|
|
result.tilt_significance = significance;
|
|
result.header_distance_mm = geom.GetDetectorDistance_mm();
|
|
JudgeCalibration(result, geom, refine_tilt);
|
|
return result;
|
|
}
|
|
|
|
void WritePoniFile(const std::string &path, const DiffractionExperiment &experiment,
|
|
const DiffractionGeometry &geom) {
|
|
std::ofstream f(path);
|
|
if (!f)
|
|
throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot write " + path);
|
|
|
|
const double pixel_m = geom.GetPixelSize_mm() * 1e-3;
|
|
// pyFAI's axis convention is the trap: Poni1 (and pixel1) is the SLOW axis - rows, our y - and
|
|
// Poni2 the FAST axis - columns, our x - both in metres from the detector origin. A transposed PONI
|
|
// file is silently wrong, so the mapping is spelled out here rather than left to the reader.
|
|
//
|
|
// DiffractionGeometry's beam_x/beam_y IS the PONI: LabCoord rotates the vector measured FROM that
|
|
// pixel, i.e. it is the point of normal incidence, so it maps straight across with no correction.
|
|
// GetDirectBeam_pxl() is a different quantity - where the direct beam lands - and parts from the
|
|
// PONI as soon as rot1/rot2 are non-zero.
|
|
//
|
|
// The half pixel is the origin convention (see docs/DETECTOR_GEOMETRY.md): our coordinates are
|
|
// pixel-centred, so beam_x = 948 means the CENTRE of pixel 948, while pyFAI measures from the edge
|
|
// of the sensor and puts the centre of pixel i at (i + 0.5) * pixel size. Without it the pattern
|
|
// pyFAI integrates sits half a pixel off ours.
|
|
const double half_pixel_m = 0.5 * pixel_m;
|
|
f << fmt::format("# Calibration done by Jungfraujoch rugnux {}\n", jfjoch_version());
|
|
// poni_version 2.1 is what pyFAI introduced "orientation" with (pyFAI 2024.01).
|
|
f << "poni_version: 2.1\n";
|
|
f << "Detector: Detector\n";
|
|
// orientation 2 is pyFAI's "origin at the top left of the image when looking FROM the sample",
|
|
// which is the MX convention Jungfraujoch assembles to. Without it pyFAI assumes its own default,
|
|
// orientation 3 (bottom left), and quietly believes increasing row means physically upwards. The
|
|
// radial integration is identical either way - a mirror preserves 2theta - but the azimuth comes
|
|
// out with the opposite sense, which matters for anything that uses chi (cake or sector
|
|
// integration, texture).
|
|
f << fmt::format("Detector_config: {{\"pixel1\": {:g}, \"pixel2\": {:g}, \"max_shape\": [{}, {}], "
|
|
"\"orientation\": 2}}\n",
|
|
pixel_m, pixel_m, experiment.GetYPixelsNumConv(), experiment.GetXPixelsNumConv());
|
|
f << fmt::format("Distance: {:.9g}\n", geom.GetDetectorDistance_mm() * 1e-3);
|
|
// Poni1 is measured from pyFAI's own origin, so declaring orientation 2 re-anchors it to the top
|
|
// edge: the same physical point is now (height - 1 - beam_y) rows down from there.
|
|
f << fmt::format("Poni1: {:.9g}\n",
|
|
(experiment.GetYPixelsNumConv() - 1 - geom.GetBeamY_pxl()) * pixel_m + half_pixel_m);
|
|
f << fmt::format("Poni2: {:.9g}\n", geom.GetBeamX_pxl() * pixel_m + half_pixel_m);
|
|
// With orientation declared, rot2 and rot3 change sign and rot1 does not, and Rot3 carries a
|
|
// further half turn:
|
|
// (Rot1, Rot2, Rot3) = (+rot1, +rot2, -rot3 + pi)
|
|
// A row flip is an improper transformation, so it reverses the sense of rotations about x and
|
|
// about the beam while leaving the one about the vertical alone. The half turn is the azimuthal
|
|
// reference: pyFAI's in-plane axes are the negatives of ours, so without it every chi comes out
|
|
// 180 degrees away. It is a rotation about the beam, so it leaves 2theta untouched - which is
|
|
// why radial integration was right all along and only the azimuth was wrong.
|
|
// Pinned empirically against pyFAI 2026.5.0 on a tilted detector (4/-6.5/13 deg, off-centre
|
|
// beam), against the lab positions of the NXmx chain: 2theta to 3.6e-15 deg and chi to 2.8e-14
|
|
// deg over the whole detector. The half turn is needed for the orientation-3 form written before
|
|
// this too, so it is not an artefact of declaring the orientation.
|
|
f << fmt::format("Rot1: {:.9g}\n", geom.GetPoniRot1_rad());
|
|
// negate() rather than a bare minus so an unrefined angle prints as 0 and not -0.
|
|
const auto negate = [](float v) { return v == 0.0f ? 0.0f : -v; };
|
|
f << fmt::format("Rot2: {:.9g}\n", geom.GetPoniRot2_rad());
|
|
f << fmt::format("Rot3: {:.9g}\n", negate(geom.GetPoniRot3_rad()) + PI);
|
|
f << fmt::format("Wavelength: {:.9g}\n", geom.GetWavelength_A() * 1e-10);
|
|
f.flush();
|
|
if (!f)
|
|
throw JFJochException(JFJochExceptionCategory::FileWriteError, "Error writing " + path);
|
|
}
|
|
|
|
namespace {
|
|
// Rounded, because a float promoted to double serialises as 765.90197753906250 and the extra figures
|
|
// are noise from the conversion rather than anything the fit measured.
|
|
double Rounded(double value, int decimals) {
|
|
const double scale = std::pow(10.0, decimals);
|
|
return std::round(value * scale) / scale;
|
|
}
|
|
} // namespace
|
|
|
|
void WriteCalibrationJson(const std::string &path, const DiffractionExperiment &experiment,
|
|
const CalibrationResult &result,
|
|
const std::string &calibrant, const std::string &method) {
|
|
const auto &g = result.geometry;
|
|
|
|
// The geometry, spelled as dataset_settings spells it and holding nothing that is not one of its
|
|
// properties. The four the schema requires are unconditional.
|
|
nlohmann::json settings;
|
|
settings["beam_x_pxl"] = Rounded(g.GetBeamX_pxl(), 3);
|
|
settings["beam_y_pxl"] = Rounded(g.GetBeamY_pxl(), 3);
|
|
settings["detector_distance_mm"] = Rounded(g.GetDetectorDistance_mm(), 4);
|
|
settings["incident_energy_keV"] = Rounded(experiment.GetDatasetSettings().GetPhotonEnergy_keV(), 6);
|
|
// The rotations ride along whenever any of them is non-zero. Omitting them would not leave the
|
|
// geometry unstated - it would state a FLAT detector, since the API's own default is 0.0 - so they
|
|
// go in together or not at all, rot3 included even though nothing here refines it.
|
|
if (g.GetPoniRot1_rad() != 0.0f || g.GetPoniRot2_rad() != 0.0f || g.GetPoniRot3_rad() != 0.0f) {
|
|
settings["poni_rot1_rad"] = Rounded(g.GetPoniRot1_rad(), 9);
|
|
settings["poni_rot2_rad"] = Rounded(g.GetPoniRot2_rad(), 9);
|
|
settings["poni_rot3_rad"] = Rounded(g.GetPoniRot3_rad(), 9);
|
|
}
|
|
|
|
nlohmann::json calibration;
|
|
// The verdict first, because everything under it is only worth reading once it is known which of the
|
|
// two this file is. A caller that reads nothing else must still not mistake a non-fit for a fit.
|
|
calibration["converged"] = result.converged;
|
|
if (!result.converged)
|
|
calibration["not_converged_reason"] = result.reason;
|
|
calibration["calibrant"] = calibrant;
|
|
calibration["method"] = method;
|
|
calibration["ring_points"] = result.ring_points;
|
|
calibration["rms_radial_pxl"] = Rounded(result.rms_radial_pxl, 4);
|
|
calibration["beam_sigma_pxl"] = Rounded(result.beam_sigma_pxl, 4);
|
|
// Where the beam actually lands, which is not the PONI above once the detector is tilted. Both are
|
|
// written because programs disagree about which one they mean by "beam centre".
|
|
const auto [direct_x, direct_y] = g.GetDirectBeam_pxl();
|
|
calibration["direct_beam_x_pxl"] = Rounded(direct_x, 3);
|
|
calibration["direct_beam_y_pxl"] = Rounded(direct_y, 3);
|
|
calibration["header_distance_mm"] = Rounded(result.header_distance_mm, 4);
|
|
calibration["tilt_refined"] = result.tilt_refined;
|
|
calibration["tilt_significance"] = Rounded(result.tilt_significance, 2);
|
|
if (result.seed_distance_mm > 0.0f)
|
|
calibration["ring_seed_distance_mm"] = Rounded(result.seed_distance_mm, 4);
|
|
|
|
if (result.uncertainty.valid) {
|
|
nlohmann::json sigma;
|
|
sigma["beam_x_pxl"] = Rounded(result.uncertainty.sigma_beam_x_pxl, 4);
|
|
sigma["beam_y_pxl"] = Rounded(result.uncertainty.sigma_beam_y_pxl, 4);
|
|
sigma["detector_distance_mm"] = Rounded(result.uncertainty.sigma_distance_mm, 5);
|
|
if (result.uncertainty.sigma_rot1_rad > 0.0 || result.uncertainty.sigma_rot2_rad > 0.0) {
|
|
sigma["poni_rot1_rad"] = Rounded(result.uncertainty.sigma_rot1_rad, 9);
|
|
sigma["poni_rot2_rad"] = Rounded(result.uncertainty.sigma_rot2_rad, 9);
|
|
sigma["correlation_beam_x_rot1"] = Rounded(result.uncertainty.corr_beam_x_rot1, 4);
|
|
sigma["correlation_beam_y_rot2"] = Rounded(result.uncertainty.corr_beam_y_rot2, 4);
|
|
}
|
|
calibration["fit_sigma"] = sigma;
|
|
}
|
|
|
|
if (result.spots_geometry) {
|
|
nlohmann::json cross;
|
|
cross["beam_x_pxl"] = Rounded(result.spots_geometry->GetBeamX_pxl(), 3);
|
|
cross["beam_y_pxl"] = Rounded(result.spots_geometry->GetBeamY_pxl(), 3);
|
|
cross["detector_distance_mm"] = Rounded(result.spots_geometry->GetDetectorDistance_mm(), 4);
|
|
cross["disagreement_pxl"] = Rounded(result.spots_disagreement_pxl, 3);
|
|
calibration["spot_cross_check"] = cross;
|
|
}
|
|
|
|
nlohmann::json out;
|
|
out["dataset_settings"] = settings;
|
|
out["calibration"] = calibration;
|
|
out["jfjoch_version"] = jfjoch_version();
|
|
|
|
std::ofstream f(path);
|
|
if (!f)
|
|
throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot write " + path);
|
|
f << out.dump(4) << "\n";
|
|
f.flush();
|
|
if (!f)
|
|
throw JFJochException(JFJochExceptionCategory::FileWriteError, "Error writing " + path);
|
|
}
|