RESULT-CHANGING: readers: a Bruker sweep that turns backwards is stated forwards, and the polarization comes from the header

Bruker frames state a negative increment (INCREME -0.2). Passed on as it was, the frame's oscillation
width was negative, and the per-frame mosaicity fit, which takes it as the frame's width, ran to its
2 degree cap on every Bruker sweep: rocking curves 25 frames long, partialities a tenth of what they
are, and on a sparse sweep no rocking event captured enough of itself to make a full, so the merge was
empty. The same rotation is now stated about the opposite axis by the opposite angles (imgcif::SweepAxis,
for the SFRM reader and the full-imgCIF CBF path). On the SMART 6000 2theta = 0 run of a twinned cubic
insulin crystal: mosaicity 2.000 -> 0.098 deg, R_meas 26.5% -> 5.9%, ISa 5.4 -> 23.8; the distance
post-refinement, which refused before, now commits 103.0 -> 106.5 mm and the cell comes out 77.5 A where
it was 75.1 (deposited 78.0) - the header distance of these frames is short, and the data now say so.
An APEX II sweep that merged nothing (empty merge, crash) now merges at ISa 19.

Polarization: a full imgCIF header's _diffrn_radiation_wavelength.polarizn_source_ratio is read as the
degree, in the plane its polarizn_source_norm states. An SFRM frame states the value APEX's conversion
writes for it, 1 - cos^2 of the monochromator angle (FILTER2: 0.0444 at 12.17 deg), so an .sfrm sweep and
its imgCIF conversion agree; with no monochromator (mirrors) that is 0. These are laboratory frames, which
before took the 0.99 synchrotron default.

--gain is accepted for a CCD imgCIF (CountsInADU), as for the other formats in detector units.

The CBF test now decodes through LoadImage: since rc181 ReadRawImage hands the byte-offset stream on
undecoded.

Tests: JFJochCBFReader_FullImgCIF_Polarization, JFJochSFRMReader_Geometry (backward sweep, FILTER2,
no-monochromator polarization), the existing Bruker and CBF reader cases.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

reader/ImgCIF.cpp

@@ -10,6 +10,7 @@
 
 #include "../common/JFJochException.h"
 #include "../common/DiffractionGeometry.h"
