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 e2c0875ec..4e4d39a7e 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -1,6 +1,9 @@ # 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. * 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..2b62eeafe 100644 --- a/docs/DETECTOR_GEOMETRY.md +++ b/docs/DETECTOR_GEOMETRY.md @@ -59,6 +59,45 @@ 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 +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..20bef0df3 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,18 +42,64 @@ 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]); + return DetectorOrientation::Match(McStasToInternal(f), 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; - } +// 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"); } - return {}; + + // 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; } @@ -721,6 +773,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..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