mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 78 Wiki Activity

rugnux: the per-reflection correction factor is named for what it is, not for what it once held

Browse Source 
The factor multiplied into each integrated intensity was called rlp, for reciprocal
Lorentz-polarization, and until this week that is all it held. It now also carries
the sensor efficiency at the angle the beam arrives, and on the stills path it holds
that efficiency and the polarization with no Lorentz term at all - correctly, since
the Lorentz factor of a still is one. Three different products under one name that
promises exactly one of them, in code where the neighbouring member is the total
correction.

Rename it prescaling_corr: multiplicative, applied before scaling, therefore not a
scale, and silent about its contents - which is the point, since the contents have
now grown twice. It is also what DIALS calls the same product. The stills refinement
member spelled "1 / rlp" becomes inv_corr, and the comments and usage text that
promised "the Lorentz-polarization factor and nothing else" now say what is actually
there.

The Lorentz term keeps its own name where it is computed, because that name is
correct. The two external spellings are untouched: the CBOR key and the reflection
dataset are a published format, and a reader that meets an unknown key would take
the factor as zero, which both the merge key and the ingest treat as a reflection to
drop - so every reflection would vanish and the run would still exit zero.

No output changes: the merged and unmerged files of two full runs are byte for byte
what the previous binary wrote, four stored files from before the efficiency
correction still re-scale identically, and the reflection datasets of the process
file are unchanged.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
This commit is contained in:
leonarski_f
2026-09-05 13:41:16 +02:00
co-authored by Claude Opus 5
parent d6ccdbdb61
commit db9cc9106f
36 changed files with 134 additions and 113 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/WriteReflections.cpp
+28 -20
View File
@@ -666,16 +666,16 @@ std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &
++j;
double sum_p = 0.0, sum_I = 0.0, sum_var = 0.0, sum_var_bkg = 0.0;
double p_rlp = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0,
double p_corr = 0.0, p_frame = 0.0, p_x = 0.0, p_y = 0.0, p_delta_phi = 0.0, p_zeta = 0.0,
p_bkg = 0.0;
for (size_t m = i; m < j; ++m) {
const Reflection &r = *parts[m].r;
const double p = r.partiality;
sum_p += p;
sum_I += static_cast<double>(r.I) * r.rlp;
sum_var += static_cast<double>(r.sigma) * r.sigma * r.rlp * r.rlp;
sum_var_bkg += static_cast<double>(r.var_bkg) * r.rlp * r.rlp;
p_rlp += p * r.rlp;
sum_I += static_cast<double>(r.I) * r.prescaling_corr;
sum_var += static_cast<double>(r.sigma) * r.sigma * r.prescaling_corr * r.prescaling_corr;
sum_var_bkg += static_cast<double>(r.var_bkg) * r.prescaling_corr * r.prescaling_corr;
p_corr += p * r.prescaling_corr;
p_frame += p * r.image_number;
p_x += p * DetectorX(r);
p_y += p * DetectorY(r);
@@ -688,13 +688,16 @@ std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &
if (sum_p < min_partiality || sum_p < min_captured_fraction)
continue;
// The Lorentz-polarization factor is applied by the writer, which multiplies I by rlp, so
// divide the event's own factor back out of the sums here. LP is the same geometry for every
// part of one event to a median 2e-4, so the file's I/LP is still the raw count sum.
full.rlp = static_cast<float>(p_rlp / sum_p);
full.I = static_cast<float>(sum_I / full.rlp);
full.sigma = static_cast<float>(std::sqrt(sum_var) / full.rlp);
full.var_bkg = static_cast<float>(sum_var_bkg / (static_cast<double>(full.rlp) * full.rlp));
// The prescaling correction is applied by the writer, which multiplies I by prescaling_corr,
// so divide the event's own factor back out of the sums here. Every term in it - Lorentz,
// polarization, sensor efficiency - is a function of the reflection's position and rocking
// geometry, which move by almost nothing between the parts of one event (median 2e-4 on the
// Lorentz term), so the file's I/LP is still the raw count sum.
full.prescaling_corr = static_cast<float>(p_corr / sum_p);
full.I = static_cast<float>(sum_I / full.prescaling_corr);
full.sigma = static_cast<float>(std::sqrt(sum_var) / full.prescaling_corr);
full.var_bkg = static_cast<float>(sum_var_bkg
/ (static_cast<double>(full.prescaling_corr) * full.prescaling_corr));
full.partiality = static_cast<float>(sum_p);
full.image_number = static_cast<float>(p_frame / sum_p);
full.observed_x = static_cast<float>(p_x / sum_p);
@@ -762,11 +765,16 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
// them there, so the index the reflection was actually measured at is recoverable - that is the
// crystal-frame information careless scales on, and what makes the file unmerged rather than a
// merge waiting to happen.
// I and SIGI are the integrated intensity with the Lorentz-polarization factor applied and
// nothing else, which is what IOBS means in every unmerged format (LP records the factor, so the
// raw counts are I/LP). LP is geometry, not a scale, and a program that reads this file has no
// way to recover it. The partiality is NOT divided out - that is a scale, FRACTIONCALC carries
// it, and every program this file is for wants to handle it its own way.
// I and SIGI are the integrated intensity with the prescaling correction applied and nothing
// else, which is what IOBS means in every unmerged format, and LP records the factor that was
// applied, so the raw counts are I/LP. That factor is geometry and instrument response, not a
// scale, and a program reading this file has no way to recompute it. The partiality is NOT
// divided out - that is a scale, FRACTIONCALC carries it, and every program this file is for
// wants to handle it its own way.
// Note that LP therefore carries slightly more than its name says: the prescaling correction now
// also holds the sensor's angle-dependent efficiency, for which DIALS writes a separate QE column
// and AIMLESS calls LP "the Lorentz/polarization correction (already applied)". A consumer that
// divides LP out to recover raw counts removes the efficiency term with it.
gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
// The batch headers are written for every image the observations span, not only for the images
// that produced one. AIMLESS starts a new run wherever the phi series jumps, so an image that
@@ -800,13 +808,13 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
mtz.data.push_back(static_cast<float>(hkl[2]));
mtz.data.push_back(static_cast<float>((partials ? 256 : 0) + isym));
mtz.data.push_back(static_cast<float>(batch));
mtz.data.push_back(r.I * r.rlp);
mtz.data.push_back(r.sigma * r.rlp);
mtz.data.push_back(r.I * r.prescaling_corr);
mtz.data.push_back(r.sigma * r.prescaling_corr);
mtz.data.push_back(r.partiality);
mtz.data.push_back(DetectorX(r));
mtz.data.push_back(DetectorY(r));
mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f);
mtz.data.push_back(r.rlp);
mtz.data.push_back(r.prescaling_corr);
mtz.data.push_back(0.0f); // FLAG: nothing here is a rejected observation
mtz.data.push_back(r.delta_phi_deg);
mtz.data.push_back(r.zeta);