From e6b550fe42bc556e57a41a76b02460b4f5382f5a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 26 Oct 2026 11:34:56 +0100 Subject: [PATCH] Tests: the stills rocking width - the window at a reflection's resolution, _process.h5 round trip StillsPartialityRefine_MeasureRunWidth_WindowAtReflectionResolution, JFJochReader_StillsRunWidth. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu (cherry picked from commit dafb8d1bf996d1000741eecf019aee780dd91ebe) --- tests/JFJochReaderTest.cpp | 55 ++++++++++++++++++++++++++++++++++++++ tests/MergeScaleTest.cpp | 42 +++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 0f7b952f1..320d434fe 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -2776,6 +2776,61 @@ TEST_CASE("JFJochReader_ReadReflections_Reindex", "[HDF5][Full]") { REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } +// A stills run's rocking curve (the one its partialities were computed from) goes to the master file and +// comes back whole, so --mode scale merges with it instead of measuring another; a file without it reads +// back without one. +TEST_CASE("JFJochReader_StillsRunWidth", "[HDF5][Full]") { + for (const bool with_width: {true, false}) { + DiffractionExperiment x(DetJF(1)); + x.FilePrefix("read_stills_run_width") + .ImagesPerTrigger(1) + .ImagesPerFile(1) + .OverwriteExistingFiles(true) + .BitDepthImage(16) + .PixelSigned(true) + .SetFileWriterFormat(FileWriterFormat::NXmxVDS) + .Compression(CompressionAlgorithm::NO_COMPRESSION); + + const StillsRunWidth w{.gamma0 = 1.25e-4, .gamma_e = 6.1e-4, .bandwidth_sigma = 1.7e-3, + .radial_scale = 1.0 + 3.3e-4, + .table_dstar = {0.1, 0.25, 0.4}, .table_s = {1.1e-4, 1.4e-4, 2.2e-4}}; + RegisterHDF5Filter(); + { + StartMessage start_message; + x.FillMessage(start_message); + FileWriter writer(start_message); + std::vector image(x.GetPixelsNum(), 0); + DataMessage message{}; + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = 0; + REQUIRE_NOTHROW(writer.WriteHDF5(message)); + EndMessage end_message; + end_message.max_image_number = 1; + if (with_width) + end_message.stills_run_width = w; + writer.WriteHDF5(end_message); + writer.Finalize(); + } + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("read_stills_run_width_master.h5")); + const auto &got = reader.GetDataset()->stills_run_width; + REQUIRE(got.has_value() == with_width); + if (with_width) { + CHECK(got->gamma0 == w.gamma0); + CHECK(got->gamma_e == w.gamma_e); + CHECK(got->bandwidth_sigma == w.bandwidth_sigma); + CHECK(got->radial_scale == w.radial_scale); + CHECK(got->table_dstar == w.table_dstar); + CHECK(got->table_s == w.table_s); + } + } + remove("read_stills_run_width_master.h5"); + remove("read_stills_run_width_data_000001.h5"); + } + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + static std::vector MakeTestSpots(int i) { return { SpotToSave{ diff --git a/tests/MergeScaleTest.cpp b/tests/MergeScaleTest.cpp index 32301970c..b7260c01f 100644 --- a/tests/MergeScaleTest.cpp +++ b/tests/MergeScaleTest.cpp @@ -1078,6 +1078,48 @@ TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WideWindowFramesIgnored") { CHECK(width->gamma_e == Catch::Approx(GAMMA_E).epsilon(0.15)); } +// Frames predicted through the run's width (2 sigma(d*) of its formula) have a window at every reflection's +// own resolution; measured with that window, the table gives the width back at each shell. +TEST_CASE("StillsPartialityRefine_MeasureRunWidth_WindowAtReflectionResolution") { + constexpr double GAMMA0 = 1.0e-4, GAMMA_E = 0.035 * M_PI / 180.0; + DiffractionExperiment x; + x.SetSpaceGroup(gemmi::get_spacegroup_by_name("P 1")); + const StillsPartialityRefine::RunWidth window{.gamma0 = GAMMA0, .gamma_e = GAMMA_E}; + auto sigma = [&](double dstar) { return std::sqrt(GAMMA0 * GAMMA0 + GAMMA_E * GAMMA_E * dstar * dstar); }; + + std::mt19937 rng(11); + std::uniform_real_distribution unit(0.0, 1.0); + std::exponential_distribution wilson(1.0); + std::normal_distribution gauss(0.0, 1.0); + std::vector outcomes(400); + for (auto &o: outcomes) { + o.image_scale_g = 0.5 + unit(rng); + for (int j = 0; j < 200; ++j) { + const double d = 1.5 + 8.0 * unit(rng); + const double s = sigma(1.0 / d); + const double de = 2.0 * s * (2.0 * unit(rng) - 1.0); + Reflection r{}; + r.h = j; r.k = 1; r.l = 1; + r.d = static_cast(d); + r.dist_ewald = static_cast(std::fabs(de)); + r.I = static_cast(*o.image_scale_g * 1000.0 * wilson(rng) * std::exp(-0.5 * de * de / (s * s)) + + 5.0 * gauss(rng)); + r.sigma = 5.0f; + r.prescaling_corr = 1.0f; + r.partiality = 1.0f; + r.image_scale_corr = static_cast(1.0 / *o.image_scale_g); + o.reflections.push_back(r); + } + } + const auto width = StillsPartialityRefine(x).MeasureRunWidth(outcomes, std::nullopt, 1.0, &window); + REQUIRE(width.has_value()); + REQUIRE(width->table_s.size() >= 3); + for (size_t i = 0; i < width->table_s.size(); ++i) { + INFO("d* = " << width->table_dstar[i]); + CHECK(width->table_s[i] == Catch::Approx(sigma(width->table_dstar[i])).epsilon(0.15)); + } +} + // Outlier rejection on the stills merge: one observation at ten times its reflection's intensity among // nineteen consistent ones is dropped, and nothing else is. TEST_CASE("MergeOnTheFly_RejectsAnOutlierInTheBand") {