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

rugnux: sum each rotation reflection's partials into one full on export

Browse Source 
--export-unmerged now writes one row per reflection with its partials
summed, which is how a rotation dataset is normally handed to a scaling
program. The previous behaviour is still available as
--export-unmerged-partials.

Partials made the file's usefulness depend on the reading program's
partial-assembly settings. On an orthorhombic case only 31% of partial
sets passed aimless's default PARTIALS test, so it merged 2% of the file:
multiplicity 4.7 against our own 12.4, CC1/2 0.788, and CC against our
merged intensities 0.806. Summed, aimless uses 89% of the file and rejects
nothing: multiplicity 12.5, CC1/2 0.931, completeness 96.1 to 98.0,
pointless Laue confidence 0.669 to 0.792, CC 0.848. The other two test
crystals are ties overall and win their outer shells, and the anomalous
signal comes out stronger from the fulls. iotbx.merging_statistics, which
never assembled partials at all, now reports real merging statistics
rather than statistics of parts.

Events are grouped by Combine()'s own rule so the file and the merge
assemble exactly the same reflections; the sum is plain with variances in
quadrature and no partiality division, which reproduces the merge's full
times FRACTIONCALC to a few percent.

Per-part metadata is combined on a partiality-weighted centroid rather
than an intensity-weighted one: 8-41% of events have parts summing to zero
or less, where intensity weights are undefined or sign-flipped, and 10-21%
of intensity-weighted centroids fall outside the event's own frame range.

