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

rugnux: write standard CCP4 anomalous MTZ (IMEAN + I(+)/I(-))

Browse Source 
In anomalous mode the merge keeps the two Friedel mates as separate rows,
which WriteMtzReflections emitted verbatim - two IMEAN rows per Bijvoet pair
that downstream tools had to re-collapse. Pair the mates into one row per
reflection with the standard CCP4 anomalous layout: IMEAN + I(+)/I(-) and the
matching F/F(+)/F(-) amplitudes. A single HKLKeyGenerator(merge_friedel=false)
yields both the shared ASU group key and which mate a row is (.plus); IMEAN/F
are the inverse-variance Friedel mean, centrics/unpaired keep one hand with the
other missing. aimless/mtz2sca/ANODE now read the file directly. The
non-anomalous MTZ path is unchanged.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-16 10:37:56 +02:00
co-authored by Claude Opus 4.8
parent ce15900471
commit ce12d69751
2 changed files with 102 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/WriteReflections.cpp
+101 -20
View File
@@ -3,9 +3,12 @@
#include "WriteReflections.h"
#include "scale_merge/Merge.h"
#include "scale_merge/HKLKey.h"
#include "scale_merge/TwinningAnalysis.h"
#include <cmath>
#include <map>
#include <tuple>
#include <fstream>
#include <iomanip>
#include <sstream>
@@ -220,36 +223,114 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
ds.crystal_name = experiment.GetSampleName();
ds.wavelength = experiment.GetWavelength_A();
int dataset_id = ds.id;
const int dataset_id = ds.id;
// In anomalous mode the merge keeps the two Friedel mates as separate rows (I+ under the ASU
// representative hkl, I- under -hkl). Emitting those verbatim gives a file with two rows per
// reflection that downstream tools have to re-collapse. Instead pair the mates into one row per
// reflection with the standard CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
// the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
if (experiment.GetScalingSettings().GetMergeFriedel()) {
mtz.add_column("H", 'H', dataset_id, -1, false);
mtz.add_column("K", 'H', dataset_id, -1, false);
mtz.add_column("L", 'H', dataset_id, -1, false);
mtz.add_column("IMEAN", 'J', dataset_id, -1, false);
mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false);
mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude
mtz.add_column("SIGF", 'Q', dataset_id, -1, false);
mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false);
mtz.nreflections = static_cast<int>(reflections.size());
mtz.data.reserve(reflections.size() * 8);
for (const auto& r : reflections) {
mtz.data.push_back(static_cast<float>(r.h));
mtz.data.push_back(static_cast<float>(r.k));
mtz.data.push_back(static_cast<float>(r.l));
mtz.data.push_back(r.I);
mtz.data.push_back(r.sigma);
mtz.data.push_back(r.F);
mtz.data.push_back(r.sigmaF);
mtz.data.push_back(r.rfree_flag ? 1.0f : 0.0f);
}
mtz.write_to_file(filename);
return;
}
// Anomalous: group the two mates by their (shared) Friedel-merged ASU representative. A single
// generator gives both the group key (its hkl, identical for +hkl and -hkl) and which mate this
// row is (.plus).
const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1));
struct AnomRow {
int h = 0, k = 0, l = 0;
float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
int rfree = 0;
};
std::map<std::tuple<int, int, int>, AnomRow> rows;
for (const auto& r : reflections) {
const HKLKey key = key_gen(r);
AnomRow& row = rows[{key.h, key.k, key.l}];
row.h = key.h; row.k = key.k; row.l = key.l;
row.rfree = r.rfree_flag ? 1 : 0;
if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; }
else { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; }
}
// Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured).
const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) {
const bool ok_a = std::isfinite(a) && sa > 0.0f;
const bool ok_b = std::isfinite(b) && sb > 0.0f;
if (ok_a && ok_b) {
const double wa = 1.0 / (static_cast<double>(sa) * sa);
const double wb = 1.0 / (static_cast<double>(sb) * sb);
val = static_cast<float>((wa * a + wb * b) / (wa + wb));
sig = static_cast<float>(1.0 / std::sqrt(wa + wb));
} else if (ok_a) { val = a; sig = sa; }
else if (ok_b) { val = b; sig = sb; }
else { val = NAN; sig = NAN; }
};
// Add columns
mtz.add_column("H", 'H', dataset_id, -1, false);
mtz.add_column("K", 'H', dataset_id, -1, false);
mtz.add_column("L", 'H', dataset_id, -1, false);
mtz.add_column("IMEAN", 'J', dataset_id, -1, false);
mtz.add_column("SIGIMEAN", 'Q', dataset_id, -1, false);
mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude
mtz.add_column("I(+)", 'K', dataset_id, -1, false);
mtz.add_column("SIGI(+)", 'M', dataset_id, -1, false);
mtz.add_column("I(-)", 'K', dataset_id, -1, false);
mtz.add_column("SIGI(-)", 'M', dataset_id, -1, false);
mtz.add_column("F", 'F', dataset_id, -1, false); // French-Wilson amplitude (mean)
mtz.add_column("SIGF", 'Q', dataset_id, -1, false);
mtz.add_column("F(+)", 'G', dataset_id, -1, false);
mtz.add_column("SIGF(+)", 'L', dataset_id, -1, false);
mtz.add_column("F(-)", 'G', dataset_id, -1, false);
mtz.add_column("SIGF(-)", 'L', dataset_id, -1, false);
mtz.add_column("FreeR_flag", 'I', dataset_id, -1, false);
// Number of rows
mtz.nreflections = reflections.size();
// Allocate data table
mtz.data.reserve(reflections.size() * 8);
for (const auto& r : reflections) {
mtz.data.push_back(static_cast<float>(r.h));
mtz.data.push_back(static_cast<float>(r.k));
mtz.data.push_back(static_cast<float>(r.l));
mtz.data.push_back(r.I);
mtz.data.push_back(r.sigma);
mtz.data.push_back(r.F);
mtz.data.push_back(r.sigmaF);
mtz.data.push_back(r.rfree_flag ? 1 : 0);
mtz.nreflections = static_cast<int>(rows.size());
mtz.data.reserve(rows.size() * 16);
for (const auto& [hkl, row] : rows) {
float i_mean, sig_i_mean, f_mean, sig_f_mean;
combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean);
combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean);
mtz.data.push_back(static_cast<float>(row.h));
mtz.data.push_back(static_cast<float>(row.k));
mtz.data.push_back(static_cast<float>(row.l));
mtz.data.push_back(i_mean);
mtz.data.push_back(sig_i_mean);
mtz.data.push_back(row.Ip);
mtz.data.push_back(row.sIp);
mtz.data.push_back(row.Im);
mtz.data.push_back(row.sIm);
mtz.data.push_back(f_mean);
mtz.data.push_back(sig_f_mean);
mtz.data.push_back(row.Fp);
mtz.data.push_back(row.sFp);
mtz.data.push_back(row.Fm);
mtz.data.push_back(row.sFm);
mtz.data.push_back(static_cast<float>(row.rfree));
}
// Write MTZ
mtz.write_to_file(filename);
}