diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index c08a7617e..0912de2a4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1076,9 +1076,19 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru const float r2 = r1 + 2.0f; const float r3 = std::sqrt(r2 * r2 + 133.0f); BraggIntegrationSettings bis = experiment_.GetBraggIntegrationSettings(); - // Remembered only where the radius actually moved, so the guard in RunAllPasses has + // Remembered where the APERTURE moved, so the starvation guard in RunAllPasses has // something to give back and cannot announce a fallback that changes nothing. - if (r1 != bis.GetR1()) bragg_before_adaptive_ = bis; + // + // The test is on all three radii and not on r1 alone. It is the same test today - r2 and + // r3 are taken from r1 just above, and the adaptive rule does not run at all when the + // radii were given explicitly - so this changes no run. It is written this way because + // what the guard is about is the aperture, and r1 is only a proxy for it. A prototype + // that widened r2 and r3 at high angle while leaving r1 at the shipped 4 took the + // neighbour-starved fraction of a dense low-symmetry crystal from 2.4% to 9.0%, eight + // times the bound this guard exists to enforce, and the guard could not fire: nothing + // had set bragg_before_adaptive_, because r1 had not moved. + if (r1 != bis.GetR1() || r2 != bis.GetR2() || r3 != bis.GetR3()) + bragg_before_adaptive_ = bis; bis.R1(r1).R2(r2).R3(r3); experiment_.ImportBraggIntegrationSettings(bis); logger.Info("Spot width: r80 = {:.2f} px at {:.0f} A ({} spots from {} images) => " diff --git a/rugnux/SpotWidth.cpp b/rugnux/SpotWidth.cpp index c31d1ccde..9095f6dbf 100644 --- a/rugnux/SpotWidth.cpp +++ b/rugnux/SpotWidth.cpp @@ -272,7 +272,11 @@ void MeasureSpotFluxCurves(const ImagePreprocessorBuffer &image, int width, int } } -std::optional spot_width::R80AtReference(const std::vector &curves) { +float spot_width::R80Fit::At(double d_A) const { + return static_cast(std::clamp(c0 + c1 / d_A, lo, hi)); +} + +std::optional spot_width::FitR80(const std::vector &curves) { if (curves.size() < MIN_SPOTS_TOTAL) return std::nullopt; // One point per band: the median curve of the band, the radius it holds 80 % of its flux at, and @@ -302,7 +306,18 @@ std::optional spot_width::R80AtReference(const std::vector &cu static_cast(band_d.size())}); } if (points.empty()) return std::nullopt; - if (points.size() == 1) return static_cast(points[0].r80); + + // Never extrapolate outside what the bands actually measured. A single band measures no slope, so + // its own value is the whole law; the bounds then bracket it and At() returns it unchanged. + R80Fit fit; + double lo = std::numeric_limits::max(), hi = 0.0; + for (const auto &p : points) { lo = std::min(lo, p.r80); hi = std::max(hi, p.r80); } + fit.lo = 0.8 * lo; + fit.hi = 1.25 * hi; + if (points.size() == 1) { + fit.c0 = points[0].r80; + return fit; + } // The mosaic contribution to the detector footprint grows as 1/d, so r80 is linear in 1/d. double sw = 0.0, sx = 0.0, sxx = 0.0, sy = 0.0, sxy = 0.0; @@ -314,15 +329,18 @@ std::optional spot_width::R80AtReference(const std::vector &cu sxy += p.weight * p.inv_d * p.r80; } const double det = sw * sxx - sx * sx; - double value = sy / sw; + fit.c0 = sy / sw; if (std::abs(det) > 1e-12) { - const double c1 = (sw * sxy - sx * sy) / det; - value = (sy - c1 * sx) / sw + c1 / D_REF_A; + fit.c1 = (sw * sxy - sx * sy) / det; + fit.c0 = (sy - fit.c1 * sx) / sw; } - // Never extrapolate outside what the bands actually measured. - double lo = std::numeric_limits::max(), hi = 0.0; - for (const auto &p : points) { lo = std::min(lo, p.r80); hi = std::max(hi, p.r80); } - return static_cast(std::clamp(value, 0.8 * lo, 1.25 * hi)); + return fit; +} + +std::optional spot_width::R80AtReference(const std::vector &curves) { + const auto fit = FitR80(curves); + if (!fit) return std::nullopt; + return fit->At(D_REF_A); } float spot_width::R1ForWidth(float r80) { diff --git a/rugnux/SpotWidth.h b/rugnux/SpotWidth.h index e5b450c7e..00fdb3d60 100644 --- a/rugnux/SpotWidth.h +++ b/rugnux/SpotWidth.h @@ -52,9 +52,28 @@ struct FluxCurve { std::array c{}; }; -// r80 at D_REF_A from a pool of curves, or nothing where too few spots were measurable. Fitted as -// r80(d) = c0 + c1/d over resolution bands, weighted by the spots in each, so crystals whose strong -// spots sit at different resolutions are compared like with like. +// The fitted width law: r80(d) = c0 + c1/d over resolution bands, weighted by the spots in each, so +// crystals whose strong spots sit at different resolutions are compared like with like. The fit is +// over 1/d because the mosaic contribution to the detector footprint grows that way, so c1 is the +// crystal's own answer to how much wider its spots get at high angle - a quantity the fit has always +// computed and nothing has ever been able to read, because only the value at one resolution was +// returned. +struct R80Fit { + // Kept in double, and evaluated in double, because At() has to reproduce the single number the + // caller used to get bit for bit: rounding the coefficients to float first and adding them after + // moves the last bit, and r1 = round(2*r80) can turn on it. + double c0 = 0.0, c1 = 0.0; + // What the bands actually measured, widened the way this fit has always widened it. The law is + // HELD at these outside the measured range rather than extrapolated - it has two coefficients and + // as few as two points, so a slope read off a short lever arm must not be allowed to run away. + double lo = 0.0, hi = 0.0; + [[nodiscard]] float At(double d_A) const; +}; + +// The width law from a pool of curves, or nothing where too few spots were measurable. +[[nodiscard]] std::optional FitR80(const std::vector &curves); + +// r80 at D_REF_A, which is the law above read at one resolution. [[nodiscard]] std::optional R80AtReference(const std::vector &curves); // The signal radius the measured width asks for: r1 = clamp(round(2*r80), 4, 6). The floor is the