diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 6913b151..826d5ba7 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "WriteReflections.h" +#include "scale_merge/Merge.h" #include #include @@ -59,6 +60,8 @@ std::string CifStr(const std::string& s) { void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, + const MergeStatistics &statistics, + const std::string &isa, const std::string &filename) { std::ofstream out(filename); @@ -110,6 +113,51 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; out << "_diffrn_detector.detector " << CifStr(experiment.GetDetectorDescription()) << "\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) has no standard CIF + // item, so it is written as a free-text pdbx value. + const auto mult = [](const MergeStatisticsShell &s) { + return s.unique_reflections > 0 ? static_cast(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(s.unique_reflections) / s.possible_unique_reflections : 0.0; }; + if (!statistics.shells.empty()) { + const auto &ov = statistics.overall; + 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"; + out << "_reflns.pdbx_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\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"; + 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) << "\n"; + } + out << "#\n"; + } // ---------- _refln loop ---------- out << "loop_\n"; @@ -199,13 +247,15 @@ void WriteHKLReflections(const std::vector &reflections, void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, + const MergeStatistics &statistics, + const std::string &isa, const std::string &filename) { switch (experiment.GetScalingSettings().GetFileFormat()) { case IntensityFormat::Text: WriteHKLReflections(reflections, filename + ".hkl"); break; case IntensityFormat::mmCIF: - WriteMmcifReflections(reflections, unitCell, experiment, filename + ".cif"); + WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, filename + ".cif"); break; case IntensityFormat::MTZ: WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 9299879d..8215dd93 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -10,9 +10,13 @@ #include "../common/UnitCell.h" #include "../common/DiffractionExperiment.h" +struct MergeStatistics; + void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, + const MergeStatistics &statistics, + const std::string &isa, const std::string &filename); void WriteMtzReflections(const std::vector &reflections, @@ -26,4 +30,6 @@ void WriteHKLReflections(const std::vector &reflections, void WriteReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, + const MergeStatistics &statistics, + const std::string &isa, const std::string &filename); \ No newline at end of file diff --git a/process/JFJochProcess.cpp b/process/JFJochProcess.cpp index e0976638..e206ea18 100644 --- a/process/JFJochProcess.cpp +++ b/process/JFJochProcess.cpp @@ -671,7 +671,9 @@ ProcessResult JFJochProcess::Run(JFJochProcessObserver *observer) { if (result.consensus_cell && write_files) { phase("Writing reflections"); - WriteReflections(sm.merged, *result.consensus_cell, experiment_, config_.output_prefix); + WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, + result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", + config_.output_prefix); // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS // comparison. The offline self-scaling result is otherwise not exposed (process.h5's // per-image arrays are only filled on the online per-image path). Sourced from the diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 4ce1e338..1f992e8a 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -406,5 +406,7 @@ int main(int argc, char **argv) { std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl; if (!output_prefix.empty()) - WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, output_prefix); + WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, + merge_engine.ErrorModelB() > 0 ? fmt::format("{:.2f}", 1.0 / merge_engine.ErrorModelB()) : "?", + output_prefix); }