diff --git a/common/DiffractionGeometry.cpp b/common/DiffractionGeometry.cpp index 6a0908dad..ed8946a3a 100644 --- a/common/DiffractionGeometry.cpp +++ b/common/DiffractionGeometry.cpp @@ -131,9 +131,9 @@ float DiffractionGeometry::CalcAzIntSolidAngleCorr(float x, float y) const { // // A vertically polarized source is expressed by a NEGATIVE coeff: flipping the sign of the cos(2phi) // term is exactly a 90 degree rotation of the plane, since cos(2(phi+90 deg)) = -cos(2phi). That is -// what the jfjoch_settings API means by "1.0 is horizontal polarization; -1.0 is vertical". It is -// not autodetected: nothing in the files rugnux reads declares a polarization plane, so the sign has -// to be given with --polarization if a beamline ever needs the other one. +// what the jfjoch_settings API means by "1.0 is horizontal polarization; -1.0 is vertical". Of the +// files rugnux reads only a d*TREK header states the plane (SOURCE_POLARZ); for every other one the +// sign has to be given with --polarization if a beamline ever needs the other one. float DiffractionGeometry::CalcAzIntPolarizationCorr(float x, float y, float coeff) const { auto cos_2theta = cosf(TwoTheta_rad(x, y)); float cos_2theta_2 = cos_2theta * cos_2theta; diff --git a/docs/CPU_DATA_ANALYSIS_IMAGE.md b/docs/CPU_DATA_ANALYSIS_IMAGE.md index 668fb8284..8ceb64824 100644 --- a/docs/CPU_DATA_ANALYSIS_IMAGE.md +++ b/docs/CPU_DATA_ANALYSIS_IMAGE.md @@ -462,12 +462,15 @@ detector tilt and a swung-out $2\theta$ arm without further work. The **polarization plane** is taken to contain the lab $x$ axis — a horizontally polarized source. A vertically polarized one is expressed by a **negative** $P$: flipping the sign of the $\cos 2\phi$ -term is exactly a $90^\circ$ rotation of the plane. The plane is **not** autodetected and is not read -from the file — no format Rugnux reads declares one — so a vertical beamline must say so with -`--polarization`. $P$ is the polarization *degree* (XDS's `FRACTION_OF_POLARIZATION` is $(1+P)/2$). -The default $P = 0.99$, an undulator value, is applied to every Rugnux run; it is not fitted, because -a single dataset does not determine it. Bending-magnet and wiggler beamlines are typically lower -(0.8–0.95), and a known value for one should be passed with `--polarization`. +term is exactly a $90^\circ$ rotation of the plane. $P$ is the polarization *degree* (XDS's +`FRACTION_OF_POLARIZATION` is $(1+P)/2$). It is not fitted, because a single dataset does not +determine it, so it comes from the file where the file states it: a miniCBF `# Polarization` line and +a d*TREK `SOURCE_POLARZ` both give the fraction $f$, read as $P = 2f - 1$ (so a laboratory tube's 0.5 +is unpolarized), and `SOURCE_POLARZ` also gives the plane's normal, which sets the sign. Where the file +says nothing (ADSC SMV, marCCD, NXmx) the default $P = 0.99$, an undulator value, applies, with the +plane horizontal; a vertical beamline must say so with `--polarization`. Bending-magnet and wiggler +beamlines are typically lower (0.8–0.95), and a known value for one should be passed with +`--polarization`. ### 2.3 Background estimate for profiles diff --git a/reader/JFJochCBFReader.cpp b/reader/JFJochCBFReader.cpp index 8cec01a4c..f3d52e8b6 100644 --- a/reader/JFJochCBFReader.cpp +++ b/reader/JFJochCBFReader.cpp @@ -286,6 +286,10 @@ void JFJochCBFReader::ReadFiles(const std::string &path) { } dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast(header0_.wavelength_A)); + // The header states the fraction f polarized in the horizontal plane, which is the plane the + // 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); dataset_->experiment.FrameTime( std::chrono::duration_cast( std::chrono::duration(header0_.period_s)), diff --git a/reader/JFJochSMVReader.cpp b/reader/JFJochSMVReader.cpp index a2eeb31cd..56369e8aa 100644 --- a/reader/JFJochSMVReader.cpp +++ b/reader/JFJochSMVReader.cpp @@ -125,6 +125,14 @@ void JFJochSMVReader::ReadFiles(const std::string &path) { } dataset_->experiment.IncidentEnergy_keV(WVL_1A_IN_KEV / static_cast(header0_.wavelength_A)); + // The polarization degree is 2f - 1 of the fraction f the header states. Its sign is the plane: + // positive where the plane contains the lab x axis (normal along y), negative where it contains y. + // A laboratory tube says 0.5, unpolarized, where the plane does not matter. + if (header0_.polarization_fraction.has_value()) { + const Coord normal = ImgCIFToInternal(*header0_.polarization_normal); + const float degree = 2.0f * static_cast(*header0_.polarization_fraction) - 1.0f; + dataset_->experiment.PolarizationFactor(std::fabs(normal.x) > std::fabs(normal.y) ? -degree : degree); + } // Only when the header states a sane one. The exposure field is not always filled in: one // deposited sweep carries -2093438692 there, which is not a time at all, and passing it on // refuses the whole dataset over a number that affects no geometry and no result. A day is a diff --git a/reader/MiniCBF.cpp b/reader/MiniCBF.cpp index 3793b0cc1..955976f6d 100644 --- a/reader/MiniCBF.cpp +++ b/reader/MiniCBF.cpp @@ -267,6 +267,8 @@ Header ParseHeader(const char *data, size_t size) { h.exposure_s = Num(t, R"(#\s*Exposure_time\s+([\d.eE+-]+))"); h.period_s = Num(t, R"(#\s*Exposure_period\s+([\d.eE+-]+))"); h.count_cutoff = Int(t, R"(#\s*Count_cutoff\s+(\d+))"); + if (Match(t, R"(#\s*Polarization\s+([\d.eE+-]+))").has_value()) + h.polarization_fraction = Num(t, R"(#\s*Polarization\s+([\d.eE+-]+))"); h.axis_name = Match(t, R"(#\s*Oscillation_axis\s+(\S+))").value_or("omega"); // "+SLOW" / "+FAST" on that same line: which of the image's two directions the spindle runs // along. Some writers state that instead of an axis name, and it is the only thing a header with diff --git a/reader/MiniCBF.h b/reader/MiniCBF.h index 6b33be351..3b82eaafb 100644 --- a/reader/MiniCBF.h +++ b/reader/MiniCBF.h @@ -45,6 +45,9 @@ struct Header { double exposure_s = 0; double period_s = 0; int64_t count_cutoff = 0; // saturation + // "# Polarization 0.990": the fraction of the beam polarized in the horizontal plane (0.5 is + // unpolarized). Absent where the line is missing or says "None". + std::optional polarization_fraction; std::string axis_name = "omega"; bool byte_offset = false; // the only conversion supported diff --git a/reader/SMV.cpp b/reader/SMV.cpp index 16c30011c..db520b463 100644 --- a/reader/SMV.cpp +++ b/reader/SMV.cpp @@ -197,6 +197,13 @@ void ParseDtrek(const std::map &kv, const std::string h.wavelength_A = !scan_wavelength.empty() ? scan_wavelength[0] : (source_wavelength.size() >= 2 ? source_wavelength[1] : 0.0); + // Read as dxtbx's FormatSMVRigaku reads it: the fraction first, then the plane normal. + const auto polarz = Nums(kv, "SOURCE_POLARZ"); + if (polarz.size() >= 4) { + h.polarization_fraction = polarz[0]; + h.polarization_normal = std::array{polarz[1], polarz[2], polarz[3]}; + } + h.overflow_ratio = static_cast(Num(kv, "RAXIS_COMPRESSION_RATIO")); const auto id = kv.find(det + "DETECTOR_IDENTIFICATION"); if (id != kv.end()) diff --git a/reader/SMV.h b/reader/SMV.h index 708a5d17a..83b4b9136 100644 --- a/reader/SMV.h +++ b/reader/SMV.h @@ -66,6 +66,10 @@ struct Header { double angle_deg; }; std::vector detector_circles; // header order: the outermost circle first + // SOURCE_POLARZ: the fraction of the beam polarized in a plane, and that plane's normal (the + // order XDS states FRACTION_OF_POLARIZATION and POLARIZATION_PLANE_NORMAL in). 0.5 is unpolarized. + std::optional polarization_fraction; + std::optional> polarization_normal; // RAXIS_COMPRESSION_RATIO: a stored value above 32767 is (value - 32768) * ratio. 0 = not used. int64_t overflow_ratio = 0; std::map raw; // every key, for anything not modelled above diff --git a/rugnux/RugnuxDefaults.cpp b/rugnux/RugnuxDefaults.cpp index 9b5d8c4e0..f827c8675 100644 --- a/rugnux/RugnuxDefaults.cpp +++ b/rugnux/RugnuxDefaults.cpp @@ -15,7 +15,10 @@ void ApplyRugnuxExperimentDefaults(DiffractionExperiment &experiment) { // azimuthal, intensity-proportional modulation (~1.3x at 2 A for a 1 A beam) that symmetry mates do // NOT share, because they sit at the same 2-theta but a different azimuth - which is exactly the // "unequal intensities forced together" signature the space-group search vetoes as pseudo-symmetry. - experiment.PolarizationFactor(0.99); + // Only where the file stated none: a miniCBF or d*TREK header says what its source was, and a + // laboratory tube there says unpolarized. + if (!experiment.GetPolarizationFactor().has_value()) + experiment.PolarizationFactor(0.99); // The background ring's radial curvature correction: auto, applied per image where the ice score // sees a smooth ring. The library default stays off (the broker keeps it off); offline, over ice diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 129487a92..1999ac00e 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -256,7 +256,7 @@ void print_usage() { std::cout << " --rot3 PONI rotation 3, about the beam (rad)" << std::endl; std::cout << " --detector-mirror-y Stored image is mirrored in Y vs the detector" << std::endl; std::cout << " --detector-quarter-turns <0-3> Stored image is turned by this many 90 deg about the beam" << std::endl; - std::cout << " --polarization Degree of polarization p (XDS FRACTION_OF_POLARIZATION = (1+p)/2), signed: positive for a horizontally polarized source, negative for a vertically polarized one (the sign rotates the polarization plane by 90 degrees). Default 0.99, an undulator value; bending-magnet and wiggler sources are typically lower (0.8-0.95), and no file format rugnux reads states it reliably, so pass a known beamline value here" << std::endl; + std::cout << " --polarization Degree of polarization p (XDS FRACTION_OF_POLARIZATION = (1+p)/2), signed: positive for a horizontally polarized source, negative for a vertically polarized one (the sign rotates the polarization plane by 90 degrees). Default: what a miniCBF \"# Polarization\" line or a d*TREK SOURCE_POLARZ states (a fraction f, read as p = 2f - 1, so 0.5 from a laboratory tube is unpolarized), and otherwise 0.99, an undulator value; bending-magnet and wiggler sources are typically lower (0.8-0.95), and neither ADSC SMV, marCCD nor NXmx files state it, so pass a known beamline value here" << std::endl; } enum { @@ -1816,14 +1816,17 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) { const auto gonio = x.GetGoniometer(); const auto [direct_x, direct_y] = x.GetDiffractionGeometry().GetDirectBeam_pxl(); logger.Info("Header: {} images, detector {} {} x {} px of {:.4f} mm, wavelength {:.5f} A, distance " - "{:.3f} mm, beam {:.2f} {:.2f} px (direct beam {:.2f} {:.2f}), tilt {:.4f} {:.4f} deg, {}", + "{:.3f} mm, beam {:.2f} {:.2f} px (direct beam {:.2f} {:.2f}), tilt {:.4f} {:.4f} deg, {}, " + "polarization {}", reader.GetNumberOfImages(), x.GetDetectorDescription(), x.GetXPixelsNum(), x.GetYPixelsNum(), x.GetPixelSize_mm(), x.GetWavelength_A(), x.GetDetectorDistance_mm(), x.GetBeamX_pxl(), x.GetBeamY_pxl(), direct_x, direct_y, x.GetPoniRot1_rad() * 180.0 / PI, x.GetPoniRot2_rad() * 180.0 / PI, gonio ? fmt::format("rotation from {:.3f} deg by {:.4f} deg per image", gonio->GetStart_deg(), gonio->GetIncrement_deg()) - : std::string("no rotation axis")); + : std::string("no rotation axis"), + x.GetPolarizationFactor() ? fmt::format("{:.3f}", *x.GetPolarizationFactor()) + : std::string("not stated")); console->OnNote(fmt::format("{} images, {} {} x {} px", reader.GetNumberOfImages(), x.GetDetectorDescription(), x.GetXPixelsNum(), x.GetYPixelsNum())); console->OnNote(fmt::format("{:.4f} A, {:.1f} mm, {}", x.GetWavelength_A(), x.GetDetectorDistance_mm(), diff --git a/tests/JFJochReaderTest.cpp b/tests/JFJochReaderTest.cpp index 32d044c49..0362cd493 100644 --- a/tests/JFJochReaderTest.cpp +++ b/tests/JFJochReaderTest.cpp @@ -4049,6 +4049,33 @@ namespace { } } +TEST_CASE("JFJochCBFReader_Polarization", "[HDF5][Full][portable]") { + const int64_t nx = 24, ny = 16; + // The line states the fraction polarized in the horizontal plane; the correction wants the degree. + WriteMiniCBFSweep("cbfpol", "# Polarization 0.990\n", nx, ny); + JFJochCBFReader reader; + REQUIRE_NOTHROW(reader.ReadFiles("cbfpol_0001.cbf")); + REQUIRE(reader.GetDataset()->experiment.GetPolarizationFactor().has_value()); + CHECK(reader.GetDataset()->experiment.GetPolarizationFactor().value() == Catch::Approx(0.98)); + reader.Close(); + RemoveMiniCBFSweep("cbfpol"); + + WriteMiniCBFSweep("cbfunpol", "# Polarization 0.500000\n", nx, ny); // a laboratory tube + JFJochCBFReader lab; + REQUIRE_NOTHROW(lab.ReadFiles("cbfunpol_0001.cbf")); + REQUIRE(lab.GetDataset()->experiment.GetPolarizationFactor().has_value()); + CHECK(lab.GetDataset()->experiment.GetPolarizationFactor().value() == Catch::Approx(0.0).margin(1e-6)); + lab.Close(); + RemoveMiniCBFSweep("cbfunpol"); + + WriteMiniCBFSweep("cbfnopol", "# Polarization None\n", nx, ny); + JFJochCBFReader none; + REQUIRE_NOTHROW(none.ReadFiles("cbfnopol_0001.cbf")); + CHECK_FALSE(none.GetDataset()->experiment.GetPolarizationFactor().has_value()); + none.Close(); + RemoveMiniCBFSweep("cbfnopol"); +} + TEST_CASE("JFJochCBFReader_AxisTableStatesTheMounting", "[HDF5][Full][portable]") { const int64_t nx = 24, ny = 16; @@ -4553,6 +4580,7 @@ TEST_CASE("JFJochSMVReader_Geometry", "[HDF5][Full][portable]") { REQUIRE(x.GetGoniometer().has_value()); CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(5.0)); CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.2)); + CHECK_FALSE(x.GetPolarizationFactor().has_value()); // an ADSC header states none reader.Close(); RemoveSMVSweep("smv_", ".img"); } @@ -4572,7 +4600,8 @@ TEST_CASE("JFJochSMVReader_Geometry", "[HDF5][Full][portable]") { << "CCD_GONIO_VECTORS=0 0 1 1 0 0 0 1 0 1 0 0 0 1 0 0 0 -1;\n" << "ROTATION=-13.500 -13.000 0.500 30.000 0.000 0.000 0.000 100.000 0.000 0.000;\n" << "ROTATION_AXIS_NAME=Omega;\nROTATION_VECTOR=1.000 0.000 0.000;\n" - << "SOURCE_WAVELENGTH=1.000000 1.541780;\nSATURATED_VALUE=4194176;\nRAXIS_COMPRESSION_RATIO=128;\n}\n"; + << "SOURCE_WAVELENGTH=1.000000 1.541780;\nSATURATED_VALUE=4194176;\nRAXIS_COMPRESSION_RATIO=128;\n" + << "SOURCE_POLARZ=0.9000 1.0000 0.0000 0.0000;\n}\n"; std::string head = h.str(); head.resize(1024, ' '); std::vector pixels(static_cast(nx) * ny, 100); @@ -4613,6 +4642,10 @@ TEST_CASE("JFJochSMVReader_Geometry", "[HDF5][Full][portable]") { CHECK(x.GetGoniometer()->GetStart_deg() == Catch::Approx(-13.5)); CHECK(x.GetGoniometer()->GetIncrement_deg() == Catch::Approx(0.5)); CHECK((x.GetGoniometer()->GetAxis() - Coord(1, 0, 0)).Length() < 1e-5f); + // A fraction of 0.9 is a degree of 0.8; the plane's normal along x puts the plane on y, which + // the correction spells as a negative degree. + REQUIRE(x.GetPolarizationFactor().has_value()); + CHECK(x.GetPolarizationFactor().value() == Catch::Approx(-0.8)); reader.Close(); remove("dtrek_0001.img"); }