General function to write reflections in multiple formats (need to better integrate into jfjoch_scale and jfjoch_process)

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

image_analysis/WriteMmcif.cpp

@@ -1,179 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
-// SPDX-License-Identifier: GPL-3.0-only
-
-#include "WriteMmcif.h"
-
-#include <cmath>
-#include <fstream>
-#include <iomanip>
-#include <sstream>
-#include <stdexcept>
-#include <ctime>
-#include <chrono>
-
-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<MergedReflection>& 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<MergedReflection>& 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();
-}

image_analysis/WriteMmcif.h

@@ -1,50 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
-// SPDX-License-Identifier: GPL-3.0-only
-
-#pragma once
-
-#include <ostream>
-#include <string>
-#include <vector>
-#include <optional>
-
-#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_<name>
-    std::string detector_name;                   // e.g. "JUNGFRAU 4M"
-    std::optional<float> wavelength_A;           // incident wavelength
-    std::optional<float> detector_distance_mm;
-    std::optional<float> sample_temperature_K;
-    std::optional<std::string> sample_name;
-    std::optional<std::string> software_version; // jfjoch version string
-
-    std::optional<std::string> source;
-    std::optional<std::string> 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<MergedReflection>& 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<MergedReflection>& reflections,
-                           const MmcifMetadata& meta);

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 <cmath>
+#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.
+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<MergedReflection> &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<MergedReflection> &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<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.rfree_flag ? 1 : 0);
+    }
+
+    // Write MTZ
+    mtz.write_to_file(filename);
+}
+
+void WriteHKLReflections(const std::vector<MergedReflection> &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<MergedReflection> &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;
+    }
+}

image_analysis/WriteReflections.h

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: 2025 Paul Scherrer Institute
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "../common/Reflection.h"
+#include "../common/UnitCell.h"
+#include "../common/DiffractionExperiment.h"
+
+void WriteMmcifReflections(const std::vector<MergedReflection> &reflections,
+                           const UnitCell &unitCell,
+                           const DiffractionExperiment &experiment,
+                           const std::string &filename);
+
+void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
+                         const UnitCell &unitCell,
+                         const DiffractionExperiment &experiment,
+                         const std::string &filename);
+
+void WriteHKLReflections(const std::vector<MergedReflection> &reflections,
+                          const std::string &filename);
+
+void WriteReflections(const std::vector<MergedReflection> &reflections,
+                         const UnitCell &unitCell,
+                         const DiffractionExperiment &experiment,
+                         const std::string &filename);

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

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);
 }