leonarski_fandClaude Opus 5
762c093b10
The export declared its rows FULL - LDTYPE=2 in the batch header, M=0 in M/ISYM - while filtering them on --min-partiality alone (0.02), where the 3D combine also applies --min-captured-fraction (0.7 on rotation) to the same summed event. AIMLESS reads the FULL declaration, reports "all runs have only fulls" and ignores FRACTIONCALC, so an event that caught a twentieth of its rocking curve entered the merge whole, with a small sigma, and was weighted heavily. There was no cut for the reader to make: the column that would let it make one is the one the reading program has been told to ignore. Dropping those events moves AIMLESS's Rmerge at 1.8 A from 1.353 to 0.694 and CC(1/2) from 0.985 to 0.993. FRACTIONCALC itself is unchanged, values above 1 included. Each part's partiality is the erf pair the predictor computed on that frame, from that frame's own refined lattice and mosaicity, so the parts of one event do not tile the rocking curve exactly - the offset steps by the wedge to within 12% of it, and the sums land in a peak at 1.000 whose 95th percentile is 1.09. It is an honest estimate of a captured fraction, and it is not the number rugnux scales on, so clamping it would only hide the spread. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
910 lines
49 KiB
C++
910 lines
49 KiB
C++
// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "WriteReflections.h"
|
|
#include "scale_merge/Merge.h"
|
|
#include "scale_merge/HKLKey.h"
|
|
#include "scale_merge/TwinningAnalysis.h"
|
|
#include "bragg_integration/SystematicAbsence.h"
|
|
#include "../common/ParallelFor.h"
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
#include <cstdio>
|
|
#include <map>
|
|
#include <tuple>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <sstream>
|
|
#include <stdexcept>
|
|
#include <ctime>
|
|
#include <chrono>
|
|
|
|
#include <gemmi/mtz.hpp>
|
|
|
|
#include "../common/GitInfo.h"
|
|
|
|
namespace {
|
|
|
|
/// Current date in ISO-8601 (YYYY-MM-DD) for the _audit block.
|
|
std::string CurrentDateISO() {
|
|
auto now = std::chrono::system_clock::now();
|
|
auto t = std::chrono::system_clock::to_time_t(now);
|
|
std::tm tm{};
|
|
#ifdef _WIN32
|
|
gmtime_s(&tm, &t);
|
|
#else
|
|
gmtime_r(&t, &tm);
|
|
#endif
|
|
char buf[32];
|
|
std::strftime(buf, sizeof(buf), "%Y-%m-%d", &tm);
|
|
return buf;
|
|
}
|
|
|
|
/// Format a double with given decimal places; returns "?" for non-finite.
|
|
/// snprintf rather than an ostringstream: the reflection loop below calls this twelve times per
|
|
/// reflection, and building a stream (and its locale) per call dominated the time spent writing a
|
|
/// merged file. Same digits - both go through the C locale's %.*f.
|
|
std::string Fmt(double val, int decimals = 4) {
|
|
if (!std::isfinite(val))
|
|
return "?";
|
|
char buf[512];
|
|
const int n = std::snprintf(buf, sizeof(buf), "%.*f", decimals, val);
|
|
return std::string(buf, n);
|
|
}
|
|
|
|
/// Quote a CIF string value; returns "?" for empty.
|
|
std::string CifStr(const std::string& s) {
|
|
if (s.empty())
|
|
return "?";
|
|
// If it contains spaces or special chars, single-quote it
|
|
if (s.find(' ') != std::string::npos ||
|
|
s.find('\'') != std::string::npos ||
|
|
s.find('#') != std::string::npos)
|
|
return "'" + s + "'";
|
|
return s;
|
|
}
|
|
|
|
// One output row per reflection in the standard CCP4 anomalous layout: the merged mean (IMEAN / F)
|
|
// plus the two Bijvoet mates (I(+)/I(-), F(+)/F(-)). The merge keeps the two mates as separate rows
|
|
// (I+ under the ASU representative hkl, I- under -hkl); this collapses them into one row so the MTZ
|
|
// and SHELX writers share one row list. has_anom is set false when NO reflection carries an anomalous
|
|
// split (e.g. the stills path) - callers then omit the +/- columns.
|
|
struct MergedOutRow {
|
|
int h = 0, k = 0, l = 0;
|
|
float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
|
|
float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
|
|
int rfree = 0;
|
|
};
|
|
|
|
std::vector<MergedOutRow> BuildMergedRows(const std::vector<MergedReflection> &reflections,
|
|
const DiffractionExperiment &experiment,
|
|
bool &has_anom) {
|
|
std::vector<MergedOutRow> out_rows;
|
|
has_anom = true;
|
|
|
|
if (experiment.GetScalingSettings().GetMergeFriedel()) {
|
|
// Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim
|
|
// from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak
|
|
// anomalous signal is preserved without reprocessing. A reflection with only one mate, or a
|
|
// centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous
|
|
// split (e.g. the stills path, which does not compute one) the anomalous columns are omitted.
|
|
has_anom = std::any_of(reflections.begin(), reflections.end(),
|
|
[](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
|
|
out_rows.reserve(reflections.size());
|
|
for (const auto& r : reflections)
|
|
out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus,
|
|
r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus,
|
|
r.sigmaF_minus, r.rfree_flag ? 1 : 0});
|
|
} else {
|
|
// Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form
|
|
// IMEAN / F as their inverse-variance Friedel mean. 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.GetSpaceGroupOrP1());
|
|
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; }
|
|
};
|
|
|
|
out_rows.reserve(rows.size());
|
|
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);
|
|
out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im,
|
|
row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree});
|
|
}
|
|
}
|
|
return out_rows;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
|
|
void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const MergeStatistics &statistics,
|
|
const ErrorModelReport &error_model,
|
|
const TwinningAnalysisResult &twinning,
|
|
const std::string &filename,
|
|
size_t nthreads) {
|
|
|
|
std::ofstream out(filename);
|
|
if (!out)
|
|
throw std::runtime_error("WriteMmcifReflections: cannot open " + filename);
|
|
|
|
out << std::fixed;
|
|
|
|
// ---------- data block ----------
|
|
out << "data_sample" << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _audit ----------
|
|
out << "_audit.revision_id 1\n";
|
|
out << "_audit.creation_date " << CurrentDateISO() << "\n";
|
|
out << "_audit.update_record 'Initial release'\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _software ----------
|
|
out << "_software.name 'Rugnux'\n";
|
|
out << "_software.version " << CifStr(jfjoch_version()) << "\n";
|
|
out << "_software.classification 'data reduction'\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _cell ----------
|
|
out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n";
|
|
out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n";
|
|
out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n";
|
|
out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n";
|
|
out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n";
|
|
out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n";
|
|
|
|
const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1();
|
|
|
|
// ---------- _symmetry ----------
|
|
out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n";
|
|
out << "_symmetry.Int_Tables_number " << sg.number << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- _diffrn_source / _diffrn_detector ----------
|
|
if (!experiment.GetSourceName().empty())
|
|
out << "_diffrn_source.pdbx_synchrotron_site " << CifStr(experiment.GetSourceName()) << "\n";
|
|
|
|
if (!experiment.GetInstrumentName().empty())
|
|
out << "_diffrn_source.pdbx_synchrotron_beamline " << CifStr(experiment.GetInstrumentName()) << "\n";
|
|
|
|
out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n";
|
|
out << "_diffrn_detector.detector " << CifStr(experiment.GetDetectorDescription()) << "\n";
|
|
// Detector geometry actually used for integration (refined, when geometry refinement ran - the rotation
|
|
// two-pass or the stills global refinement update it on experiment_ before the written pass). jfjoch_
|
|
// local-data-name items: feedback of the distance / beam centre the data was reduced with.
|
|
out << "_diffrn_detector.jfjoch_distance_mm " << Fmt(experiment.GetDetectorDistance_mm(), 4) << "\n";
|
|
out << "_diffrn_detector.jfjoch_beam_center_x_pxl " << Fmt(experiment.GetBeamX_pxl(), 2) << "\n";
|
|
out << "_diffrn_detector.jfjoch_beam_center_y_pxl " << Fmt(experiment.GetBeamY_pxl(), 2) << "\n";
|
|
out << "#\n";
|
|
|
|
// ---------- merging statistics (_reflns overall + _reflns_shell loop) ----------
|
|
// cc_half and r_meas are stored as fractions (0-1), which is the mmCIF convention. ISa (the
|
|
// Diederichs asymptotic I/sigma, 1/b of the a*sigma^2 + (b*I)^2 error model) and the twinning
|
|
// indicators below have no standard mmCIF item. They are written under the "jfjoch" reserved
|
|
// prefix (_reflns.jfjoch_*), the IUCr-sanctioned local-data-name extension for private items -
|
|
// NOT the "pdbx_" prefix, which is owned by the wwPDB PDBx/mmCIF dictionary and must not label
|
|
// items that dictionary does not define. (The other pdbx_ items here are genuine PDBx items.)
|
|
const auto mult = [](const MergeStatisticsShell &s) {
|
|
return s.unique_reflections > 0 ? static_cast<double>(s.total_observations) / s.unique_reflections : 0.0; };
|
|
const auto compl_pct = [](const MergeStatisticsShell &s) {
|
|
return s.possible_unique_reflections > 0
|
|
? 100.0 * static_cast<double>(s.unique_reflections) / s.possible_unique_reflections : 0.0; };
|
|
if (!statistics.shells.empty()) {
|
|
const auto &ov = statistics.overall;
|
|
// Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a
|
|
// non-anomalous merge keeps its previous stats block / shell-loop columns unchanged.
|
|
const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous);
|
|
out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n";
|
|
out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n";
|
|
out << "_reflns.number_obs " << ov.unique_reflections << "\n";
|
|
out << "_reflns.pdbx_number_measured_all " << ov.total_observations << "\n";
|
|
out << "_reflns.pdbx_redundancy " << Fmt(mult(ov), 2) << "\n";
|
|
out << "_reflns.percent_possible_obs " << Fmt(compl_pct(ov), 1) << "\n";
|
|
out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n";
|
|
out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n";
|
|
out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n";
|
|
if (has_anom)
|
|
out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3)
|
|
<< " # SigAno = <|dano|>/<sigma(dano)>\n";
|
|
// ISa in XDS's sense: the whole-range 1/sqrt(a*b) of the error model below, so this item can
|
|
// be read straight against a CORRECT.LP. The strong-reflection asymptote - a tier XDS does not
|
|
// have, and always the more optimistic of the two - is written separately rather than here.
|
|
out << "_reflns.jfjoch_diffrn_ISa " << CifStr(error_model.isa)
|
|
<< " # 1/sqrt(a*b), the XDS convention\n";
|
|
if (!error_model.isa_asymptotic.empty())
|
|
out << "_reflns.jfjoch_diffrn_ISa_asymptotic " << CifStr(error_model.isa_asymptotic)
|
|
<< " # strong-reflection asymptote (Diederichs); rotation path only\n";
|
|
if (!error_model.a.empty())
|
|
out << "_reflns.jfjoch_error_model_a " << CifStr(error_model.a)
|
|
<< " # sigma^2 = a*(sigma0^2 + b*I^2), XDS convention\n";
|
|
if (!error_model.b.empty())
|
|
out << "_reflns.jfjoch_error_model_b " << CifStr(error_model.b) << "\n";
|
|
// Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's
|
|
// "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded.
|
|
if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0)
|
|
out << "_reflns.B_iso_Wilson_estimate " << Fmt(statistics.wilson_b, 2) << "\n";
|
|
// Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above).
|
|
if (twinning.l_test_pairs > 0) {
|
|
out << "_reflns.jfjoch_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3)
|
|
<< " # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n";
|
|
out << "_reflns.jfjoch_L_test_mean_L_squared " << Fmt(twinning.mean_l_squared, 3)
|
|
<< " # <L^2> (untwinned 0.333, perfect twin 0.200)\n";
|
|
}
|
|
if (twinning.moment_reflections > 0)
|
|
out << "_reflns.jfjoch_second_moment_I " << Fmt(twinning.second_moment, 3)
|
|
<< " # <I^2>/<I>^2 (untwinned 2.00, perfect twin 1.50)\n";
|
|
// Radiation-damage monitor (rotation): the relative Debye-Waller B change from the first to the last
|
|
// frame (A^2). A large magnitude flags a dose-dependent resolution-scale change = radiation damage;
|
|
// positive is the typical direction (high-resolution intensity fades with dose). No standard mmCIF item.
|
|
if (std::isfinite(statistics.radiation_damage_delta_b))
|
|
out << "_reflns.jfjoch_radiation_damage_relative_B " << Fmt(statistics.radiation_damage_delta_b, 2)
|
|
<< " # relative-B first->last over the run (A^2); + = high-res fades with dose\n";
|
|
// Diffraction anisotropy. The eigen-decomposition of the anisotropy tensor has standard PDBx
|
|
// items; the eigenvalues there must be non-negative, and only the deviatoric part of the tensor
|
|
// is determined at all (its isotropic part is degenerate with the overall scale), so they are
|
|
// written relative to the weakest direction - eigenvalue_3 is 0 by construction and
|
|
// eigenvalue_1 is the anisotropic deltaB. The eigenvectors are in the PDB orthogonalisation
|
|
// convention, which is the one gemmi (and hence rugnux) uses throughout.
|
|
//
|
|
// The directional diffraction LIMITS are deliberately NOT written as
|
|
// _reflns.pdbx_aniso_diffraction_limit_*: the dictionary defines those as the semi-axes of an
|
|
// ellipsoid fitted to a diffraction cut-off surface, which is a different construction from the
|
|
// one below and one rugnux does not perform - it cuts nothing on a directional criterion. They
|
|
// go under the jfjoch local prefix with their own definition instead.
|
|
const auto &an = statistics.anisotropy;
|
|
if (an.n_reflections > 0 && an.n_cells > 0 && std::isfinite(an.delta_b)) {
|
|
out << "_reflns.pdbx_orthogonalization_convention pdb\n";
|
|
for (int i = 0; i < 3; ++i) {
|
|
out << "_reflns.pdbx_aniso_B_tensor_eigenvalue_" << (i + 1) << " "
|
|
<< Fmt(an.eigenvalue[i] - an.eigenvalue[2], 2) << "\n";
|
|
for (int j = 0; j < 3; ++j)
|
|
out << "_reflns.pdbx_aniso_B_tensor_eigenvector_" << (i + 1) << "_ortho[" << (j + 1)
|
|
<< "] " << Fmt(an.eigenvector[i][j], 4) << "\n";
|
|
}
|
|
out << "_reflns.jfjoch_aniso_delta_B " << Fmt(an.delta_b, 2)
|
|
<< " # range of the principal components (A^2), fitted on intensities\n";
|
|
if (std::isfinite(an.delta_b_linear))
|
|
out << "_reflns.jfjoch_aniso_delta_B_linear " << Fmt(an.delta_b_linear, 2)
|
|
<< " # deltaB implied by the s^2 slope alone: what a Debye-Waller B accounts for\n";
|
|
out << "_reflns.jfjoch_aniso_shape " << AnisotropyShapeCode(an.shape)
|
|
<< " # resolution signature of the directional deficit\n";
|
|
if (std::isfinite(an.floor))
|
|
out << "_reflns.jfjoch_aniso_floor " << Fmt(an.floor, 3)
|
|
<< " # deltaB this data set's own systematic error could manufacture (A^2)\n";
|
|
if (std::isfinite(an.significance))
|
|
out << "_reflns.jfjoch_aniso_significance " << Fmt(an.significance, 2)
|
|
<< " # deltaB(linear) / floor\n";
|
|
out << "_reflns.jfjoch_aniso_verdict " << AnisotropyVerdictCode(an.verdict)
|
|
<< " # DETECTED / NOT_DETECTED / CANNOT_DETERMINE, at deltaB(linear)/floor > 3.5\n";
|
|
for (int i = 0; i < 3; ++i)
|
|
if (std::isfinite(an.d_min_axis[i]))
|
|
out << "_reflns.jfjoch_aniso_d_min_" << (i + 1) << " "
|
|
<< Fmt(an.d_min_axis[i], 2)
|
|
<< " # <I/sigma(I)> = 2 in a 20 deg cone about eigenvector " << (i + 1)
|
|
<< (an.d_min_censored[i] ? "; at the edge of the measured data, so a bound" : "")
|
|
<< "\n";
|
|
}
|
|
out << "#\n";
|
|
|
|
// Per-batch relative-B curve (the radiation-damage monitor, rotation): one relative Debye-Waller B
|
|
// per rotation-range batch, measured before any correction. rotation_start_deg = id * batch_deg.
|
|
if (!statistics.radiation_damage_b_batch.empty()) {
|
|
out << "loop_\n";
|
|
out << "_jfjoch_radiation_damage_batch.id\n";
|
|
out << "_jfjoch_radiation_damage_batch.rotation_start_deg\n";
|
|
out << "_jfjoch_radiation_damage_batch.relative_B\n";
|
|
for (size_t i = 0; i < statistics.radiation_damage_b_batch.size(); ++i)
|
|
out << " " << (i + 1) << " "
|
|
<< Fmt(static_cast<double>(i) * statistics.radiation_damage_batch_deg, 1) << " "
|
|
<< Fmt(statistics.radiation_damage_b_batch[i], 2) << "\n";
|
|
out << "#\n";
|
|
}
|
|
|
|
out << "loop_\n";
|
|
out << "_reflns_shell.d_res_high\n";
|
|
out << "_reflns_shell.d_res_low\n";
|
|
out << "_reflns_shell.number_measured_obs\n";
|
|
out << "_reflns_shell.number_unique_obs\n";
|
|
out << "_reflns_shell.pdbx_redundancy\n";
|
|
out << "_reflns_shell.percent_possible_obs\n";
|
|
out << "_reflns_shell.meanI_over_sigI_obs\n";
|
|
out << "_reflns_shell.pdbx_Rrim_I_all\n";
|
|
out << "_reflns_shell.pdbx_CC_half\n";
|
|
if (has_anom)
|
|
out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n";
|
|
for (const auto &s : statistics.shells) {
|
|
if (s.unique_reflections == 0)
|
|
continue;
|
|
out << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " "
|
|
<< s.total_observations << " " << s.unique_reflections << " "
|
|
<< Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " "
|
|
<< Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4);
|
|
if (has_anom)
|
|
out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3);
|
|
out << "\n";
|
|
}
|
|
out << "#\n";
|
|
}
|
|
|
|
// ---------- _refln loop ----------
|
|
out << "loop_\n";
|
|
out << "_refln.index_h\n";
|
|
out << "_refln.index_k\n";
|
|
out << "_refln.index_l\n";
|
|
out << "_refln.intensity_meas\n";
|
|
out << "_refln.intensity_sigma\n";
|
|
out << "_refln.pdbx_I_plus\n";
|
|
out << "_refln.pdbx_I_plus_sigma\n";
|
|
out << "_refln.pdbx_I_minus\n";
|
|
out << "_refln.pdbx_I_minus_sigma\n";
|
|
out << "_refln.F_meas_au\n";
|
|
out << "_refln.F_meas_sigma_au\n";
|
|
out << "_refln.pdbx_F_plus\n";
|
|
out << "_refln.pdbx_F_plus_sigma\n";
|
|
out << "_refln.pdbx_F_minus\n";
|
|
out << "_refln.pdbx_F_minus_sigma\n";
|
|
out << "_refln.status_free\n";
|
|
out << "_refln.status\n";
|
|
|
|
// One row per unique reflection, twelve formatted floats each - tens of megabytes on a crowded
|
|
// crystal, and the largest single-threaded stretch left in a run. Nothing about a row depends on
|
|
// any other, so each worker formats its own block into its own string and the blocks go to the
|
|
// file in order. The columns are written exactly as the stream wrote them: the same "%.4f" (or
|
|
// "?" where the value is not finite), right-aligned in the same width.
|
|
{
|
|
const size_t nrow = reflections.size();
|
|
const size_t nw = std::max<size_t>(nthreads, 1);
|
|
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
|
|
std::vector<std::string> block(nch);
|
|
ParallelChunks(nch, nw, [&](int tlo, int thi) {
|
|
for (int t = tlo; t < thi; ++t) {
|
|
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
|
|
std::string &s = block[t];
|
|
s.reserve((hi - lo) * 208);
|
|
const auto column = [&s](const std::string &v, size_t width) {
|
|
if (v.size() < width) s.append(width - v.size(), ' ');
|
|
s.append(v);
|
|
s.push_back(' ');
|
|
};
|
|
for (size_t i = lo; i < hi; ++i) {
|
|
const auto &r = reflections[i];
|
|
column(std::to_string(r.h), 5);
|
|
column(std::to_string(r.k), 5);
|
|
column(std::to_string(r.l), 5);
|
|
column(Fmt(r.I, 4), 14);
|
|
column(Fmt(r.sigma, 4), 14);
|
|
column(Fmt(r.I_plus, 4), 14);
|
|
column(Fmt(r.sigma_plus, 4), 14);
|
|
column(Fmt(r.I_minus, 4), 14);
|
|
column(Fmt(r.sigma_minus, 4), 14);
|
|
column(Fmt(r.F, 4), 14);
|
|
column(Fmt(r.sigmaF, 4), 14);
|
|
column(Fmt(r.F_plus, 4), 14);
|
|
column(Fmt(r.sigmaF_plus, 4), 14);
|
|
column(Fmt(r.F_minus, 4), 14);
|
|
column(Fmt(r.sigmaF_minus, 4), 14);
|
|
s.push_back(r.rfree_flag ? '1' : '0');
|
|
s.append(" o\n"); // 'o' = observed
|
|
}
|
|
}
|
|
});
|
|
for (const std::string &s : block)
|
|
out.write(s.data(), static_cast<std::streamsize>(s.size()));
|
|
}
|
|
|
|
out << "#\n";
|
|
out << "# End of reflections\n";
|
|
out.close();
|
|
}
|
|
|
|
void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const std::string &filename) {
|
|
gemmi::Mtz mtz;
|
|
|
|
// Optional but recommended metadata
|
|
mtz.spacegroup = &experiment.GetSpaceGroupOrP1();
|
|
mtz.set_cell_for_all(unitCell);
|
|
|
|
// Producing-software provenance in the MTZ header (title + HISTORY, the CCP4 convention).
|
|
mtz.title = "Rugnux merged reflections";
|
|
mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
|
|
|
|
// Add dataset
|
|
gemmi::Mtz::Dataset& ds = mtz.add_dataset("native");
|
|
ds.crystal_name = experiment.GetSampleName();
|
|
ds.wavelength = experiment.GetWavelength_A();
|
|
|
|
const int dataset_id = ds.id;
|
|
|
|
// One row per reflection in the CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
|
|
// the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
|
|
bool has_anom = true;
|
|
const std::vector<MergedOutRow> out_rows = BuildMergedRows(reflections, experiment, has_anom);
|
|
|
|
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);
|
|
if (has_anom) {
|
|
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
|
|
mtz.add_column("SIGF", 'Q', dataset_id, -1, false);
|
|
if (has_anom) {
|
|
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);
|
|
|
|
mtz.nreflections = static_cast<int>(out_rows.size());
|
|
mtz.data.reserve(out_rows.size() * (has_anom ? 16 : 8));
|
|
for (const auto& row : out_rows) {
|
|
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(row.Imean);
|
|
mtz.data.push_back(row.sImean);
|
|
if (has_anom) {
|
|
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(row.Fmean);
|
|
mtz.data.push_back(row.sFmean);
|
|
if (has_anom) {
|
|
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));
|
|
}
|
|
mtz.write_to_file(filename);
|
|
}
|
|
|
|
void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
|
|
const DiffractionExperiment &experiment,
|
|
const std::string &filename,
|
|
size_t nthreads) {
|
|
bool has_anom = true;
|
|
const std::vector<MergedOutRow> rows = BuildMergedRows(reflections, experiment, has_anom);
|
|
|
|
// SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as
|
|
// h k l I sigma(I). The Bijvoet mates are written separately - I(+) at +hkl, I(-) at -hkl - so the
|
|
// anomalous differences survive; a reflection with no anomalous split is written once as its mean.
|
|
// Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale
|
|
// is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the
|
|
// anomalous signal is untouched. The file ends with a 0 0 0 terminator record.
|
|
const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; };
|
|
|
|
double max_abs = 0.0;
|
|
for (const auto& r : rows) {
|
|
if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)});
|
|
if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)});
|
|
if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean))
|
|
max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)});
|
|
}
|
|
const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0;
|
|
|
|
std::ofstream out(filename);
|
|
if (!out)
|
|
throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename);
|
|
// Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker
|
|
// blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in
|
|
// the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write.
|
|
const auto column = [](std::string &s, const std::string &v, size_t width) {
|
|
if (v.size() < width) s.append(width - v.size(), ' ');
|
|
s.append(v);
|
|
};
|
|
const auto num2 = [](double v) {
|
|
char b[64];
|
|
const int n = std::snprintf(b, sizeof b, "%.2f", v);
|
|
return std::string(b, static_cast<size_t>(std::clamp(n, 0, static_cast<int>(sizeof b) - 1)));
|
|
};
|
|
const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) {
|
|
column(s, std::to_string(h), 4);
|
|
column(s, std::to_string(k), 4);
|
|
column(s, std::to_string(l), 4);
|
|
column(s, num2(scale * I), 8);
|
|
column(s, num2(scale * sigma), 8);
|
|
s.push_back('\n');
|
|
};
|
|
{
|
|
const size_t nrow = rows.size();
|
|
const size_t nw = std::max<size_t>(nthreads, 1);
|
|
const int nch = static_cast<int>(ThreadsForWork(nrow, nw, 4096));
|
|
std::vector<std::string> block(nch);
|
|
ParallelChunks(nch, nw, [&](int tlo, int thi) {
|
|
for (int t = tlo; t < thi; ++t) {
|
|
const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch;
|
|
std::string &s = block[t];
|
|
s.reserve((hi - lo) * 2 * 29);
|
|
for (size_t i = lo; i < hi; ++i) {
|
|
const auto &r = rows[i];
|
|
const bool plus = usable(r.Ip, r.sIp);
|
|
const bool minus = usable(r.Im, r.sIm);
|
|
if (plus) emit(s, r.h, r.k, r.l, r.Ip, r.sIp);
|
|
if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm);
|
|
if (!plus && !minus && usable(r.Imean, r.sImean))
|
|
emit(s, r.h, r.k, r.l, r.Imean, r.sImean);
|
|
}
|
|
}
|
|
});
|
|
for (const std::string &s : block)
|
|
out.write(s.data(), static_cast<std::streamsize>(s.size()));
|
|
}
|
|
std::string tail;
|
|
emit(tail, 0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker
|
|
out.write(tail.data(), static_cast<std::streamsize>(tail.size()));
|
|
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.
|
|
// The floor is the SAME pair the 3D combine applies to the same event - --min-partiality (default
|
|
// 0.02) and --min-captured-fraction (default 0.7 on rotation) - so what is exported as a full is
|
|
// exactly what rugnux merges as a full. The captured-fraction floor is what makes the FULL
|
|
// declaration true: the batch header says LDTYPE=2 and M/ISYM carries M=0, AIMLESS reads that as
|
|
// "all runs have only fulls" and ignores FRACTIONCALC entirely, so an event that caught a twentieth
|
|
// of its rocking curve does not enter the merge at a twentieth of its intensity - it enters whole,
|
|
// with a small sigma, and is weighted heavily. There is no cut for the reader to make: the column
|
|
// that would let it make one is one the reading program has been told to ignore. Measured, dropping
|
|
// those events moves AIMLESS's Rmerge at 1.8 A from 1.353 to 0.694 and CC(1/2) from 0.985 to 0.993.
|
|
//
|
|
// FRACTIONCALC itself is written as summed, above 1 included. Each part's partiality is the erf pair
|
|
// the predictor computed on THAT frame, from that frame's own refined lattice and its own mosaicity
|
|
// estimate, so the parts of one event do not tile the rocking curve exactly - measured, the offset
|
|
// from one part to the next steps by the wedge to within 12% of it, and the sums land in a peak at
|
|
// 1.000 whose 95th percentile is 1.09. It is the honest estimate of a captured fraction and it is
|
|
// not the number rugnux scales on (the merge recomputes every partiality from a frame-order-smoothed
|
|
// mosaicity, RotationScaleMerge::SmoothMosaicity), so clamping it to 1 would hide the spread and buy
|
|
// a reading program nothing.
|
|
std::vector<Reflection> SumRockingEvents(const std::vector<IntegrationOutcome> &outcomes,
|
|
double min_partiality, double min_captured_fraction) {
|
|
constexpr float MAX_FRAME_GAP = 2.0f; // == RotationScaleMerge's: what makes one rocking event
|
|
|
|
// The sort key travels with the part instead of being read back through the pointer, the way the
|
|
// merge's own ingest sort carries it (RotationScaleMerge's SortKey): there are millions of parts
|
|
// and an indirect compare is a cache miss on every one of them. The keys are the same values in
|
|
// the same order, so introsort makes the same comparisons and the same swaps and leaves the same
|
|
// order - which matters, because two parts can genuinely share (h,k,l) and image_number and the
|
|
// event sums below are floating point.
|
|
struct Part {
|
|
int32_t h, k, l;
|
|
float image_number;
|
|
const Reflection *r;
|
|
};
|
|
std::vector<Part> parts;
|
|
size_t n_parts = 0;
|
|
for (const auto &outcome : outcomes)
|
|
n_parts += outcome.reflections.size();
|
|
parts.reserve(n_parts);
|
|
for (const auto &outcome : outcomes)
|
|
for (const auto &r : outcome.reflections)
|
|
parts.push_back({r.h, r.k, r.l, r.image_number, &r});
|
|
std::sort(parts.begin(), parts.end(), [](const Part &a, const Part &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].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;
|
|
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].r;
|
|
i = j;
|
|
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));
|
|
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 = &experiment.GetSpaceGroupOrP1();
|
|
mtz.set_cell_for_all(unitCell);
|
|
mtz.title = "Rugnux unmerged observations";
|
|
mtz.history.push_back("From Rugnux " + jfjoch_version() + ", data reduction");
|
|
mtz.add_base(); // the HKL_base dataset and the H K L columns
|
|
gemmi::Mtz::Dataset &ds = mtz.add_dataset("unmerged");
|
|
ds.crystal_name = experiment.GetSampleName();
|
|
ds.wavelength = experiment.GetWavelength_A();
|
|
|
|
// Every data column belongs to HKL_base and only the batches carry a dataset id, as in a
|
|
// POINTLESS file. Down to FLAG this is POINTLESS's own column set; the four after it are what
|
|
// rugnux measures beyond it - the offset of the reflection from the centre of its rocking curve,
|
|
// the Lorentz geometry of that curve, and the background that was subtracted.
|
|
mtz.add_column("M/ISYM", 'Y', 0, -1, false);
|
|
mtz.add_column("BATCH", 'B', 0, -1, false);
|
|
mtz.add_column("I", 'J', 0, -1, false);
|
|
mtz.add_column("SIGI", 'Q', 0, -1, false);
|
|
mtz.add_column("FRACTIONCALC", 'R', 0, -1, false);
|
|
mtz.add_column("XDET", 'R', 0, -1, false);
|
|
mtz.add_column("YDET", 'R', 0, -1, false);
|
|
mtz.add_column("ROT", 'R', 0, -1, false);
|
|
mtz.add_column("LP", 'R', 0, -1, false);
|
|
mtz.add_column("FLAG", 'I', 0, -1, false);
|
|
mtz.add_column("DELPHI", 'R', 0, -1, false);
|
|
mtz.add_column("ZETA", 'R', 0, -1, false);
|
|
mtz.add_column("BGMEAN", 'R', 0, -1, false);
|
|
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 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;
|
|
};
|
|
|
|
// H K L are the ASU indices and M/ISYM says which symmetry operation (and which Friedel hand) got
|
|
// 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.
|
|
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
|
|
// indexed nothing leaves a hole that cuts the sweep in two: measured on a 360 deg sweep with 67
|
|
// unindexed images, AIMLESS made 23 runs of it and its scale model diverged - CC(1/2) 0.811,
|
|
// Rmeas overflowed to -1266, and the result no longer correlated with rugnux's own merge
|
|
// (Pearson 0.0018). Filling the gaps gives one run and the numbers the sweep deserves.
|
|
int first_batch = 0, last_batch = -1;
|
|
|
|
// Lattice-centring absences are integrated on purpose - prediction runs in P so the space-group
|
|
// search can confirm or disprove the centring - but they are not written here. POINTLESS reads
|
|
// the interstitial nodes, takes the lattice for primitive and demotes the space group, and they
|
|
// roughly halve AIMLESS's <I/sigma>. Screw and glide absences ARE written: they are the evidence
|
|
// the space group was chosen on, so a program reading this file can check that choice, and
|
|
// removing them turns its test into an assumption. XDS and DIALS draw the line in the same place.
|
|
const char centering = mtz.spacegroup ? mtz.spacegroup->hm[0] : 'P';
|
|
|
|
const auto add_row = [&](const Reflection &r) {
|
|
if (systematic_absence(r.h, r.k, r.l, centering))
|
|
return;
|
|
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));
|
|
if (last_batch < first_batch) { first_batch = batch; last_batch = batch; }
|
|
first_batch = std::min(first_batch, batch);
|
|
last_batch = std::max(last_batch, 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(),
|
|
experiment.GetScalingSettings().GetMinCapturedFraction()))
|
|
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());
|
|
|
|
// The batch header's orientation matrix is the crystal at rotation angle zero - each batch's own
|
|
// PHISTT is applied on top of it - but the lattice stored with an outcome is the crystal as it
|
|
// stood on that image. Turn the first indexed one back by its own angle to get the orientation of
|
|
// the sweep, which is the one matrix POINTLESS also writes into every batch.
|
|
std::optional<CrystalLattice> lattice_at_zero;
|
|
std::optional<float> mosaicity_deg;
|
|
for (const auto &outcome : outcomes) {
|
|
if (outcome.reflections.empty() || outcome.latt.CalcVolume() <= 1.0f)
|
|
continue;
|
|
const float mid_deg = phi_start_deg(outcome.reflections.front().image_number) + wedge_deg / 2.0f;
|
|
lattice_at_zero = gon ? outcome.latt.Multiply(gon->GetTransformationAngle(mid_deg)) : outcome.latt;
|
|
mosaicity_deg = outcome.mosaicity_deg;
|
|
break;
|
|
}
|
|
|
|
// The batch header is written in the "Cambridge" frame - z along the principal rotation axis, x
|
|
// along the beam - while the jfjoch lab frame has the beam along +z, so the two are related by a
|
|
// rotation. These three lab-frame vectors are the Cambridge axes; a lab vector's components in
|
|
// that frame are its dot products with them. A still has no rotation axis, and any axis across
|
|
// the beam then defines the frame just as consistently.
|
|
// The axis is NEGATED: rugnux turns an observation made at angle phi back to phi = 0 by rotating
|
|
// it by +phi about the goniometer axis, so the crystal itself turns by -phi about it, while the
|
|
// MTZ batch header's scan axis is the one a batch's own increasing PHI turns the crystal about.
|
|
const Coord beam = experiment.GetDiffractionGeometry().GetScatteringVector().Normalize();
|
|
const Coord z_cam = gon ? -gon->GetAxis().Normalize() : Coord(0, 1, 0);
|
|
const Coord y_cam = (z_cam % beam).Normalize();
|
|
const Coord x_cam = (y_cam % z_cam).Normalize();
|
|
|
|
gemmi::Mtz::Batch batch;
|
|
batch.title = "Rugnux";
|
|
batch.axes.emplace_back("PHI");
|
|
batch.set_dataset_id(ds.id);
|
|
batch.ints[12] = 1; // ncryst
|
|
batch.ints[14] = partials ? 1 : 2; // ldtype: oscillation data (2D spots) / area detector (3D)
|
|
batch.ints[15] = 1; // jsaxs: the goniostat scan axis
|
|
batch.ints[17] = 1; // ngonax
|
|
batch.ints[19] = 1; // ndet
|
|
batch.set_cell(unitCell);
|
|
if (lattice_at_zero) {
|
|
// Orientation matrix U, built from the reciprocal axes and stored column by column in
|
|
// Cambridge components, as gemmi's XDS_ASCII converter builds it (gemmi/xds2mtz.hpp).
|
|
const Coord ar = lattice_at_zero->Astar().Normalize();
|
|
const Coord cr = (ar % lattice_at_zero->Bstar()).Normalize();
|
|
const Coord u[3] = {ar, cr % ar, cr};
|
|
for (int i = 0; i < 3; i++) {
|
|
batch.floats[6 + 3 * i] = u[i] * x_cam;
|
|
batch.floats[7 + 3 * i] = u[i] * y_cam;
|
|
batch.floats[8 + 3 * i] = u[i] * z_cam;
|
|
}
|
|
}
|
|
batch.floats[21] = mosaicity_deg.value_or(0.0f); // crydat(0), the reflecting range
|
|
batch.floats[40] = 1.0f; // scanax = [0, 0, 1]: the rotation axis IS z in the Cambridge frame
|
|
batch.floats[47] = wedge_deg;
|
|
batch.floats[61] = 1.0f; // e1 = scanax, the only goniostat axis
|
|
batch.floats[80] = -1.0f; // idealised source vector, antiparallel to the beam
|
|
batch.floats[83] = -(beam * x_cam); // s0, the source vector of the geometry as it really stands
|
|
batch.floats[84] = -(beam * y_cam);
|
|
batch.floats[85] = -(beam * z_cam);
|
|
batch.set_wavelength(experiment.GetWavelength_A());
|
|
batch.floats[111] = experiment.GetDetectorDistance_mm();
|
|
batch.floats[113] = 1.0f; // detector limits, in pixels
|
|
batch.floats[114] = static_cast<float>(experiment.GetXPixelsNum());
|
|
batch.floats[115] = 1.0f;
|
|
batch.floats[116] = static_cast<float>(experiment.GetYPixelsNum());
|
|
|
|
for (int number = first_batch; number <= last_batch; number++) {
|
|
batch.number = number;
|
|
batch.floats[36] = phi_start_deg(static_cast<float>(number - 1)); // phistt
|
|
batch.floats[37] = batch.floats[36] + wedge_deg; // phiend
|
|
mtz.batches.push_back(batch);
|
|
}
|
|
|
|
mtz.sort(5); // by H K L M/ISYM BATCH, the order POINTLESS leaves an unmerged file in
|
|
mtz.write_to_file(filename);
|
|
}
|
|
|
|
void WriteReflections(const std::vector<MergedReflection> &reflections,
|
|
const UnitCell &unitCell,
|
|
const DiffractionExperiment &experiment,
|
|
const MergeStatistics &statistics,
|
|
const ErrorModelReport &error_model,
|
|
const TwinningAnalysisResult &twinning,
|
|
const std::string &filename,
|
|
size_t nthreads) {
|
|
// Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 /
|
|
// phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as
|
|
// the SHELXC / ANODE substructure-solution input).
|
|
WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz");
|
|
WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning,
|
|
filename + ".cif", nthreads);
|
|
WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads);
|
|
}
|