Events below --min-partiality are not written. Declaring a reflection whose
rocking curve caught a thousandth of itself as a complete measurement gives
a scaling program noise with a small sigma and it weights it heavily; there
is no partial test to catch it once it is a full. This is the same constant
already applied to the same quantity in the merge, and it is what takes the
orthorhombic case from CC 0.219 to 0.848. --min-partiality 0 writes every
event.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
leonarski_f
2026-08-25 16:18:14 +02:00
co-authored by Claude Opus 5
parent 86b7e07131
commit bb5aaca68e
6 changed files with 186 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/WriteReflections.cpp
+137 -33
View File
@@ -538,9 +538,104 @@ void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
out.close();
}
namespace {
// The detector position an observation is written at: a reflection too weak to have a measured
// centroid still has a predicted one.
float DetectorX(const Reflection &r) { return std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x; }
float DetectorY(const Reflection &r) { return std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y; }
// Sum each rocking event into one full observation. A rotation reflection is integrated image by
// image, so it arrives here as a run of partials over consecutive frames; the run is cut where the
// 3D combine cuts it - same raw hkl, frames no further apart than MAX_FRAME_GAP - so the exported
// file and rugnux's own merge see exactly the same events.
// The parts are added, plainly, with their variances in quadrature, which is what every other
// rotation program writes as a full. Nothing is divided by the partiality: FRACTIONCALC carries the
// summed rocking-curve fraction, and an event the sweep cut short says so there rather than being
// scaled up or dropped here.
// The metadata a full carries is the partiality-weighted mean of its parts - the centroid of the
// rocking curve, which is what XD/YD/ZD mean in an XDS file and what a scale model downstream reads
// these columns for. An intensity weighting is the other candidate and it does not survive weak
// data: the parts of a weak reflection scatter about zero, so on the crystals this was measured on
// 8-41% of events have parts summing to zero or less and 10-21% of the intensity-weighted centroids
// fall outside the event's own frame range. Partialities are the rocking curve's own weights and are
// positive by construction.
// min_partiality is the combine's own floor on the assembled full (--min-partiality, default 0.02,
// 0 = off). An event that caught a thousandth of its rocking curve is not a measurement of that
// reflection - rugnux does not merge one either - and writing it as a full hands the reading program
// a whole observation whose intensity is noise and whose sigma is small, so it is weighted heavily.
// On a crystal whose rocking curves span twenty frames those events are 7% of the file and cost
// POINTLESS the point group. Everything above the floor is written with its honest FRACTIONCALC,
// truncated edge-of-sweep events included, for the reader to cut where it wants.
std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
double min_partiality) {
constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event
std::vector<const Reflection *> parts;
for (const auto &outcome : outcomes)
for (const auto &r : outcome.reflections)
parts.push_back(&r);
std::sort(parts.begin(), parts.end(), [](const Reflection *a, const Reflection *b) {
return std::tie(a->h, a->k, a->l, a->image_number)
< std::tie(b->h, b->k, b->l, b->image_number);
});
std::vector<Reflection> fulls;
for (size_t i = 0; i < parts.size(); ) {
size_t j = i + 1;
while (j < parts.size() && parts[j]->h == parts[i]->h && parts[j]->k == parts[i]->k
&& parts[j]->l == parts[i]->l
&& parts[j]->image_number - parts[j - 1]->image_number <= MAX_FRAME_GAP)
++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,
p_bkg = 0.0;
for (size_t m = i; m < j; ++m) {
const Reflection &r = *parts[m];
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;
p_frame += p * r.image_number;
p_x += p * DetectorX(r);
p_y += p * DetectorY(r);
p_delta_phi += p * r.delta_phi_deg;
p_zeta += p * r.zeta;
p_bkg += p * r.bkg;
}
Reflection full = *parts[i];
i = j;
if (sum_p < min_partiality)
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));
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);
full.observed_y = static_cast<float>(p_y / sum_p);
full.delta_phi_deg = static_cast<float>(p_delta_phi / sum_p);
full.zeta = static_cast<float>(p_zeta / sum_p);
full.bkg = static_cast<float>(p_bkg / sum_p);
fulls.push_back(full);
}
return fulls;
}
} // namespace
void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes,
const UnitCell &unitCell,
const DiffractionExperiment &experiment,
bool sum_partials,
const std::string &filename) {
gemmi::Mtz mtz;
mtz.spacegroup = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1));
@@ -572,13 +667,15 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
mtz.add_column("BGVAR", 'R', 0, -1, false);
const auto gon = experiment.GetGoniometer();
// A scanning goniometer means rugnux integrated a rocking curve image by image, so every
// observation here is one PART of a reflection; on a still each is a whole measurement. The
// difference has to be declared in two places, because that is where POINTLESS and AIMLESS look
// for it: the batch header's data type, and the M flag that M/ISYM packs as 256*M + ISYM. Told
// neither, they take each partial for a whole reflection and neither the symmetry determination
// nor the scaling survives it.
const bool partials = gon && gon->IsScanning();
// A scanning goniometer means rugnux integrated a rocking curve image by image, so each
// observation is one PART of a reflection unless the parts are summed here; on a still each is a
// whole measurement either way, and consecutive stills are different crystals, so there is
// nothing to sum there. Which of the two the file holds has to be declared in two places, because
// that is where POINTLESS and AIMLESS look for it: the batch header's data type, and the M flag
// that M/ISYM packs as 256*M + ISYM. Told nothing, they take each partial for a whole reflection
// and neither the symmetry determination nor the scaling survives it.
const bool scanning = gon && gon->IsScanning();
const bool partials = scanning && !sum_partials;
const float wedge_deg = gon ? gon->GetWedge_deg() : 0.0f;
const auto phi_start_deg = [&](float image_number) {
return gon ? gon->GetAngle_deg(image_number) : 0.0f;
@@ -595,32 +692,39 @@ void WriteUnmergedMtzReflections(const std::vector<IntegrationOutcome> &outcomes
// it, and every program this file is for wants to handle it its own way.
gemmi::UnmergedHklMover hkl_mover(mtz.spacegroup);
std::set<int> batch_numbers;
for (const auto &outcome : outcomes) {
for (const auto &r : outcome.reflections) {
std::array<int, 3> hkl{r.h, r.k, r.l};
const int isym = hkl_mover.move_to_asu(hkl);
const int batch = 1 + static_cast<int>(std::lround(r.image_number));
batch_numbers.insert(batch);
mtz.data.push_back(static_cast<float>(hkl[0]));
mtz.data.push_back(static_cast<float>(hkl[1]));
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.partiality);
// A reflection too weak to have a measured centroid still has a predicted position, and a
// detector coordinate is what the scale models downstream read this column for.
mtz.data.push_back(std::isfinite(r.observed_x) ? r.observed_x : r.predicted_x);
mtz.data.push_back(std::isfinite(r.observed_y) ? r.observed_y : r.predicted_y);
mtz.data.push_back(phi_start_deg(r.image_number) + wedge_deg / 2.0f);
mtz.data.push_back(r.rlp);
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);
mtz.data.push_back(r.bkg);
mtz.data.push_back(r.var_bkg);
}
const auto add_row = [&](const Reflection &r) {
std::array<int, 3> hkl{r.h, r.k, r.l};
const int isym = hkl_mover.move_to_asu(hkl);
// A summed full's image_number is its rocking-curve centroid, so this is the batch the
// reflection is centred on - which is what a batch means for a full everywhere else.
const int batch = 1 + static_cast<int>(std::lround(r.image_number));
batch_numbers.insert(batch);
mtz.data.push_back(static_cast<float>(hkl[0]));
mtz.data.push_back(static_cast<float>(hkl[1]));
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.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(0.0f); // FLAG: nothing here is a rejected observation
mtz.data.push_back(r.delta_phi_deg);
mtz.data.push_back(r.zeta);
mtz.data.push_back(r.bkg);
mtz.data.push_back(r.var_bkg);
};
if (scanning && sum_partials) {
for (const auto &r : SumRockingEvents(outcomes,
experiment.GetScalingSettings().GetMinPartiality()))
add_row(r);
} else {
for (const auto &outcome : outcomes)
for (const auto &r : outcome.reflections)
add_row(r);
}
mtz.nreflections = static_cast<int>(mtz.data.size() / mtz.columns.size());