diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d652cc312..a4af137b6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -96,6 +96,7 @@ ADD_EXECUTABLE(jfjoch_test FrenchWilsonTest.cpp ReindexAmbiguityTest.cpp LoadReferenceMtzFreeFlagsTest.cpp + WriteReflectionsTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp new file mode 100644 index 000000000..1508145f3 --- /dev/null +++ b/tests/WriteReflectionsTest.cpp @@ -0,0 +1,72 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include + +#include + +#include "../common/DiffractionExperiment.h" +#include "../image_analysis/WriteReflections.h" +#include "SyntheticMergedReflections.h" + +namespace { + // The wavelength CCP4's mtzlib substitutes when the data columns belong to dataset 0, which it + // takes for the reserved HKL_base: Cu K-alpha, and the wrong edge for anything reading f'/f'' + // out of the file. + constexpr double CU_KALPHA_A = 1.54187; + + constexpr UnitCell TETRAGONAL_CELL{47.0f, 47.0f, 63.0f, 90.0f, 90.0f, 90.0f}; + + DiffractionExperiment TestExperiment() { + DiffractionExperiment x; + x.IncidentEnergy_keV(12.7f); // ~0.976 A, nowhere near the Cu K-alpha default + x.SpaceGroupNumber(96); // P 43 21 2 + x.SetUnitCell(TETRAGONAL_CELL); + return x; + } +} + +TEST_CASE("Merged MTZ: the data dataset is id 1 and carries the wavelength", "[write_reflections]") { + jfjoch_test::SyntheticMergeParams params; + params.true_space_group = "P 43 21 2"; + params.twin_supergroup = "P 43 21 2"; + params.d_min_A = 5.0; + const auto reflections = jfjoch_test::GenerateSyntheticMerged(params); + REQUIRE(!reflections.empty()); + + const auto experiment = TestExperiment(); + const auto path = (std::filesystem::temp_directory_path() / "rugnux_merged_wavelength.mtz").string(); + WriteMtzReflections(reflections, TETRAGONAL_CELL, experiment, path); + const gemmi::Mtz mtz = gemmi::read_mtz_file(path); + std::filesystem::remove(path); + + // HKL_base at id 0, the data at id 1. A data dataset written at id 0 occupies the id MTZ + // reserves for the base, and mtzlib then reports CU_KALPHA_A instead of the real wavelength. + REQUIRE(mtz.datasets.size() == 2); + CHECK(mtz.datasets[0].id == 0); + CHECK(mtz.datasets[0].dataset_name == "HKL_base"); + CHECK(mtz.datasets[1].id == 1); + CHECK(mtz.datasets[1].wavelength == Catch::Approx(experiment.GetWavelength_A()).epsilon(1e-5)); + CHECK(mtz.datasets[1].wavelength != Catch::Approx(CU_KALPHA_A).epsilon(1e-3)); + + // The wavelength is read off the dataset the data columns belong to, so they have to be on the + // data dataset and not on the base. + for (const char *label : {"IMEAN", "SIGIMEAN", "F", "SIGF", "FreeR_flag"}) { + const gemmi::Mtz::Column *col = mtz.column_with_label(label); + REQUIRE(col != nullptr); + CHECK(col->dataset_id == 1); + } + + // Cell and space group travel in the same header. + CHECK(mtz.spacegroup != nullptr); + CHECK(mtz.spacegroup->number == 96); + CHECK(mtz.cell.a == Catch::Approx(TETRAGONAL_CELL.a).epsilon(1e-5)); + CHECK(mtz.cell.c == Catch::Approx(TETRAGONAL_CELL.c).epsilon(1e-5)); + CHECK(mtz.cell.gamma == Catch::Approx(90.0).epsilon(1e-5)); + CHECK(mtz.datasets[1].cell.a == Catch::Approx(TETRAGONAL_CELL.a).epsilon(1e-5)); + + CHECK(mtz.nreflections == static_cast(reflections.size())); +}