diff --git a/reader/ImgCIF.cpp b/reader/ImgCIF.cpp index ab8eb6699..9d1800b36 100644 --- a/reader/ImgCIF.cpp +++ b/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 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 PolarizationFactor(const Description &d) { + if (!d.polarization_ratio) + return {}; + return static_cast(*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 &v) { return -ToInternal(v); } +GoniometerAxis SweepAxis(const std::string &name, double start_deg, double increment_deg, + const std::array &v) { + const double sign = increment_deg < 0.0 ? -1.0 : 1.0; + return GoniometerAxis(name, static_cast(sign * start_deg), static_cast(sign * increment_deg), + GoniometerAxisOf(v) * static_cast(sign), {}); +} + DetectorPlane PlaneOf(const Geometry &g, const std::string &path) { const Coord origin = ToInternal(g.origin); const Coord fast = ToInternal(g.fast); diff --git a/reader/ImgCIF.h b/reader/ImgCIF.h index d621b8ff1..fedce33c8 100644 --- a/reader/ImgCIF.h +++ b/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 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 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 Parse(const std::string &text); @@ -94,6 +106,13 @@ Coord ToInternal(const std::array &v); // master's vector carries, so a right-handed imgCIF rotation about v is the internal -v. Coord GoniometerAxisOf(const std::array &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 &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 diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index 59ec50222..09c3d5be1 100644 --- a/reader/JFJochCBFReader.cpp +++ b/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(*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::duration(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(layout.start_deg), - static_cast(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(layout.start_deg), + static_cast(layout.increment_deg), RotationAxis(header0_), {})); dataset_->error_value = -1; dataset_->experiment.ImagesPerTrigger(static_cast(files_.size())); diff --git a/reader/JFJochSFRMReader.cpp b/reader/JFJochSFRMReader.cpp index 665336d9e..5e92928c1 100644 --- a/reader/JFJochSFRMReader.cpp +++ b/reader/JFJochSFRMReader.cpp @@ -3,6 +3,8 @@ #include "JFJochSFRMReader.h" +#include + #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(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::duration(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(layout.start_deg), - static_cast(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(files_.size())); diff --git a/reader/SFRM.cpp b/reader/SFRM.cpp index d3ba9bf64..f8da7bf4d 100644 --- a/reader/SFRM.cpp +++ b/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) diff --git a/reader/SFRM.h b/reader/SFRM.h index 9ac9a256b..47795c8db 100644 --- a/reader/SFRM.h +++ b/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 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 raw; // every key; repeated keys joined with a space }; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index e5eeafd01..0e2156f43 100644 --- a/rugnux/rugnux_cli.cpp +++ b/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; diff --git a/tests/BrukerReaderTest.cpp b/tests/BrukerReaderTest.cpp index 23b7b82e1..d1e25ea1d 100644 --- a/tests/BrukerReaderTest.cpp +++ b/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(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 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 pixels(nx * ny, 1); @@ -252,6 +281,7 @@ struct SFRMFrame { double start = 0, increment = 0.5; std::array angles{0, 0, 0, 0}; // 2theta omega phi chi double distance_cm = 5.0; + double monochromator_2theta = 0; // FILTER2 std::vector pixels; std::vector over16; std::vector 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;