+#include "../common/JFJochMath.h"
 
 namespace imgcif {
 
@@ -272,9 +273,18 @@ std::optional<Description> Parse(const std::string &text) {
     d.exposure_s = NumberOr(Item(text, "_diffrn_scan_frame.integration_time").value_or(""), 0.0);
     d.detector = Item(text, "_diffrn_detector.type").value_or(Item(text, "_diffrn_detector.detector").value_or(""));
     d.ccd = Item(text, "_diffrn_detector.detector").value_or("").find("CCD") != std::string::npos;
+    if (const auto ratio = Item(text, "_diffrn_radiation_wavelength.polarizn_source_ratio"))
+        d.polarization_ratio = NumberOr(*ratio, 0.0);
+    d.polarization_norm_deg = NumberOr(Item(text, "_diffrn_radiation_wavelength.polarizn_source_norm").value_or(""), 0.0);
     return d;
 }
 
+std::optional<float> PolarizationFactor(const Description &d) {
+    if (!d.polarization_ratio)
+        return {};
+    return static_cast<float>(*d.polarization_ratio * std::cos(2.0 * d.polarization_norm_deg * PI / 180.0));
+}
+
 Geometry Evaluate(const Description &d) {
     Geometry g;
 
@@ -351,6 +361,13 @@ Coord GoniometerAxisOf(const std::array<double, 3> &v) {
     return -ToInternal(v);
 }
 
+GoniometerAxis SweepAxis(const std::string &name, double start_deg, double increment_deg,
+                         const std::array<double, 3> &v) {
+    const double sign = increment_deg < 0.0 ? -1.0 : 1.0;
+    return GoniometerAxis(name, static_cast<float>(sign * start_deg), static_cast<float>(sign * increment_deg),
+                          GoniometerAxisOf(v) * static_cast<float>(sign), {});
+}
+
 DetectorPlane PlaneOf(const Geometry &g, const std::string &path) {
     const Coord origin = ToInternal(g.origin);
     const Coord fast = ToInternal(g.fast);

reader/ImgCIF.h

@@ -12,6 +12,7 @@
 
 #include "../common/Coord.h"
 #include "../common/DetectorOrientation.h"
+#include "../common/GoniometerAxis.h"
 
 // The parts of a full imgCIF header (the form CBFlib writes, and Bruker's APEX writes when it converts
 // .sfrm frames) that say where the detector and the goniometer stand. Unlike a PILATUS miniCBF, which
@@ -62,8 +63,19 @@ struct Description {
     std::string detector;
     // _diffrn_detector.detector names a CCD: its pixel values are ADU, not counted photons.
     bool ccd = false;
+    // _diffrn_radiation_wavelength.polarizn_source_ratio, (Ip - In) / (Ip + In), and
+    // polarizn_source_norm, the angle (deg) from the laboratory Y axis to the normal of the plane Ip
+    // is polarized in. Nothing where the header states no ratio.
+    std::optional<double> polarization_ratio;
+    double polarization_norm_deg = 0;
 };
 
+// The polarization as DiffractionExperiment::PolarizationFactor takes it - the degree, positive for a
+// beam polarized in the horizontal (laboratory X-Z) plane - or nothing where the header states none.
+// A plane at an angle psi from the horizontal enters as its projection, the degree times cos(2 psi),
+// as the correction knows only those two planes.
+std::optional<float> PolarizationFactor(const Description &d);
+
 // The description, or nothing when the header is not a full imgCIF one - no axis table, or no
 // statement of where the axes stand. Throws on a table that is there but cannot be read.
 std::optional<Description> Parse(const std::string &text);
@@ -94,6 +106,13 @@ Coord ToInternal(const std::array<double, 3> &v);
 // master's vector carries, so a right-handed imgCIF rotation about v is the internal -v.
 Coord GoniometerAxisOf(const std::array<double, 3> &v);
 
+// The sweep as GoniometerAxis takes it, turning the scanned axis v from start_deg by increment_deg per
+// frame. A sweep that turns backwards (Bruker frames state a negative increment) is the same rotation
+// about -v by the opposite angles, and is stated that way: the analysis takes the oscillation width of a
+// frame to be positive, as every other format states it.
+GoniometerAxis SweepAxis(const std::string &name, double start_deg, double increment_deg,
+                         const std::array<double, 3> &v);
+
 // The detector plane in the terms the geometry is kept in: the point of normal incidence (PONI) in
 // pixels, the distance along the normal, and the rotation that takes the detector frame to the
 // laboratory, split into the discrete way the image is stored and the continuous rest. The table

reader/JFJochCBFReader.cpp

@@ -326,6 +326,10 @@ void JFJochCBFReader::ReadFiles(const std::string &path) {
     // correction takes to contain the lab x axis; the degree it wants is 2f - 1.
     if (header0_.polarization_fraction.has_value())
         dataset_->experiment.PolarizationFactor(2.0f * static_cast<float>(*header0_.polarization_fraction) - 1.0f);
+    // A full imgCIF header states the degree itself, and the plane.
+    if (header0_.imgcif.has_value())
+        if (const auto p = imgcif::PolarizationFactor(*header0_.imgcif))
+            dataset_->experiment.PolarizationFactor(*p);
     dataset_->experiment.FrameTime(
         std::chrono::duration_cast<std::chrono::nanoseconds>(
             std::chrono::duration<double>(header0_.period_s)),
@@ -370,11 +374,12 @@ void JFJochCBFReader::ReadFiles(const std::string &path) {
     // series can be missing frames, and those are gaps in the rotation rather than images to close up.
     const auto layout = sweep::Place(frames, "CBFReader");
     files_ = layout.files;
-    dataset_->experiment.Goniometer(GoniometerAxis(geometry0.has_value() ? geometry0->scan_axis : header0_.axis_name,
-                                                   static_cast<float>(layout.start_deg),
-                                                   static_cast<float>(layout.increment_deg),
-                                                   geometry0.has_value() ? imgcif::GoniometerAxisOf(geometry0->scan_direction)
-                                                                         : RotationAxis(header0_), {}));
+    if (geometry0.has_value())
+        dataset_->experiment.Goniometer(imgcif::SweepAxis(geometry0->scan_axis, layout.start_deg, layout.increment_deg,
+                                                        geometry0->scan_direction));
+    else
+        dataset_->experiment.Goniometer(GoniometerAxis(header0_.axis_name, static_cast<float>(layout.start_deg),
+                                                       static_cast<float>(layout.increment_deg), RotationAxis(header0_), {}));
 
     dataset_->error_value = -1;
     dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(files_.size()));

reader/JFJochSFRMReader.cpp

@@ -3,6 +3,8 @@
 
 #include "JFJochSFRMReader.h"
 
+#include <cmath>
+
 #include "../common/JFJochException.h"
 #include "../common/JFJochMath.h"
 #include "../common/Logger.h"
@@ -68,6 +70,12 @@ imgcif::Description Chain(const sfrm::Header &h, const std::string &path) {
                               + ", which is not a rotation this reader supports");
     d.axes[AXIS_NAMES[h.axis]].increment = h.increment_deg;
     d.wavelength_A = h.wavelength_A;
+    // The polarization APEX's conversion states for the same frame: 1 - cos^2 of the monochromator's
+    // scattering angle (FILTER2), in the horizontal plane. A tube with mirrors and no monochromator
+    // states 0, unpolarized - these are laboratory frames, and the synchrotron default does not apply.
+    const double cos_mono = std::cos(h.monochromator_2theta_deg * PI / 180.0);
+    d.polarization_ratio = 1.0 - cos_mono * cos_mono;
+    d.polarization_norm_deg = 0;
     d.overload = h.saturation;
     d.exposure_s = h.exposure_s;
     d.detector = h.detector;
@@ -119,6 +127,7 @@ void JFJochSFRMReader::ReadFiles(const std::string &path) {
     dataset_->experiment.DetectorDistance_mm(plane0.distance_mm);
     dataset_->experiment.PoniRot1_rad(plane0.rot1).PoniRot2_rad(plane0.rot2).PoniRot3_rad(plane0.rot3);
     dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast<float>(header0_.wavelength_A));
+    dataset_->experiment.PolarizationFactor(imgcif::PolarizationFactor(Chain(header0_, files_[0])));
     if (header0_.exposure_s > 0.0)
         dataset_->experiment.FrameTime(
             std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::duration<double>(header0_.exposure_s)),
@@ -145,10 +154,8 @@ void JFJochSFRMReader::ReadFiles(const std::string &path) {
 
     const auto layout = sweep::Place(frames, "SFRMReader");
     files_ = layout.files;
-    dataset_->experiment.Goniometer(GoniometerAxis(geometry0.scan_axis,
-                                                   static_cast<float>(layout.start_deg),
-                                                   static_cast<float>(layout.increment_deg),
-                                                   imgcif::GoniometerAxisOf(geometry0.scan_direction), {}));
+    dataset_->experiment.Goniometer(imgcif::SweepAxis(geometry0.scan_axis, layout.start_deg, layout.increment_deg,
+                                                    geometry0.scan_direction));
 
     dataset_->error_value = -1;
     dataset_->experiment.ImagesPerTrigger(static_cast<int64_t>(files_.size()));

reader/SFRM.cpp

@@ -162,6 +162,7 @@ Header ParseHeader(const char *data, size_t size) {
     h.saturation = Num(kv, "CCDPARM", 4);
     for (size_t i = 0; i < 6; i++)
         h.detpar[i] = Num(kv, "DETPAR", i);
+    h.monochromator_2theta_deg = Num(kv, "FILTER2", 0);
     h.raw = kv;
 
     if (h.header_blocks <= 0 || h.nrows <= 0 || h.ncols <= 0)

reader/SFRM.h

@@ -54,6 +54,7 @@ struct Header {
     double exposure_s = 0;              // CUMULAT
     double saturation = 0;              // CCDPARM[4]; 0 where not stated
     std::array<double, 6> detpar{};     // DETPAR: three offsets, then pitch, roll, yaw (deg)
+    double monochromator_2theta_deg = 0;   // FILTER2[0]: the monochromator's scattering angle; 0 for none
     std::string detector;               // DETTYPE[0]
     std::map<std::string, std::string> raw;   // every key; repeated keys joined with a space
 };

rugnux/rugnux_cli.cpp

@@ -1926,8 +1926,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
                            && JFJochXTCReader::CanRead(input_file);
             if (input_is_xtc && !all_inputs(JFJochXTCReader::CanRead))
                 refuse_mixed();
-            if (gain_per_photon && !reader_ptr && !input_is_xtc && !input_is_smv && !input_is_marccd && !input_is_sfrm) {
-                logger.Error("--gain is for input in detector units (SACLA, CXI, XTC, SMV, marCCD, SFRM or "
+            if (gain_per_photon && !reader_ptr && !input_is_xtc && !input_is_smv && !input_is_marccd && !input_is_sfrm
+                && !input_is_cbf) {
+                logger.Error("--gain is for input in detector units (SACLA, CXI, XTC, SMV, marCCD, SFRM, a CCD imgCIF or "
                              "--crystfel-geometry input); {} is read in the units it states", input_file);
                 exit(EXIT_FAILURE);
             }
@@ -2045,6 +2046,11 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     const bool counts_in_adu = input_is_smv || input_is_marccd || input_is_sfrm || input_is_xtc || input_is_panel_table
                                || input_is_cxi || input_is_crystfel || (input_is_cbf && cbf_reader.CountsInADU());
     const bool xfel_input = input_is_panel_table || input_is_cxi || input_is_crystfel;
+    // A CBF is in detector units only where it is a full imgCIF that names a CCD (CountsInADU).
+    if (gain_per_photon && input_is_cbf && !counts_in_adu) {
+        logger.Error("--gain is for input in detector units; {} is a CBF of counted photons", input_file);
+        exit(EXIT_FAILURE);
+    }
     if (counts_in_adu) {
         const bool float_frames = xfel_input && xfel_reader.StoredFloat();
         const float unit = float_frames ? 1.0f / 16.0f : 1.0f;

tests/BrukerReaderTest.cpp

@@ -43,6 +43,7 @@ struct ImgCIFFrame {
     double omega_inc = 0, kappa_inc = 0, phi_inc = 0.5;
     double two_theta = 0;
     double distance = 50;
+    std::string polarization;    // extra _diffrn_radiation_wavelength lines, if any
 };
 
 // A full imgCIF frame of a kappa head scanning phi, with a detector on a 2theta arm. The image's
@@ -54,7 +55,8 @@ void WriteImgCIF(const std::string &name, const ImgCIFFrame &f, int64_t nx, int6
     std::ostringstream h;
     h << "###CBF: VERSION 1.7.11\n# CBF file written by CBFlib v0.9.6\n\ndata_image\n\n"
       << "_diffrn_scan_frame.integration_time 5\n"
-      << "_diffrn_radiation_wavelength.id SOURCE\n_diffrn_radiation_wavelength.wavelength 0.71073\n\n"
+      << "_diffrn_radiation_wavelength.id SOURCE\n_diffrn_radiation_wavelength.wavelength 0.71073\n"
+      << f.polarization << "\n"
       << "loop_\n_axis.id\n_axis.depends_on\n_axis.equipment\n_axis.type\n"
       << "_axis.vector[1]\n_axis.vector[2]\n_axis.vector[3]\n_axis.offset[1]\n_axis.offset[2]\n_axis.offset[3]\n"
       << " OMEGA . goniometer rotation 1 0 0 0 0 0\n"
@@ -148,17 +150,18 @@ TEST_CASE("JFJochCBFReader_FullImgCIF_KappaChain", "[HDF5][Full][portable]") {
         CHECK(x.GetPixelSize_mm() == Catch::Approx(0.12));
         CHECK(x.GetWavelength_A() == Catch::Approx(0.71073));
         CHECK(x.GetSaturationLimit() == 600000);
+        CHECK_FALSE(x.GetPolarizationFactor().has_value());   // the header states none
         // fast along laboratory +Y (internal -y), slow along +X
         const auto geom = x.GetDiffractionGeometry();
         CHECK((geom.GetFastAxis() - Coord(0, -1, 0)).Length() < 1e-5f);
         CHECK((geom.GetSlowAxis() - Coord(1, 0, 0)).Length() < 1e-5f);
         CHECK((geom.GetNormalAxis() - Coord(0, 0, 1)).Length() < 1e-5f);
 
-        JFJochReaderRawImage raw;
-        REQUIRE(reader.ReadRawImage(1, raw));
-        const auto *v = reinterpret_cast<const int32_t *>(raw.image.GetCompressed());
-        CHECK(v[0] == 5);
-        CHECK(v[17] == 1000);
+        // ReadRawImage hands the byte-offset stream on undecoded; LoadImage decodes it
+        const auto image = reader.LoadImage(1);
+        REQUIRE(image);
+        CHECK(image->Image()[0] == 5);
+        CHECK(image->Image()[17] == 1000);
         CHECK(reader.GetDataset()->pixel_mask->GetMask()[3] != 0);
         CHECK(reader.GetDataset()->pixel_mask->GetMask()[0] == 0);
         reader.Close();
@@ -191,6 +194,32 @@ TEST_CASE("JFJochCBFReader_FullImgCIF_KappaChain", "[HDF5][Full][portable]") {
 
 // What the full imgCIF path must never do: read a series it cannot describe as one rotation about
 // one axis as if it could.
+TEST_CASE("JFJochCBFReader_FullImgCIF_Polarization", "[HDF5][Full][portable]") {
+    const int64_t nx = 24, ny = 16;
+    const std::vector<int32_t> pixels(nx * ny, 5);
+    auto read = [&](const std::string &prefix, const std::string &lines) {
+        for (int frame = 1; frame <= 2; frame++) {
+            ImgCIFFrame f;
+            f.phi = 10 + 0.5 * (frame - 1);
+            f.polarization = lines;
+            WriteImgCIF(FrameName(prefix, frame, ".cbf"), f, nx, ny, pixels);
+        }
+        JFJochCBFReader reader;
+        reader.ReadFiles(FrameName(prefix, 1, ".cbf"));
+        return reader.GetDataset()->experiment.GetPolarizationFactor();
+    };
+    // the degree as stated, in the horizontal plane when the plane's normal is the laboratory Y axis
+    const auto horizontal = read("imgcif_pol_h", "_diffrn_radiation_wavelength.polarizn_source_ratio 0.0444\n"
+                                                 "_diffrn_radiation_wavelength.polarizn_source_norm 0\n");
+    REQUIRE(horizontal.has_value());
+    CHECK(*horizontal == Catch::Approx(0.0444));
+    // a normal turned by 90 degrees puts the plane vertical, which the correction takes as the opposite sign
+    const auto vertical = read("imgcif_pol_v", "_diffrn_radiation_wavelength.polarizn_source_ratio 0.5\n"
+                                               "_diffrn_radiation_wavelength.polarizn_source_norm 90\n");
+    REQUIRE(vertical.has_value());
+    CHECK(*vertical == Catch::Approx(-0.5));
+}
+
 TEST_CASE("JFJochCBFReader_FullImgCIF_refuses_what_it_cannot_describe", "[HDF5][Full][portable]") {
     const int64_t nx = 8, ny = 8;
     const std::vector<int32_t> pixels(nx * ny, 1);
@@ -252,6 +281,7 @@ struct SFRMFrame {
     double start = 0, increment = 0.5;
     std::array<double, 4> angles{0, 0, 0, 0};    // 2theta omega phi chi
     double distance_cm = 5.0;
+    double monochromator_2theta = 0;   // FILTER2
     std::vector<uint8_t> pixels;
     std::vector<uint16_t> over16;
     std::vector<int32_t> over32;
@@ -295,6 +325,7 @@ void WriteSFRM(const std::string &name, const SFRMFrame &f, int64_t n) {
     h += Line("NEXP", "2 0 32 0 1");
     h += Line("CCDPARM", "14 10 170 0 900000");
     h += Line("DETPAR", "0 0 0 0 0 0");
+    h += Line("FILTER2", std::to_string(f.monochromator_2theta) + " 0");
     h.resize(15 * 512, ' ');
     std::string body(f.pixels.begin(), f.pixels.end());
     Pad(body, 16);
@@ -402,6 +433,9 @@ TEST_CASE("JFJochSFRMReader_Geometry", "[HDF5][Full][portable]") {
         CHECK(x.GetBeamY_pxl() == Catch::Approx((n - 1) / 2.0).margin(1e-3));
         CHECK(x.GetSaturationLimit() == 900000);
         CHECK(x.GetWavelength_A() == Catch::Approx(1.54184));
+        // laboratory frames with no monochromator: unpolarized, not the synchrotron default
+        REQUIRE(x.GetPolarizationFactor().has_value());
+        CHECK(*x.GetPolarizationFactor() == Catch::Approx(0.0).margin(1e-6));
         CHECK(reader.GetDataset()->pixel_mask->GetMask()[0] != 0);
         CHECK(reader.GetDataset()->pixel_mask->GetMask()[1] == 0);
 
@@ -411,6 +445,29 @@ TEST_CASE("JFJochSFRMReader_Geometry", "[HDF5][Full][portable]") {
         reader.Close();
     }
 
+    SECTION("a sweep that turns backwards is stated forwards about the opposite axis") {
+        for (int frame = 1; frame <= 3; frame++) {
+            f.start = 10.0 - 0.5 * (frame - 1);
+            f.increment = -0.5;
+            f.angles = {0, f.start, 0, 54.74};
+            f.monochromator_2theta = 12.17;
+            WriteSFRM(FrameName("sfrm_back", frame, ".sfrm"), f, n);
+        }
+        JFJochSFRMReader reader;
+        REQUIRE_NOTHROW(reader.ReadFiles("sfrm_back_0001.sfrm"));
+        CHECK(reader.GetNumberOfImages() == 3);
+        const auto x = reader.GetDataset()->experiment;
+        REQUIRE(x.GetGoniometer().has_value());
+        CHECK((x.GetGoniometer()->GetAxis() - Coord(-1, 0, 0)).Length() < 1e-5f);
+        CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5));
+        CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(-10.0));
+        // a graphite monochromator: the 1 - cos^2(2theta_m) APEX's conversion states, horizontal
+        REQUIRE(x.GetPolarizationFactor().has_value());
+        const double c = std::cos(12.17 * PI / 180.0);
+        CHECK(*x.GetPolarizationFactor() == Catch::Approx(1.0 - c * c));
+        reader.Close();
+    }
+
     SECTION("phi scan on a fixed-chi head turns about the inclined phi axis") {
         for (int frame = 1; frame <= 2; frame++) {
             f.axis = 3;