From 39823ddb67fdc52589e2fc80378e43b4fd333bfb Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Tue, 19 May 2026 22:12:42 +0200 Subject: [PATCH] General function to write reflections in multiple formats (need to better integrate into jfjoch_scale and jfjoch_process) --- image_analysis/CMakeLists.txt | 4 +- image_analysis/WriteMmcif.cpp | 179 ----------------------- image_analysis/WriteMmcif.h | 50 ------- image_analysis/WriteReflections.cpp | 214 ++++++++++++++++++++++++++++ image_analysis/WriteReflections.h | 29 ++++ tools/jfjoch_process.cpp | 37 +---- tools/jfjoch_scale.cpp | 33 +---- 7 files changed, 251 insertions(+), 295 deletions(-) delete mode 100644 image_analysis/WriteMmcif.cpp delete mode 100644 image_analysis/WriteMmcif.h create mode 100644 image_analysis/WriteReflections.cpp create mode 100644 image_analysis/WriteReflections.h diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 40bd209f..66537115 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -29,8 +29,8 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC RotationIndexer.h RotationParameters.cpp RotationParameters.h - WriteMmcif.cpp - WriteMmcif.h + WriteReflections.cpp + WriteReflections.h LoadFCalcFromMtz.cpp LoadFCalcFromMtz.h UpdateReflectionResolution.cpp diff --git a/image_analysis/WriteMmcif.cpp b/image_analysis/WriteMmcif.cpp deleted file mode 100644 index db5aa29b..00000000 --- a/image_analysis/WriteMmcif.cpp +++ /dev/null @@ -1,179 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#include "WriteMmcif.h" - -#include -#include -#include -#include -#include -#include -#include - -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. -std::string Fmt(double val, int decimals = 4) { - if (!std::isfinite(val)) - return "?"; - std::ostringstream ss; - ss << std::fixed << std::setprecision(decimals) << val; - return ss.str(); -} - -/// 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; -} - -} // namespace - -void WriteMmcifReflections(std::ostream& out, - const std::vector& reflections, - const MmcifMetadata& meta) { - - out << std::fixed; - - // ---------- data block ---------- - out << "data_" << meta.data_block_name << "\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 'Jungfraujoch'\n"; - if (meta.software_version.has_value()) - out << "_software.version " << CifStr(meta.software_version.value()) << "\n"; - out << "_software.classification reduction\n"; - out << "#\n"; - - // ---------- _cell ---------- - out << "_cell.length_a " << Fmt(meta.unit_cell.a, 3) << "\n"; - out << "_cell.length_b " << Fmt(meta.unit_cell.b, 3) << "\n"; - out << "_cell.length_c " << Fmt(meta.unit_cell.c, 3) << "\n"; - out << "_cell.angle_alpha " << Fmt(meta.unit_cell.alpha, 2) << "\n"; - out << "_cell.angle_beta " << Fmt(meta.unit_cell.beta, 2) << "\n"; - out << "_cell.angle_gamma " << Fmt(meta.unit_cell.gamma, 2) << "\n"; - - // ---------- _symmetry ---------- - out << "_symmetry.space_group_name_H-M " << CifStr(meta.space_group_name) << "\n"; - out << "_symmetry.Int_Tables_number " << meta.space_group_number << "\n"; - out << "#\n"; - - // ---------- _diffrn_source / _diffrn_detector ---------- - if (meta.source.has_value()) { - out << "_diffrn_source.pdbx_synchrotron_site " - << CifStr(meta.source.value()) << "\n"; - } - - if (meta.beamline.has_value()) { - out << "_diffrn_source.pdbx_synchrotron_beamline " - << CifStr(meta.beamline.value()) << "\n"; - } - - if (meta.wavelength_A.has_value()) { - out << "_diffrn_radiation_wavelength.wavelength " - << Fmt(meta.wavelength_A.value(), 5) << "\n"; - } - - if (!meta.detector_name.empty()) { - out << "_diffrn_detector.detector " << CifStr(meta.detector_name) << "\n"; - } - - if (meta.detector_distance_mm.has_value()) { - out << "_diffrn_detector.distance " - << Fmt(meta.detector_distance_mm.value(), 2) << "\n"; - } - - if (meta.sample_temperature_K.has_value()) { - out << "_diffrn.ambient_temp " - << Fmt(meta.sample_temperature_K.value(), 1) << "\n"; - } - - if (meta.sample_name.has_value() && !meta.sample_name->empty()) { - out << "_entity.id 1\n"; - out << "_entity.type polymer\n"; - out << "_entity.pdbx_description " << CifStr(meta.sample_name.value()) << "\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.status\n"; - - for (const auto& r : reflections) { - out << std::setw(5) << r.h << " " - << std::setw(5) << r.k << " " - << std::setw(5) << r.l << " " - << std::setw(14) << Fmt(r.I, 4) << " " - << std::setw(14) << Fmt(r.sigma, 4) << " " - << "o" // 'o' = observed - << "\n"; - } - - out << "#\n"; - out << "# End of reflections\n"; -} - -void WriteMmcifReflections(const std::string& path, - const std::vector& reflections, - const MmcifMetadata& meta) { - std::ofstream file(path); - if (!file) - throw std::runtime_error("WriteMmcifReflections: cannot open " + path); - WriteMmcifReflections(file, reflections, meta); - file.close(); - if (!file) - throw std::runtime_error("WriteMmcifReflections: I/O error writing " + path); -} - -void MmcifMetadata::Fill(const DiffractionExperiment &experiment) { - if (experiment.GetUnitCell().has_value()) - unit_cell = experiment.GetUnitCell().value(); - - auto sg = experiment.GetGemmiSpaceGroup(); - if (sg.has_value()) { - space_group_name = sg->hm; - space_group_number = sg->number; - } - - detector_name = experiment.GetDetectorDescription(); - wavelength_A = experiment.GetWavelength_A(); - detector_distance_mm = experiment.GetDetectorDistance_mm(); - sample_temperature_K = experiment.GetSampleTemperature_K(); - sample_name = experiment.GetSampleName(); - - beamline = experiment.GetInstrumentName(); - source = experiment.GetSourceName(); -} diff --git a/image_analysis/WriteMmcif.h b/image_analysis/WriteMmcif.h deleted file mode 100644 index ade273fc..00000000 --- a/image_analysis/WriteMmcif.h +++ /dev/null @@ -1,50 +0,0 @@ -// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute -// SPDX-License-Identifier: GPL-3.0-only - -#pragma once - -#include -#include -#include -#include - -#include "../common/Reflection.h" -#include "../common/UnitCell.h" -#include "../common/DiffractionExperiment.h" - -/// Metadata needed to write a meaningful mmCIF reflection file. -struct MmcifMetadata { - // Required - UnitCell unit_cell{}; - std::string space_group_name; // e.g. "P 21 21 21" - int space_group_number = 1; - - // Optional but recommended - std::string data_block_name = "jfjoch"; // CIF data_ - std::string detector_name; // e.g. "JUNGFRAU 4M" - std::optional wavelength_A; // incident wavelength - std::optional detector_distance_mm; - std::optional sample_temperature_K; - std::optional sample_name; - std::optional software_version; // jfjoch version string - - std::optional source; - std::optional beamline; - - void Fill(const DiffractionExperiment& experiment); -}; - -/// Write a PDBx/mmCIF reflection file (containing _refln loop with F/sigmaF/I/sigmaI) -/// to the given output stream. -/// -/// The file follows the conventions expected by CCP4/CCTBX/Phenix for -/// structure-factor mmCIF files. -void WriteMmcifReflections(std::ostream& out, - const std::vector& reflections, - const MmcifMetadata& meta); - -/// Convenience overload that writes to a file path. -/// Throws std::runtime_error on I/O failure. -void WriteMmcifReflections(const std::string& path, - const std::vector& reflections, - const MmcifMetadata& meta); \ No newline at end of file diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp new file mode 100644 index 00000000..6913b151 --- /dev/null +++ b/image_analysis/WriteReflections.cpp @@ -0,0 +1,214 @@ +// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "WriteReflections.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#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. +std::string Fmt(double val, int decimals = 4) { + if (!std::isfinite(val)) + return "?"; + std::ostringstream ss; + ss << std::fixed << std::setprecision(decimals) << val; + return ss.str(); +} + +/// 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; +} + +} // namespace + + +void WriteMmcifReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename) { + + 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 'Jungfraujoch'\n"; + + out << "_software.version " << CifStr(jfjoch_version()) << "\n"; + out << "_software.classification 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"; + + auto *sg = gemmi::find_spacegroup_by_number(experiment.GetSpaceGroupNumber().value_or(1)); + if (sg == nullptr) + throw std::runtime_error("WriteMmcifReflections: invalid space group number"); + + // ---------- _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"; + + // ---------- _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.status_free\n"; + out << "_refln.status\n"; + + for (const auto& r : reflections) { + out << std::setw(5) << r.h << " " + << std::setw(5) << r.k << " " + << std::setw(5) << r.l << " " + << std::setw(14) << Fmt(r.I, 4) << " " + << std::setw(14) << Fmt(r.sigma, 4) << " " + << (r.rfree_flag ? 1 : 0) << " " + << "o" // 'o' = observed + << "\n"; + } + + out << "#\n"; + out << "# End of reflections\n"; + out.close(); +} + +void WriteMtzReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename) { + gemmi::Mtz mtz; + + // Optional but recommended metadata + mtz.spacegroup = gemmi::find_spacegroup_by_number( + experiment.GetSpaceGroupNumber().value_or(1)); + mtz.set_cell_for_all(unitCell); + + // Add dataset + gemmi::Mtz::Dataset& ds = mtz.add_dataset("native"); + ds.crystal_name = experiment.GetSampleName(); + ds.wavelength = experiment.GetWavelength_A(); + + int dataset_id = ds.id; + + // 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("FreeR_flag", 'I', dataset_id, -1, false); + + // Number of rows + mtz.nreflections = reflections.size(); + + // Allocate data table + mtz.data.reserve(reflections.size() * 6); + + for (const auto& r : reflections) { + mtz.data.push_back(static_cast(r.h)); + mtz.data.push_back(static_cast(r.k)); + mtz.data.push_back(static_cast(r.l)); + mtz.data.push_back(r.I); + mtz.data.push_back(r.sigma); + mtz.data.push_back(r.rfree_flag ? 1 : 0); + } + + // Write MTZ + mtz.write_to_file(filename); +} + +void WriteHKLReflections(const std::vector &reflections, + const std::string &filename) { + std::ofstream hkl_file(filename); + if (!hkl_file) + throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot open file " + filename + " for writing"); + + for (const auto &r: reflections) + hkl_file << r.h << " " << r.k << " " << r.l << " " + << r.I << " " << r.sigma << " " + << (r.rfree_flag ? 1 : 0) << std::endl; + + hkl_file.close(); +} + +void WriteReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + 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"); + break; + case IntensityFormat::MTZ: + WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); + break; + } +} diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h new file mode 100644 index 00000000..9299879d --- /dev/null +++ b/image_analysis/WriteReflections.h @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../common/Reflection.h" +#include "../common/UnitCell.h" +#include "../common/DiffractionExperiment.h" + +void WriteMmcifReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename); + +void WriteMtzReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename); + +void WriteHKLReflections(const std::vector &reflections, + const std::string &filename); + +void WriteReflections(const std::vector &reflections, + const UnitCell &unitCell, + const DiffractionExperiment &experiment, + const std::string &filename); \ No newline at end of file diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index ae09db0e..10cc300c 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -28,7 +28,7 @@ #include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" -#include "../image_analysis/WriteMmcif.h" +#include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" void print_usage(Logger &logger) { @@ -719,39 +719,8 @@ int main(int argc, char **argv) { merged_statistics.Print(logger); - { - const std::string hkl_path = output_prefix + "_intensities.hkl"; - std::ofstream hkl_file(hkl_path); - if (!hkl_file) { - logger.Error("Cannot open {} for writing", hkl_path); - } else { - for (const auto &r: merged_reflections) { - hkl_file << r.h << " " << r.k << " " << r.l << " " - << r.I << " " << r.sigma - << "\n"; - } - hkl_file.close(); - logger.Info("Wrote {} reflections to {}", merged_reflections.size(), hkl_path); - } - } - - try { - const std::string cif_path = output_prefix + "_intensities.cif"; - - MmcifMetadata cif_meta; - - cif_meta.Fill(experiment); - - if (rotation_indexer_ret.has_value()) - cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell(); - cif_meta.data_block_name = output_prefix; - - WriteMmcifReflections(cif_path, merged_reflections, cif_meta); - logger.Info("Wrote mmCIF reflections to {}", cif_path); - } catch (const std::exception &e) { - logger.Error("Failed to write mmCIF: {}", e.what()); - } - + if (consensus_cell && !output_prefix.empty()) + WriteReflections(merged_reflections, *consensus_cell, experiment, output_prefix); } // Write End Message diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 8a3678ff..d2b975e5 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -28,7 +28,7 @@ #include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/scale_merge/Merge.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" -#include "../image_analysis/WriteMmcif.h" +#include "../image_analysis/WriteReflections.h" #include "../image_analysis/UpdateReflectionResolution.h" void print_usage(Logger &logger) { @@ -333,33 +333,6 @@ int main(int argc, char **argv) { // Print resolution-shell statistics table merge_stats.Print(logger); - { - const std::string hkl_path = output_prefix + "_intensities.hkl"; - std::ofstream hkl_file(hkl_path); - if (!hkl_file) { - logger.Error("Cannot open {} for writing", hkl_path); - } else { - for (const auto &r: merge_result) { - hkl_file << r.h << " " << r.k << " " << r.l << " " - << r.I << " " << r.sigma - << "\n"; - } - hkl_file.close(); - logger.Info("Wrote {} reflections to {}", merge_result.size(), hkl_path); - } - } - - try { - const std::string cif_path = output_prefix + "_intensities.cif"; - - MmcifMetadata cif_meta; - - cif_meta.Fill(experiment); - cif_meta.data_block_name = output_prefix; - - WriteMmcifReflections(cif_path, merge_result, cif_meta); - logger.Info("Wrote mmCIF reflections to {}", cif_path); - } catch (const std::exception &e) { - logger.Error("Failed to write mmCIF: {}", e.what()); - } + if (!output_prefix.empty()) + WriteReflections(merge_result, *experiment.GetUnitCell(), experiment, output_prefix); }