From 5c44e544dc219561dd4c0b0d610e120ad366331b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 30 Aug 2026 12:19:26 +0200 Subject: [PATCH 1/2] reader: place a detector swung out on a 2theta arm where the file says it stands Chemical crystallography reaches high angle by swinging the detector out on a 2theta arm. Both readers had the number and neither used it: the miniCBF header's Detector_2theta was parsed into a struct member nothing ever read, and on the NXmx side the rotation was in the depends_on chain, which was not followed at all. A sweep taken at 30 degrees was therefore processed with its detector plane 30 degrees from where it stood, and nothing indexed. The geometry could already express it, and needed no change: the arm turns the detector about the sample, so the distance is still measured along the detector normal and the beam centre is still the point of normal incidence - which is exactly the PONI convention, and a swung detector is one PONI rotation. What moves is the direct beam, by distance*tan(2theta), off the beam centre and often off the detector. NXmx is the harder half, because the swing has no field of its own: it is one rotation in the chain the detector's position depends on, and "two_theta" is only one beamline's name for that dataset. So the chain is followed and its rotations composed, rather than a field of one name being looked for - each transformation states its vector in the frame of the one it depends on, which is why the product is the whole placement. Translations are skipped; they are the distance and the beam centre, which the file states separately in the square-on frame. Vectors come from McStas through the same 180-degree turn about z the module directions already use, a proper rotation, so an axis carried through it turns the same way. The three rotations a file this system writes ARE that chain, and are also read as the PONI angles - so those three paths are skipped, or every tilted file we have ever written would come back tilted twice. That is the one way this change could have broken existing data, and the test for it writes a tilted file and reads it back. For miniCBF the arm turns about the base spindle axis: on the four-circle geometry those headers describe the two are one axis, and the imgCIF axis table such a header carries states them with the same vector. Both now come from one constant, so a later correction to the frame moves them together. Measured. On a swung NXmx sweep the chain gives rot2 = -0.34907 rad for the 20 degrees it states, and the sweep goes from "nothing was integrated" to 25000 reflections at 82.2% completeness and CC(1/2) 0.9993, in the same space group and the same cell to 0.03 A as the square-on sweep of that crystal; the opposite sign indexes nothing. A miniCBF sweep at 30 degrees goes the same way, to 0.585 A, and a second sweep of that crystal at 55 degrees reaches 0.476 A and reproduces the cell again - with a low-resolution limit of 2.36 A rather than 13 A, which is what a detector swung that far records. On all of them post-refinement recovers the header's own beam centre and distance, and the beam stop shadow sits within four pixels of where the swung geometry puts the direct beam, 417 and 537 pixels from where the unswung one does. Seven sets whose detector is square to the beam, three of them carrying a chain whose 2theta is zero, are byte-identical in .hkl, .mtz, .cif and the image statistics. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f --- docs/CHANGELOG.md | 1 + docs/DETECTOR_GEOMETRY.md | 15 +++ reader/HDF5MetadataSource.cpp | 84 +++++++++++++- reader/JFJochCBFReader.cpp | 23 +++- tests/DiffractionGeometryTest.cpp | 45 ++++++++ tests/JFJochReaderTest.cpp | 179 ++++++++++++++++++++++++++++++ 6 files changed, 343 insertions(+), 4 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index e2c0875ec..862358d7a 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,7 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* A detector swung out on a 2theta arm is placed where it stands, from the depends_on transformation chain of an NXmx master or the `Detector_2theta` line of a miniCBF header; both were previously read and then ignored. * `jfjoch_viewer` opens PILATUS miniCBF sweeps - naming any frame opens the whole sweep - and can run a processing job on one. * A detector whose stored image is mirrored in Y or mounted at a multiple of 90 degrees can be described as such, in the detector configuration or with `--detector-mirror-y` / `--detector-quarter-turns`, rather than having to be expressed as a detector rotation. * The rotation first pass refines twelve candidate lattices rather than four, so a correct cell that the pre-refinement ranking put fifth is still reached. diff --git a/docs/DETECTOR_GEOMETRY.md b/docs/DETECTOR_GEOMETRY.md index c59b41e77..3cc216692 100644 --- a/docs/DETECTOR_GEOMETRY.md +++ b/docs/DETECTOR_GEOMETRY.md @@ -59,6 +59,21 @@ rot2 = asin(-slow.z) rot1 = atan2(-fast.z, normal.z) rot3 = atan with `rot2` in [-90°, 90°]. The angles are what is stored and what is written out, so a geometry given as angles comes back exactly as it was given. +## A detector swung out on a 2theta arm + +Chemical crystallography reaches high angle by swinging the detector out on a 2theta arm rather than by +moving it closer. The arm turns the detector about the sample, so it changes nothing else: the distance +is still measured along the detector normal, and the beam centre is still the point of normal incidence, +which is where the arm's own axis meets the detector and does not move. The swing is therefore exactly a +PONI rotation, and the direct beam is what moves - by `distance * tan(2theta)`, off the beam centre and +often off the detector altogether. + +Nothing has to be given for this: rugnux takes it from the file. An NXmx master states the detector's +position as a `depends_on` chain of transformations, and the arm is one rotation in that chain - so the +chain is followed, rather than a field of one particular name being looked for. A PILATUS miniCBF states +it as `# Detector_2theta`, which turns about the same axis as the base spindle, the two being one axis on +the four-circle geometry those headers describe. + ## Mirrored and quarter-turned detectors On top of the continuous tilt the detector setup carries a **discrete image orientation**: whether the diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 47686e91b..ac68aaa60 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -14,8 +14,14 @@ #include "../common/Logger.h" #include "../common/ROIDefinition.h" +// A McStas direction in the internal frame. The two differ by a 180 degree turn about z, which is a +// rotation and not a mirror - so an axis carried through it turns the same way by the same angle. +static Coord McStasToInternal(const std::vector &v) { + return {static_cast(-v[0]), static_cast(-v[1]), static_cast(v[2])}; +} + // The image orientation the file itself states, in its NXdetector_module pixel directions. NXmx gives -// those in the McStas frame, which is the internal frame turned 180 degrees about z. +// those in the McStas frame. // // Only an exact match against one of the eight discrete orientations is taken. Anything else is a // continuous rotation of the detector in its own plane, which belongs in rot1/rot2/rot3 and cannot be @@ -36,8 +42,8 @@ static std::optional ReadModuleOrientation(HDF5Object *file if ((f.size() != 3) || (s.size() != 3)) return {}; - const Coord fast(-f[0], -f[1], f[2]); - const Coord slow(-s[0], -s[1], s[2]); + const Coord fast = McStasToInternal(f); + const Coord slow = McStasToInternal(s); for (int64_t quarter_turns = 0; quarter_turns < 4; quarter_turns++) { for (bool mirror_y: {false, true}) { @@ -50,6 +56,63 @@ static std::optional ReadModuleOrientation(HDF5Object *file return {}; } +// Where the detector stands, from the chain of transformations the file says it depends on. +// +// NXmx has no field for a detector swung out on a 2theta arm. It states the detector's position as a +// depends_on chain and the arm is one rotation in that chain, so following the chain is the only way +// to find it: "two_theta" is one beamline's name for that dataset and the next spells it otherwise. +// +// Only the rotations are taken, composed from the detector outwards. Each transformation states its +// vector in the frame of the one it depends on, so the product is the rotation that carries a +// detector square to the beam to where this one stands. The translations in the chain are the +// detector distance and the beam centre, which the file states separately in that square-on frame - +// the arm turns the detector about the sample and moves neither, and a Diamond master writes the same +// beam_center_x/y for a swung sweep as for the square-on one beside it. Nothing comes back when no +// rotation in the chain turns, which is every detector square to the beam. +static std::optional ReadDetectorRotationChain(HDF5Object *file) { + std::string node = file->GetString("/entry/instrument/detector/depends_on"); + if (node.empty() && file->IsDataSet("/entry/instrument/detector/module/module_offset")) { + HDF5DataSet module_offset(*file, "/entry/instrument/detector/module/module_offset"); + if (module_offset.AttrExists("depends_on")) + node = module_offset.ReadAttrStr("depends_on"); + } + + // A file this system wrote states its PONI angles in the chain as well, and they are read from + // these three paths just before this is called. Taking them here too would apply the tilt twice. + static const std::set poni_angles = {"/entry/instrument/detector/transformations/rot1", + "/entry/instrument/detector/transformations/rot2", + "/entry/instrument/detector/transformations/rot3"}; + + RotMatrix chain; + bool turns = false; + std::set seen; + while ((node != ".") && !node.empty() && file->IsDataSet(node) && seen.insert(node).second) { + HDF5DataSet axis(*file, node); + const std::string current = node; + node = axis.AttrExists("depends_on") ? axis.ReadAttrStr("depends_on") : "."; + + if (poni_angles.contains(current) || !axis.AttrExists("transformation_type") || !axis.AttrExists("vector") + || (axis.ReadAttrStr("transformation_type") != "rotation")) + continue; + + std::vector value; + axis.ReadVector(value); + const auto vec = axis.ReadAttrVec("vector"); + if (value.empty() || (value[0] == 0.0) || (vec.size() != 3)) + continue; + + // NXmx states a rotation in degrees unless it says otherwise. + const bool radians = axis.AttrExists("units") && (axis.ReadAttrStr("units") == "rad"); + const auto angle_rad = static_cast(radians ? value[0] : value[0] * PI / 180.0); + chain = RotMatrix(angle_rad, McStasToInternal(vec)) * chain; + turns = true; + } + + if (!turns) + return {}; + return chain; +} + inline std::pair parse_bravais_lattice(const std::string &val) { if (val.empty()) @@ -721,6 +784,21 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen master_file->GetOptFloat("/entry/instrument/detector/transformations/rot2").value_or(0.0)); dataset->experiment.PoniRot3_rad( master_file->GetOptFloat("/entry/instrument/detector/transformations/rot3").value_or(0.0)); + // A detector swung out on a 2theta arm - routine in chemical crystallography - and any other + // rotation the file puts in the detector's chain. It turns the detector about the sample, so + // it carries the whole square-on geometry with it and composes on the left of the PONI + // rotation the file states directly. + if (const auto chain = ReadDetectorRotationChain(master_file.get())) { + float rot1 = 0, rot2 = 0, rot3 = 0; + PoniAnglesFromMatrix(chain.value() + * PoniRotMatrix(dataset->experiment.GetPoniRot1_rad(), + dataset->experiment.GetPoniRot2_rad(), + dataset->experiment.GetPoniRot3_rad()), + rot1, rot2, rot3); + dataset->experiment.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3); + Logger("HDF5Reader").Info("Detector placed by its NXmx transformation chain: " + "rot1 {:.5f} rot2 {:.5f} rot3 {:.5f} rad", rot1, rot2, rot3); + } dataset->experiment.SampleTemperature_K(master_file->GetOptFloat("/entry/sample/temperature")); dataset->experiment.BeamX_pxl(master_file->GetFloat("/entry/instrument/detector/beam_center_x")); diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index 1bd4e038e..a08639663 100644 --- a/reader/JFJochCBFReader.cpp +++ b/reader/JFJochCBFReader.cpp @@ -90,6 +90,11 @@ std::vector CollectSweep(const std::string &path) { return out; } +// The base rotation axis of the instrument, in the internal frame (x along increasing detector +// column, y along increasing row, z along the beam). Both the spindle and the detector arm turn +// about it: the spindle in RotationAxis below, and the 2theta arm in ReadFiles. +const Coord BASE_AXIS(-1.0f, 0.0f, 0.0f); + // The axis a miniCBF sweep turns about, in the internal frame (x along increasing detector column, // y along increasing row, z along the beam). // @@ -106,7 +111,7 @@ std::vector CollectSweep(const std::string &path) { // source; internal z points the other way, hence the minus. A kappa arm cannot be expressed at all: // its inclination is a property of the hardware that no miniCBF header states. Coord RotationAxis(const minicbf::Header &h) { - const Coord base(-1.0f, 0.0f, 0.0f); + const Coord base = BASE_AXIS; if (!minicbf::ScansPhi(h)) return base; @@ -161,6 +166,22 @@ void JFJochCBFReader::ReadFiles(const std::string &path) { dataset_->experiment.BeamX_pxl(static_cast(header0_.beam_x_px)); dataset_->experiment.BeamY_pxl(static_cast(header0_.beam_y_px)); dataset_->experiment.DetectorDistance_mm(static_cast(header0_.distance_m * 1000.0)); + + // A detector swung out on a 2theta arm, which small-molecule collection uses routinely. The arm + // turns the detector about the sample, so it carries the square-on geometry with it: the header's + // Detector_distance stays the distance along the detector normal and Beam_xy stays the point of + // normal incidence, neither of which the swing moves - which is exactly what the PONI convention + // wants, so the swing is a PONI rotation and nothing else in the header changes. It turns about + // the base spindle axis, the four-circle geometry these headers describe having the arm and the + // spindle on one axis; the imgCIF axis table such a header carries states the two with the same + // vector. + if (header0_.two_theta_deg != 0.0) { + float rot1 = 0, rot2 = 0, rot3 = 0; + PoniAnglesFromMatrix(RotMatrix(static_cast(header0_.two_theta_deg * PI / 180.0), BASE_AXIS), + rot1, rot2, rot3); + dataset_->experiment.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3); + } + dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast(header0_.wavelength_A)); dataset_->experiment.FrameTime( std::chrono::duration_cast( diff --git a/tests/DiffractionGeometryTest.cpp b/tests/DiffractionGeometryTest.cpp index 81f7651f5..690807168 100644 --- a/tests/DiffractionGeometryTest.cpp +++ b/tests/DiffractionGeometryTest.cpp @@ -831,3 +831,48 @@ TEST_CASE("DetectorOrientation_recip_roundtrip") { } } } + +// A detector swung out on a 2theta arm, which is how chemical crystallography reaches high angle. +// The arm turns the detector about the sample, so the geometry that describes it is the PONI rotation +// and nothing else moves: the distance stays the distance along the detector normal and the beam +// centre stays the point of normal incidence. What DOES move is the direct beam, which is no longer +// at the beam centre - the two coincide only on a detector square to the beam. +TEST_CASE("DiffractionGeometry_TwoThetaArm", "[LinearAlgebra][Coord]") { + const float two_theta = 30.0f * PI / 180.0f; + const float distance_mm = 160.0f, pixel_mm = 0.172f, wavelength = 0.6889f; + const float bx = 740.0f, by = 866.0f; + + DiffractionGeometry geom; + geom.BeamX_pxl(bx).BeamY_pxl(by).DetectorDistance_mm(distance_mm) + .PixelSize_mm(pixel_mm).Wavelength_A(wavelength); + // The arm turns about the internal x axis; a rotation of +2theta about it is rot2 = -2theta. + geom.PoniRot2_rad(-two_theta); + + // The beam centre pixel is the PONI: still on the detector normal through the sample, and now + // 2theta away from the beam. + CHECK(geom.TwoTheta_rad(bx, by) == Catch::Approx(two_theta)); + CHECK(geom.LabCoord(bx, by).Length() == Catch::Approx(distance_mm)); + CHECK(geom.GetNormalAxis() * Coord(0, 0, 1) == Catch::Approx(cosf(two_theta))); + // The plane turned about x, so the fast axis - along +x - did not move, and the slow one tipped + // out of the detector plane by the full 2theta. + CHECK((geom.GetFastAxis() - Coord(1, 0, 0)).Length() < 1e-6f); + CHECK(geom.GetSlowAxis() * Coord(0, 0, 1) == Catch::Approx(sinf(two_theta))); + + // The direct beam is off the PONI by D*tan(2theta), along the direction the arm swung. + auto [direct_x, direct_y] = geom.GetDirectBeam_pxl(); + CHECK(direct_x == Catch::Approx(bx)); + CHECK(direct_y == Catch::Approx(by + distance_mm * tanf(two_theta) / pixel_mm)); + + // Resolution at the PONI is the Bragg spacing of 2theta, not of a pixel at zero distance from + // the beam centre - the reason a swung detector reaches so much further than a square-on one. + CHECK(geom.PxlToRes(bx, by) == Catch::Approx(wavelength / (2.0f * sinf(two_theta / 2.0f)))); + + // Round trip through reciprocal space, at the PONI and away from it in both directions. + const std::vector> probes = + {{bx, by}, {bx + 300.0f, by - 500.0f}, {bx - 700.0f, by + 200.0f}}; + for (const auto &[x, y]: probes) { + auto [back_x, back_y] = geom.RecipToDetector(geom.DetectorToRecip(x, y)); + CHECK(back_x == Catch::Approx(x)); + CHECK(back_y == Catch::Approx(y)); + } +} diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 97410fb40..89c0c3a1c 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -3594,3 +3594,182 @@ TEST_CASE("JFJochReader_ThirdPartyNXmxMaster", "[HDF5][Full]") { // No leftover HDF5 objects REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); } + +// A detector swung out on a 2theta arm. NXmx has no field for it: the swing is one rotation in the +// depends_on chain the detector's position is stated as, and "two_theta" is only one beamline's name +// for that dataset. So the chain is what the reader follows, and the chain here carries two rotations +// about different axes, outboard of the translation that sets the distance - a file that stated only +// the innermost one, or composed them the other way round, gives a different plane. +// +// Both axes are stated in McStas, which is the internal frame turned half a turn about z: a reader +// that takes the vector as it stands swings the detector the wrong way, which is twice the error of +// not reading it at all. +TEST_CASE("JFJochReader_DetectorTwoThetaArm", "[HDF5][Full]") { + const hsize_t nx = 8, ny = 6; + const double two_theta_deg = 20.0, tilt_deg = 7.0; + std::vector image(nx * ny, 5); + WriteThirdPartyDataFile("two_theta_000001.h5", image, 2, ny, nx); + + { + HDF5File master("two_theta_master.h5"); + HDF5Group entry(master, "entry"); + entry.SaveScalar("definition", "NXmx"); + HDF5Group instrument(entry, "instrument"); + HDF5Group beam(instrument, "beam"); + beam.SaveScalar("incident_wavelength", 0.6889)->Units("angstrom"); + HDF5Group transformations(instrument, "transformations"); + // Outermost first in the file, innermost first along the chain: det_z -> two_theta -> tilt + transformations.SaveVector("tilt", std::vector{tilt_deg}) + ->Transformation("deg", ".", "detector", "", "rotation", {0, 1, 0}); + transformations.SaveVector("two_theta", std::vector{two_theta_deg}) + ->Transformation("deg", "/entry/instrument/transformations/tilt", + "detector", "", "rotation", {-1, 0, 0}); + transformations.SaveVector("det_z", std::vector{160.0}) + ->Transformation("mm", "/entry/instrument/transformations/two_theta", + "detector", "", "translation", {0, 0, 1}); + HDF5Group detector(instrument, "detector"); + detector.SaveScalar("depends_on", "/entry/instrument/transformations/det_z"); + detector.SaveScalar("description", "PILATUS 2M"); + detector.SaveScalar("beam_center_x", 4.0)->Units("pixels"); + detector.SaveScalar("beam_center_y", 3.0)->Units("pixels"); + detector.SaveScalar("distance", 0.160)->Units("m"); + detector.SaveScalar("x_pixel_size", 0.172)->Units("mm"); + detector.SaveScalar("y_pixel_size", 0.172)->Units("mm"); + detector.SaveScalar("sensor_thickness", 0.32)->Units("mm"); + detector.SaveScalar("count_time", 0.2); + detector.SaveScalar("saturation_value", static_cast(65535)); + HDF5Group data(entry, "data"); + data.ExternalLink("two_theta_000001.h5", "/data", "data_000001"); + } + + DiffractionGeometry geom; + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("two_theta_master.h5")); + geom = reader.GetDataset()->experiment.GetDiffractionGeometry(); + } + + // The chain as it stands in the internal frame: McStas (-1,0,0) is internal (1,0,0) and McStas + // (0,1,0) is internal (0,-1,0), and the outer rotation multiplies on the left. + const auto to_rad = [](double deg) { return static_cast(deg * PI / 180.0); }; + const RotMatrix expected = RotMatrix(to_rad(tilt_deg), {0, -1, 0}) + * RotMatrix(to_rad(two_theta_deg), {1, 0, 0}); + for (int64_t column = 0; column < 3; column++) + CHECK((geom.GetDetectorMatrix().Column(column) - expected.Column(column)).Length() < 1e-5f); + + // Distance and beam centre are the ones the file states: the arm turns the detector about the + // sample and moves neither. + CHECK(geom.GetDetectorDistance_mm() == Catch::Approx(160.0)); + CHECK(geom.GetBeamX_pxl() == Catch::Approx(4.0)); + CHECK(geom.GetBeamY_pxl() == Catch::Approx(3.0)); + // And the beam centre pixel is now that far from the beam - the whole point of a 2theta arm. + CHECK(geom.TwoTheta_rad(4.0f, 3.0f) * 180.0f / PI + == Catch::Approx(angle_deg(expected * Coord(0, 0, 1), Coord(0, 0, 1)))); + CHECK(geom.TwoTheta_rad(4.0f, 3.0f) * 180.0f / PI > two_theta_deg); + + remove("two_theta_000001.h5"); + remove("two_theta_master.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + +// The same file with the arm parked at zero: the chain is there, nothing in it turns, and the +// geometry must be exactly the square-on one. This is nearly every file, so it has to cost nothing. +TEST_CASE("JFJochReader_DetectorTwoThetaZeroIsSquareOn", "[HDF5][Full]") { + const hsize_t nx = 8, ny = 6; + std::vector image(nx * ny, 5); + WriteThirdPartyDataFile("two_theta_zero_000001.h5", image, 2, ny, nx); + + { + HDF5File master("two_theta_zero_master.h5"); + HDF5Group entry(master, "entry"); + entry.SaveScalar("definition", "NXmx"); + HDF5Group instrument(entry, "instrument"); + HDF5Group beam(instrument, "beam"); + beam.SaveScalar("incident_wavelength", 0.6889)->Units("angstrom"); + HDF5Group transformations(instrument, "transformations"); + transformations.SaveVector("two_theta", std::vector{0.0}) + ->Transformation("deg", ".", "detector", "", "rotation", {-1, 0, 0}); + transformations.SaveVector("det_z", std::vector{160.0}) + ->Transformation("mm", "/entry/instrument/transformations/two_theta", + "detector", "", "translation", {0, 0, 1}); + HDF5Group detector(instrument, "detector"); + detector.SaveScalar("depends_on", "/entry/instrument/transformations/det_z"); + detector.SaveScalar("description", "PILATUS 2M"); + detector.SaveScalar("beam_center_x", 4.0)->Units("pixels"); + detector.SaveScalar("beam_center_y", 3.0)->Units("pixels"); + detector.SaveScalar("distance", 0.160)->Units("m"); + detector.SaveScalar("x_pixel_size", 0.172)->Units("mm"); + detector.SaveScalar("y_pixel_size", 0.172)->Units("mm"); + detector.SaveScalar("sensor_thickness", 0.32)->Units("mm"); + detector.SaveScalar("count_time", 0.2); + detector.SaveScalar("saturation_value", static_cast(65535)); + HDF5Group data(entry, "data"); + data.ExternalLink("two_theta_zero_000001.h5", "/data", "data_000001"); + } + + DiffractionGeometry geom; + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("two_theta_zero_master.h5")); + geom = reader.GetDataset()->experiment.GetDiffractionGeometry(); + } + + CHECK(geom.GetPoniRot1_rad() == 0.0f); + CHECK(geom.GetPoniRot2_rad() == 0.0f); + CHECK(geom.GetPoniRot3_rad() == 0.0f); + CHECK(geom.TwoTheta_rad(4.0f, 3.0f) == 0.0f); + + remove("two_theta_zero_000001.h5"); + remove("two_theta_zero_master.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} + +// A file this system wrote states its own PONI angles twice: as the three scalars the reader takes +// them from, and as three rotations in the detector's depends_on chain. Following the chain must +// therefore skip them - applied on top of the scalars they would tilt the detector twice, which is +// how a correct 2theta reader breaks every tilted file this system has ever written. A test that only +// wrote an untilted detector could not see it. +TEST_CASE("JFJochReader_DetectorChainDoesNotDoubleTheTilt", "[HDF5][Full]") { + const float rot1 = 0.031f, rot2 = -0.047f, rot3 = 0.019f; + + DiffractionExperiment x(DetJF(1)); + x.ImagesPerTrigger(2).OverwriteExistingFiles(true).FilePrefix("test_ponichain"); + x.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(150) + .IncidentEnergy_keV(WVL_1A_IN_KEV).PixelSigned(false).BitDepthImage(16) + .FrameTime(std::chrono::microseconds(500), std::chrono::microseconds(10)); + x.PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3); + + RegisterHDF5Filter(); + std::vector image(x.GetPixelsNum(), 0); + + StartMessage start_message; + x.FillMessage(start_message); + FileWriter file_set(start_message); + DataMessage message{}; + for (int i = 0; i < x.GetImageNum(); i++) { + message.image = CompressedImage(image, x.GetXPixelsNum(), x.GetYPixelsNum()); + message.number = i; + REQUIRE_NOTHROW(file_set.WriteHDF5(message)); + } + EndMessage end_message; + end_message.max_image_number = x.GetImageNum(); + file_set.WriteHDF5(end_message); + file_set.Finalize(); + + DiffractionGeometry geom; + { + JFJochHDF5Reader reader; + REQUIRE_NOTHROW(reader.ReadFile("test_ponichain_master.h5")); + geom = reader.GetDataset()->experiment.GetDiffractionGeometry(); + } + CHECK(geom.GetPoniRot1_rad() == Catch::Approx(rot1).margin(1e-6)); + CHECK(geom.GetPoniRot2_rad() == Catch::Approx(rot2).margin(1e-6)); + CHECK(geom.GetPoniRot3_rad() == Catch::Approx(rot3).margin(1e-6)); + for (int64_t column = 0; column < 3; column++) + CHECK((geom.GetDetectorMatrix().Column(column) + - PoniRotMatrix(rot1, rot2, rot3).Column(column)).Length() < 1e-5f); + + remove("test_ponichain_master.h5"); + remove("test_ponichain_data_000001.h5"); + REQUIRE(H5Fget_obj_count(H5F_OBJ_ALL, H5F_OBJ_ALL) == 0); +} From a27c4cf26fc86c4e9922b8c66535ef541967b1c8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 30 Aug 2026 14:03:32 +0200 Subject: [PATCH 2/2] reader: take a miniCBF's mounting from the imgCIF axis table its header states A miniCBF header states three things about how the instrument is put together that the reader was assuming instead: which laboratory direction the image's columns run along, which its rows run along, and which the spindle turns about. Some beamlines append a CBF template block holding the full imgCIF axis table, which says all three outright. Two instruments in the corpus are not what was assumed, in two different ways. One mounts its detector a quarter turn round, so the image's columns run vertically. Another turns its spindle about the VERTICAL, with the image mounted the usual way; its table says so, and its "# Oscillation_axis" line says so a second way, by naming the image direction the spindle runs along rather than a vector. Either error leaves the spindle 90 degrees from the image. That is not a sign, so the run's axis-sign rescue cannot reach it, and no refinement recovers it: all three affected sweeps indexed nothing usable. So the table is read. The element axes give the image orientation, matched against the eight discrete mountings exactly as the NXmx module directions already are - the match itself moves to DetectorOrientation, so both readers share one definition rather than two copies. The goniometer axis with no parent gives the spindle DIRECTION; its sign stays the rescue's business, which is the part a convention can legitimately differ on. The detector axis with no parent gives the 2theta arm, replacing the assumption that the arm shares the spindle's axis - the one header stating both states them with the same vector, so this changes no answer, only what it rests on. imgCIF's frame differs from the internal one by a half turn about x, a rotation and not a mirror, as writer/HDF5NXmx.cpp already records from the other side. Where a header carries no table, a "+SLOW" on the Oscillation_axis line still says the spindle runs along the image's slow direction. That is the only thing one of the three affected sets says about it. The axis NAME on that line stays unusable - the header that carries both says "X.CW" where its own table says Y - but the direction token is not: where both are present they agree, which is what makes reading it evidence rather than a guess. Also: naming a frame with no directory at all now finds its sweep. parent_path() of a bare filename is empty and iterating an empty path finds nothing, so running from inside the data directory reported that no images were found. Measured, with nothing on the command line. The vertical-spindle protein set goes from no usable lattice to 100% indexed, P 6(3) 2 2 with a cell 0.43% from deposited, 87846 reflections at 86.3% completeness and CC(1/2) 0.995. Its companion from the same detector, which has no table and only the +SLOW token, goes from a spurious monoclinic cell at 2.3% completeness and I/sigma 0.21 to the right orthorhombic lattice, 97.7% indexed, 59.7% complete, CC(1/2) 0.996. The quarter-turned set's three sweeps, at three arm positions, now all index without the hand-passed quarter turn they needed and agree on one cell to 0.03 A. Six miniCBF sets that state no table and no +SLOW - including one whose Oscillation_axis line names an axis in a third dialect - are byte-identical in .hkl, .mtz, .cif and the image statistics, as are two NXmx sets, which is the shared orientation matcher moving nothing on that path either. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01T3yNBXk4wKdMZy1ak2NY7f --- common/DetectorOrientation.cpp | 12 +++ common/DetectorOrientation.h | 8 ++ docs/CHANGELOG.md | 2 + docs/DETECTOR_GEOMETRY.md | 24 +++++ reader/HDF5MetadataSource.cpp | 13 +-- reader/JFJochCBFReader.cpp | 75 ++++++++++++--- reader/MiniCBF.cpp | 134 +++++++++++++++++++++++++++ reader/MiniCBF.h | 16 ++++ tests/JFJochReaderTest.cpp | 162 +++++++++++++++++++++++++++++++++ 9 files changed, 419 insertions(+), 27 deletions(-) diff --git a/common/DetectorOrientation.cpp b/common/DetectorOrientation.cpp index 4e7095c4c..e11ec56e7 100644 --- a/common/DetectorOrientation.cpp +++ b/common/DetectorOrientation.cpp @@ -46,6 +46,18 @@ RotMatrix DetectorOrientation::Matrix() const { return {rz_x, mirror_y ? -rz_y : rz_y, {0, 0, 1}}; } +std::optional DetectorOrientation::Match(const Coord &fast, const Coord &slow) { + for (int64_t quarter_turns = 0; quarter_turns < 4; quarter_turns++) { + for (const bool mirror_y: {false, true}) { + const DetectorOrientation orientation(mirror_y, quarter_turns); + const RotMatrix m = orientation.Matrix(); + if (((m.Column(0) - fast).Length() < 1e-3f) && ((m.Column(1) - slow).Length() < 1e-3f)) + return orientation; + } + } + return {}; +} + bool DetectorOrientation::operator==(const DetectorOrientation &other) const { return (mirror_y == other.mirror_y) && (quarter_turns == other.quarter_turns); } diff --git a/common/DetectorOrientation.h b/common/DetectorOrientation.h index e9a48a2ce..d6a25d43a 100644 --- a/common/DetectorOrientation.h +++ b/common/DetectorOrientation.h @@ -3,6 +3,8 @@ #pragma once +#include + #include "Coord.h" // How the stored image is laid out in the detector plane: mirrored in Y, and/or turned by a multiple @@ -38,5 +40,11 @@ public: // Rz(quarter_turns * 90 deg) * diag(1,-1,1)^mirror_y. Entries are exactly 0 and +-1. [[nodiscard]] RotMatrix Matrix() const; + // The orientation whose fast and slow axes are the two directions given, or nothing when they are + // not one of the eight. Nothing means the image is turned in its own plane by something that is + // not a multiple of 90 degrees, which is a continuous rotation of the detector and belongs in the + // PONI angles - it cannot be told apart from the tilt by looking at the two directions alone. + [[nodiscard]] static std::optional Match(const Coord &fast, const Coord &slow); + bool operator==(const DetectorOrientation &other) const; }; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 862358d7a..4e4d39a7e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,8 @@ # Changelog ## 1.0.0 ### 1.0.0-rc.166 +* A miniCBF sweep takes the mounting from the imgCIF axis table its header carries: which laboratory direction the image's columns and rows run along, and which the spindle turns about. Where there is no table, a `+SLOW` on the `Oscillation_axis` line still says the spindle runs along the image's slow direction. All three were previously assumed, and two instruments are not what was assumed. +* Naming a miniCBF frame with no directory - a frame in the working directory - finds its sweep instead of reporting that no images were found. * A detector swung out on a 2theta arm is placed where it stands, from the depends_on transformation chain of an NXmx master or the `Detector_2theta` line of a miniCBF header; both were previously read and then ignored. * `jfjoch_viewer` opens PILATUS miniCBF sweeps - naming any frame opens the whole sweep - and can run a processing job on one. * A detector whose stored image is mirrored in Y or mounted at a multiple of 90 degrees can be described as such, in the detector configuration or with `--detector-mirror-y` / `--detector-quarter-turns`, rather than having to be expressed as a detector rotation. diff --git a/docs/DETECTOR_GEOMETRY.md b/docs/DETECTOR_GEOMETRY.md index 3cc216692..2b62eeafe 100644 --- a/docs/DETECTOR_GEOMETRY.md +++ b/docs/DETECTOR_GEOMETRY.md @@ -59,6 +59,30 @@ rot2 = asin(-slow.z) rot1 = atan2(-fast.z, normal.z) rot3 = atan with `rot2` in [-90°, 90°]. The angles are what is stored and what is written out, so a geometry given as angles comes back exactly as it was given. +## What a miniCBF header states about the mounting + +A PILATUS miniCBF gives the geometry twice. The `# ` lines every writer produces carry the distance, the +beam centre and the angles; some beamlines then append a CBF template block holding a full **imgCIF axis +table**, which states the laboratory direction of the image's fast and slow pixel directions, of the base +goniometer axis, and of a 2theta arm where there is one. Where that table is present it is read, in +preference to any assumption - it is the same information NXmx puts in `fast_pixel_direction` / +`slow_pixel_direction` and the goniometer `vector`, in the form this format states it. + +imgCIF's laboratory frame has Z from the sample towards the source and Y opposite gravity, so it differs +from the internal frame by a half turn about x - a rotation, not a mirror, so an axis carried through it +turns the same way by the same angle. + +Two things a header can state that an assumption gets wrong by 90 degrees, which no refinement recovers +and which the run's axis-sign rescue cannot reach either, a quarter turn not being a sign: + +* the image mounted a quarter turn round, so its columns run vertically; +* a spindle that turns about the **vertical** rather than the horizontal. + +Where a header carries no axis table, a `+SLOW` on its `# Oscillation_axis` line still says the spindle +runs along the image's slow direction rather than its fast one. The axis *name* on that line is not +usable - one header says `X.CW +SLOW` where its own table says the axis is Y - but the direction token is, +and on the header that states both they agree. + ## A detector swung out on a 2theta arm Chemical crystallography reaches high angle by swinging the detector out on a 2theta arm rather than by diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index ac68aaa60..20bef0df3 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -42,18 +42,7 @@ static std::optional ReadModuleOrientation(HDF5Object *file if ((f.size() != 3) || (s.size() != 3)) return {}; - const Coord fast = McStasToInternal(f); - const Coord slow = McStasToInternal(s); - - for (int64_t quarter_turns = 0; quarter_turns < 4; quarter_turns++) { - for (bool mirror_y: {false, true}) { - const DetectorOrientation orientation(mirror_y, quarter_turns); - const RotMatrix m = orientation.Matrix(); - if (((m.Column(0) - fast).Length() < 1e-3f) && ((m.Column(1) - slow).Length() < 1e-3f)) - return orientation; - } - } - return {}; + return DetectorOrientation::Match(McStasToInternal(f), McStasToInternal(s)); } // Where the detector stands, from the chain of transformations the file says it depends on. diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index a08639663..c60033e11 100644 --- a/reader/JFJochCBFReader.cpp +++ b/reader/JFJochCBFReader.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -57,7 +58,12 @@ std::optional