leonarski_fandClaude Opus 5
1388b16d7a
The spot-finding resolution estimate was clamped so it could never beat the detector corner. On a crystal that diffracts past the corner that reports where the DETECTOR stops, which is the one thing this number is not for - it is meant to say how far a merge of data like these would reach, a property of the crystal and the exposure. The clamp also hid the interesting case: an estimate finer than what the run actually merged is the statement "this run was detector-limited", and there was no way to make it. The statistic already extrapolates. Its quantile sits in the middle of the fall-off, well inside what the detector records, so it goes on measuring the crystal's own decay when the detector cuts that decay short. Measured by truncating the spot lists of 31 battery crystals at an artificial detector edge and scoring the unclamped answer against each crystal's own measured CC1/2 = 0.30 crossing, it holds its 8-9% floor out to about 1.7x past the cut and only then drifts pessimistic, which is the safe direction. Every genuinely detector-limited crystal in the battery needs between 1.10x and 1.63x. Against a truth corrected for censoring - the six crystals whose merge is cut off by their own detector cannot have a measured crossing, so theirs is extrapolated from multiplicity-corrected <I/sigma> and anchored on the 25 where both exist: symmetric-log RMS 13.5 -> 9.4% over 37 crystals, 25 -> 28 within 0.2 A. On the six detector-limited ones 26.1 -> 9.8% and the bias goes +19 -> -3%; on the 31 that are not, 9.23 -> 9.32%, i.e. it costs them nothing. The 0.30 tail fraction and the 2.25 reach were refit by leave-one-out against that truth and did not move. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01FBumeJVx4oeXxiBRpkrE5H
62 lines
3.1 KiB
C++
62 lines
3.1 KiB
C++
// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include "../image_analysis/spot_finding/SpotUtils.h"
|
|
|
|
TEST_CASE("FilterSpuriousHighResolutionSpots") {
|
|
std::vector<SpotToSave> spots;
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 18.0, .indexed = false});
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 20.0, .indexed = false});
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 30.0, .indexed = false});
|
|
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 6.0, .indexed = false});
|
|
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 2.0, .indexed = false});
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 1.9, .indexed = false});
|
|
spots.push_back(SpotToSave{.x = 1, .y = 2, .intensity = 3, .d_A = 1.3, .indexed = false});
|
|
|
|
FilterSpuriousHighResolutionSpots(spots, 1.57); // roughly 0.25 in 1/d
|
|
|
|
REQUIRE(spots.size() == 4);
|
|
// Spots are sorted by resolution
|
|
CHECK(spots[0].d_A == Catch::Approx(30.0));
|
|
CHECK(spots[1].d_A == Catch::Approx(20.0));
|
|
CHECK(spots[2].d_A == Catch::Approx(18.0));
|
|
CHECK(spots[3].d_A == Catch::Approx(6.0));
|
|
}
|
|
|
|
|
|
TEST_CASE("GetResolution") {
|
|
// Eleven equally strong spots at 1/d^2 = 0.1, 0.2, ... 1.1. Walking in from the highest-resolution
|
|
// one, four of the eleven are the first to carry 30% of the weight, so the quantile is the fourth
|
|
// spot in, 1/d^2 = 0.8. The estimate is that resolution taken 2.25x further in 1/d.
|
|
std::vector<SpotToSave> spots;
|
|
for (int i = 1; i <= 11; i++)
|
|
spots.push_back(SpotToSave{.intensity = 100.0f, .d_A = 1.0f / std::sqrt(0.1f * static_cast<float>(i))});
|
|
const auto d = GetResolution(spots);
|
|
REQUIRE(d.has_value());
|
|
CHECK(*d == Catch::Approx(1.0 / (2.25 * std::sqrt(0.8))).epsilon(1e-4));
|
|
|
|
// The answer is not limited to what a detector records. Keeping only the five spots a detector
|
|
// reaching 1/d^2 = 0.5 would have recorded leaves the quantile at 0.4, and the estimate still
|
|
// extrapolates 2.25x past it instead of stopping at the cut.
|
|
const std::vector<SpotToSave> cut(spots.begin(), spots.begin() + 5);
|
|
CHECK(*GetResolution(cut) == Catch::Approx(1.0 / (2.25 * std::sqrt(0.4))).epsilon(1e-4));
|
|
|
|
// Ice-flagged spots take no part, however strong they are.
|
|
std::vector<SpotToSave> with_ice = spots;
|
|
with_ice.push_back(SpotToSave{.intensity = 1e6f, .d_A = 0.5f, .ice_ring = true});
|
|
CHECK(*GetResolution(with_ice) == Catch::Approx(*d));
|
|
|
|
// A weak high-resolution spot moves the answer far less than a strong one, which is the point of
|
|
// weighting by sqrt(I) rather than counting: the old order statistic would follow it entirely.
|
|
std::vector<SpotToSave> with_spur = spots;
|
|
with_spur.push_back(SpotToSave{.intensity = 1.0f, .d_A = 0.5f});
|
|
CHECK(*GetResolution(with_spur) == Catch::Approx(*d).epsilon(0.02));
|
|
|
|
// Too few spots to have a fall-off at all.
|
|
CHECK_FALSE(GetResolution(std::vector<SpotToSave>(3)).has_value());
|
|
}
|