leonarski_fandClaude Opus 5
46c50ac744
WriteReflections had no test at all, and the defect it hid was invisible to two of the three libraries that read an MTZ: with the data dataset written at id 0 - the id reserved for HKL_base - gemmi and iotbx still returned the real wavelength, and only CCP4's mtzinfo fell back to the 1.54187 A Cu K-alpha default. iotbx reading such a file reports two datasets both numbered 0. The test writes a merged MTZ through the real writer at a wavelength well away from that default, reads the file back, and asserts the layout the report depends on: HKL_base at id 0, the data dataset at id 1 carrying the wavelength, and every data column owned by the data dataset rather than the base - that ownership is what the reported wavelength is read off. Cell and space group are checked with it, since they travel in the same header. Verified against the pre-fix writer (base absent, H K L added on the data dataset): it fails on the dataset count. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
73 lines
3.1 KiB
C++
73 lines
3.1 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
|
|
#include <filesystem>
|
|
#include <string>
|
|
|
|
#include <gemmi/mtz.hpp>
|
|
|
|
#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<int>(reflections.size()));
|
|
}
|