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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu
(cherry picked from commit dafb8d1bf996d1000741eecf019aee780dd91ebe)

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<int16_t> 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<SpotToSave> MakeTestSpots(int i) {
     return {
         SpotToSave{

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<double> unit(0.0, 1.0);
+    std::exponential_distribution<double> wilson(1.0);
+    std::normal_distribution<double> gauss(0.0, 1.0);
+    std::vector<IntegrationOutcome> 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<float>(d);
+            r.dist_ewald = static_cast<float>(std::fabs(de));
+            r.I = static_cast<float>(*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<float>(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") {