rugnux: the unmerged MTZ carries the sensor efficiency in its own column, not inside LP
LP is a column other programs read and un-apply. It is documented as the Lorentz-polarization factor, and until the sensor efficiency was folded into the same product that is what it held. Measured on our own unmerged file, LP spanned a factor of 1.1343 across the detector where pure L/P spans nothing of the sort - the excess is the efficiency, 13% end to end at 13 keV and 34% at 18 keV. Both reference implementations keep it out. Recomputing pure L/P from a stored XDS file's own geometry over its 124k reflections leaves RLP flat to 0.1% from 8.6 to 33.7 degrees, where a folded-in efficiency would have shown a 7% trend - and XDS has the sensor numbers in hand. DIALS fills LP from lorentz and polarization alone and writes QE as a separate column, even a column of ones. Split them: LP is L/P again, QE is the efficiency, and the intensity is unchanged, so a reader that un-applies LP recovers what it expects and one that wants raw counts divides by LP and multiplies by QE. Only the unmerged file moves - every other column is bit-identical and the file grows by exactly one float per reflection plus one header record. The merged files are byte for byte what they were. The process file gains an optional qe dataset beside the existing one rather than changing what that one means, so a file written before the efficiency existed still loads, and reads back as a correction of exactly one - which is what it was. It is stored rather than recomputed on read because the writer has no geometry to recompute it from, and because recomputing would have written a radial trend into every stored file that never had one. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
@@ -10,6 +10,8 @@
|
||||
|
||||
#include "../common/DiffractionExperiment.h"
|
||||
#include "../image_analysis/WriteReflections.h"
|
||||
#include "../image_analysis/IntegrationOutcome.h"
|
||||
#include <cmath>
|
||||
#include "SyntheticMergedReflections.h"
|
||||
|
||||
namespace {
|
||||
@@ -70,3 +72,59 @@ TEST_CASE("Merged MTZ: the data dataset is id 1 and carries the wavelength", "[w
|
||||
|
||||
CHECK(mtz.nreflections == static_cast<int>(reflections.size()));
|
||||
}
|
||||
|
||||
TEST_CASE("Unmerged MTZ: LP is Lorentz-polarization and QE carries the sensor efficiency",
|
||||
"[write_reflections]") {
|
||||
// The whole point of the split: LP must mean what XDS and DIALS mean by it, and the raw count
|
||||
// sum must still be recoverable from the file alone, as I / LP * QE.
|
||||
auto experiment = TestExperiment();
|
||||
experiment.Goniometer(GoniometerAxis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}));
|
||||
|
||||
IntegrationOutcome outcome;
|
||||
const float raw[3] = {1000.0f, 250.0f, 40.0f};
|
||||
const float corr[3] = {1.75f, 2.50f, 0.90f}; // the whole prescaling product
|
||||
const float qe[3] = {0.9375f, 0.8125f, 1.0f}; // 1.0 = the sensor said nothing to correct
|
||||
for (int i = 0; i < 3; ++i) {
|
||||
Reflection r{};
|
||||
r.h = 4 + i; r.k = 2; r.l = 6;
|
||||
r.image_number = static_cast<float>(i);
|
||||
r.d = 5.0f + i;
|
||||
r.I = raw[i]; // the writer is what applies the factor
|
||||
r.sigma = std::sqrt(raw[i]);
|
||||
r.prescaling_corr = corr[i];
|
||||
r.qe_corr = qe[i];
|
||||
r.partiality = 1.0f;
|
||||
r.predicted_x = 100.0f + i; r.predicted_y = 200.0f + i;
|
||||
r.observed_x = NAN; r.observed_y = NAN;
|
||||
outcome.reflections.push_back(r);
|
||||
}
|
||||
|
||||
const auto path = (std::filesystem::temp_directory_path() / "rugnux_unmerged_qe.mtz").string();
|
||||
WriteUnmergedMtzReflections({outcome}, TETRAGONAL_CELL, experiment, false, path);
|
||||
const gemmi::Mtz mtz = gemmi::read_mtz_file(path);
|
||||
std::filesystem::remove(path);
|
||||
|
||||
const gemmi::Mtz::Column *c_I = mtz.column_with_label("I");
|
||||
const gemmi::Mtz::Column *c_lp = mtz.column_with_label("LP");
|
||||
const gemmi::Mtz::Column *c_qe = mtz.column_with_label("QE");
|
||||
REQUIRE(c_I != nullptr);
|
||||
REQUIRE(c_lp != nullptr);
|
||||
REQUIRE(c_qe != nullptr); // DIALS writes this column even when there is nothing in it
|
||||
REQUIRE(mtz.nreflections == 3);
|
||||
|
||||
for (int i = 0; i < 3; ++i) {
|
||||
const float I = mtz.data[i * mtz.columns.size() + c_I->idx];
|
||||
const float LP = mtz.data[i * mtz.columns.size() + c_lp->idx];
|
||||
const float QE = mtz.data[i * mtz.columns.size() + c_qe->idx];
|
||||
INFO("row " << i);
|
||||
// LP holds Lorentz x polarization alone: the sensor term has been divided out of it.
|
||||
CHECK(LP == Catch::Approx(corr[i] / qe[i]).epsilon(1e-5));
|
||||
// QE is a divisor normalised to 1 at normal incidence, so it never drops below 1.
|
||||
CHECK(QE == Catch::Approx(1.0f / qe[i]).epsilon(1e-5));
|
||||
CHECK(QE >= 1.0f);
|
||||
// ... and the two together put the raw counts back.
|
||||
CHECK(I / LP * QE == Catch::Approx(raw[i]).epsilon(1e-4));
|
||||
// The intensity itself did not move: it is still the fully corrected value.
|
||||
CHECK(I == Catch::Approx(raw[i] * corr[i]).epsilon(1e-5));
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user