leonarski_fandClaude Opus 5
abfdb89b79
The centre in the file is often a placeholder, and nothing measures it until post-refinement has already indexed the sweep - by which time a wrong centre has chosen the lattice. Two exact facts about a rotation sweep give it from spot positions alone, with no cell, no orientation matrix and nothing indexed. Rotating 180 degrees about the spindle and taking -h negates a reflection's component along the spindle and leaves the rest, so with the spindle perpendicular to the beam the Laue condition is preserved and the spots recorded half a turn apart are mirror images along the spindle. Those are Friedel mates, not the same reflection. The same reflection appears twice for a different reason: it meets the Ewald sphere on two crossings, generally not half a turn apart, differing only in the sign of the component perpendicular to both the spindle and the beam. The first observable gives the coordinate along the spindle, the second the coordinate across it. Each candidate pairing votes and the true value accumulates while wrong pairings scatter. Both observables need guarding, because a vote is a comb and the tallest tooth is not always the right one. Along the spindle a false pairing cannot fake the equality of Friedel amplitudes. Across it, the two crossings of one reflection are separated by a sweep angle its own position fixes, which no accidental pair reproduces. The mirror is exact in the laboratory frame, so it is only as good as the rotation axis. Every file here states an ideal axis and none of them has one; a skew about the beam spreads the vote instead of shifting it, and past a milliradian it moves an otherwise correct answer by pixels while every internal statistic still looks healthy. It is therefore fitted, not assumed. A tilt of the axis towards the beam is measured and reported but not applied, being confounded with the detector rotation until that is fitted too. Nothing inside the fit can see a wrong tooth - when the vote flips, every frame pair flips with it - so the answer is checked from outside, by asking whether it depends on where the search began. That, and a floor on the angular span the pairs cover, are what refuse the cases this cannot measure: a sweep barely past half a turn is the dangerous one, not the short one, because at exactly half a turn there is nothing to fit and just past it there is almost nothing. Where the sweep is too short for any of this the radial background profile gives a coarser centre from a handful of images, and where neither can measure it the file's value is kept. The beam-stop projection now takes its own frames rather than sharing the sample, so turning this on cannot change the mask; and both samples keep away from the ends of the sweep, where shutter synchronisation spoils an image. Reading twice as many frames as before costs a few seconds once, and is what makes the answer independent of which frames were drawn. Off by default. Over the 38-crystal rotation battery it serves every dataset, agrees with XDS's refined direct beam to 0.116 px in the median against 0.135 for the value in the file, and changes no space group. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
38 lines
1.9 KiB
C++
38 lines
1.9 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include <optional>
|
|
#include <vector>
|
|
|
|
#include "../../common/DiffractionExperiment.h"
|
|
#include "../../common/PixelMask.h"
|
|
|
|
struct BeamCenterEstimate {
|
|
float beam_x_pxl = 0.0f;
|
|
float beam_y_pxl = 0.0f;
|
|
float sigma_pxl = 0.0f; // 1 sigma on the fitted shift, the larger of the two axes
|
|
};
|
|
|
|
// Beam centre from the isotropy of the scattered background, before anything is indexed.
|
|
//
|
|
// The solvent and air scatter is isotropic in 2-theta about the beam, so a centre that is off
|
|
// shifts each azimuthal sector's radial profile by a different amount. Sector k's profile is
|
|
// m_k * g(2theta + d_k), with the shift d_k = Jx_k*dx + Jy_k*dy and m_k an amplitude that
|
|
// absorbs anything multiplicative and azimuthal - a holder arm, a cryostream shadow, a
|
|
// flat-field gradient. Fitting the amplitude alongside the shift is what makes this usable:
|
|
// a 50% shadow over one sextant otherwise reads as several tens of pixels of centre error.
|
|
//
|
|
// The leverage comes from the CURVATURE of the radial profile - the water ring - because for a
|
|
// pure exponential decay g' is proportional to g and shift and amplitude are indistinguishable.
|
|
// The sigma measures that leverage, so it grows as the curvature weakens, and the caller's gate on
|
|
// it is what keeps an ill-determined centre out. It is a precision and not an accuracy: on a
|
|
// background with NO curvature at all there is nothing to separate the two parameters, the fit
|
|
// follows the noise in g' instead, and it does so confidently.
|
|
//
|
|
// `mean` is a per-pixel projection over a few tens of frames, NAN where no frame contributed.
|
|
std::optional<BeamCenterEstimate>
|
|
FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const PixelMask &mask,
|
|
const std::vector<float> &mean);
|