From 5405ee929e0c16000a362cccd7b6c080e6c7a8d9 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Tue, 3 Mar 2026 14:44:54 +0100
Subject: [PATCH] jfjoch_process: Add unit cell parameter

---
 tools/jfjoch_process.cpp | 79 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 78 insertions(+), 1 deletion(-)

diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp
index db751f45..6d86ac3f 100644
--- a/tools/jfjoch_process.cpp
+++ b/tools/jfjoch_process.cpp
@@ -44,8 +44,68 @@ void print_usage(Logger &logger) {
     logger.Info("   -M         Scale and merge (refine mosaicity) and write scaled.hkl + image.dat");
     logger.Info("   -P<txt>    Partiality refinement fixed|rot|unity (default: fixed)");
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
+    logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
 }
 
+void trim_in_place(std::string& t) {
+    size_t b = 0;
+    while (b < t.size() && std::isspace(static_cast<unsigned char>(t[b]))) b++;
+    size_t e = t.size();
+    while (e > b && std::isspace(static_cast<unsigned char>(t[e - 1]))) e--;
+    t = t.substr(b, e - b);
+};
+
+std::optional<UnitCell> parse_unit_cell_arg(const char* arg)  {
+    if (!arg)
+        return std::nullopt;
+
+    std::string s(arg);
+
+
+    trim_in_place(s);
+
+    if (s.size() >= 2 && ((s.front() == '"' && s.back() == '"') || (s.front() == '\'' && s.back() == '\''))) {
+        s = s.substr(1, s.size() - 2);
+        trim_in_place(s);
+    }
+
+    std::vector<std::string> parts;
+    parts.reserve(6);
+    size_t start = 0;
+    while (true) {
+        size_t pos = s.find(',', start);
+        if (pos == std::string::npos) {
+            parts.push_back(s.substr(start));
+            break;
+        }
+        parts.push_back(s.substr(start, pos - start));
+        start = pos + 1;
+    }
+
+    if (parts.size() != 6)
+        return std::nullopt;
+
+    auto parse_float_strict = [](const std::string& t, float& out) -> bool {
+        try {
+            size_t idx = 0;
+            out = std::stof(t, &idx);
+            return idx == t.size();
+        } catch (...) {
+            return false;
+        }
+    };
+
+    UnitCell uc{};
+    if (!parse_float_strict(parts[0], uc.a)) return std::nullopt;
+    if (!parse_float_strict(parts[1], uc.b)) return std::nullopt;
+    if (!parse_float_strict(parts[2], uc.c)) return std::nullopt;
+    if (!parse_float_strict(parts[3], uc.alpha)) return std::nullopt;
+    if (!parse_float_strict(parts[4], uc.beta)) return std::nullopt;
+    if (!parse_float_strict(parts[5], uc.gamma)) return std::nullopt;
+
+    return uc;
+};
+
 int main(int argc, char **argv) {
     RegisterHDF5Filter();
 
@@ -66,6 +126,7 @@ int main(int argc, char **argv) {
     bool run_scaling = false;
     bool anomalous_mode = false;
     std::optional<int> space_group_number;
+    std::optional<UnitCell> fixed_reference_unit_cell;
 
     ScaleMergeOptions::PartialityModel partiality_model = ScaleMergeOptions::PartialityModel::Fixed;
 
@@ -78,7 +139,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MP:AD:")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MP:AD:C:")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -121,6 +182,18 @@ int main(int argc, char **argv) {
             case 'A':
                 anomalous_mode = true;
                 break;
+            case 'C': {
+                auto uc = parse_unit_cell_arg(optarg);
+                if (!uc.has_value()) {
+                    logger.Error("Invalid -C unit cell. Expected: -C\"a,b,c,alpha,beta,gamma\" (6 floats, comma-separated, no spaces). Got: {}", optarg ? optarg : "<null>");
+                    print_usage(logger);
+                    exit(EXIT_FAILURE);
+                }
+                fixed_reference_unit_cell = uc;
+                logger.Info("Fixed reference unit cell set: a={:.3f} b={:.3f} c={:.3f} alpha={:.3f} beta={:.3f} gamma={:.3f}",
+                            uc->a, uc->b, uc->c, uc->alpha, uc->beta, uc->gamma);
+                break;
+            }
             case 'P':
                 if (strcmp(optarg, "unity") == 0)
                     partiality_model = ScaleMergeOptions::PartialityModel::Unity;
@@ -188,6 +261,10 @@ int main(int argc, char **argv) {
     experiment.OverwriteExistingFiles(true);
     experiment.PolarizationFactor(0.99);
 
+    if (fixed_reference_unit_cell.has_value()) {
+        experiment.SetUnitCell(*fixed_reference_unit_cell);
+    }
+
     // Configure Indexing
     IndexingSettings indexing_settings;
     if (use_fft)