From 1a70c1987a39d4aa15985bfaa31653d975ab196a Mon Sep 17 00:00:00 2001
From: leonarski_f <filip.leonarski@psi.ch>
Date: Mon, 18 May 2026 21:28:09 +0200
Subject: [PATCH] jfjoch_scale: A bit nicer handling of reflection statistics
 and preparation

---
 image_analysis/UpdateReflectionResolution.cpp | 22 +++++++++++---
 image_analysis/UpdateReflectionResolution.h   | 12 +++++++-
 tools/jfjoch_scale.cpp                        | 30 ++++++-------------
 3 files changed, 38 insertions(+), 26 deletions(-)

diff --git a/image_analysis/UpdateReflectionResolution.cpp b/image_analysis/UpdateReflectionResolution.cpp
index b803c70e..456d56f3 100644
--- a/image_analysis/UpdateReflectionResolution.cpp
+++ b/image_analysis/UpdateReflectionResolution.cpp
@@ -4,20 +4,34 @@
 #include "UpdateReflectionResolution.h"
 #include "../common/CrystalLattice.h"
 
-void UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections) {
+ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections) {
+    ResolutionStats ret;
+
     CrystalLattice latt(cell);
     const auto astar = latt.Astar();
     const auto bstar = latt.Bstar();
     const auto cstar = latt.Cstar();
 
     for (auto &image: reflections) {
+        if (!image.empty()) {
+            ret.n_images++;
+            ret.n_reflections += image.size();
+        }
+
         for (auto &r: image) {
             Coord q = r.h * astar + r.k * bstar + r.l * cstar;
             auto qlen = q.Length();
-            if (qlen > 1e-6)
-                r.d = 1/qlen;
-            else
+            if (qlen > 1e-6) {
+                const float d = 1/qlen;
+                if (ret.d_high > d)
+                    ret.d_high = d;
+                if (ret.d_low < d)
+                    ret.d_low = d;
+                r.d = d;
+                r.image_scale_corr = r.rlp / r.partiality;
+            } else
                 r.d = NAN;
         }
     }
+    return ret;
 }
diff --git a/image_analysis/UpdateReflectionResolution.h b/image_analysis/UpdateReflectionResolution.h
index 6c3c26ad..53c26e04 100644
--- a/image_analysis/UpdateReflectionResolution.h
+++ b/image_analysis/UpdateReflectionResolution.h
@@ -1,8 +1,18 @@
 // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
+#pragma once
+
 #include <vector>
 #include "../common/Reflection.h"
 #include "../common/UnitCell.h"
 
-void UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections);
+struct ResolutionStats {
+    float d_low = std::numeric_limits<float>::max();
+    float d_high = 0.0f;
+    int n_reflections = 0;
+    int n_images = 0;
+};
+
+// Returns the highest (best) resolution in Angstrom observed by reflection
+ResolutionStats UpdateReflectionResolution(const UnitCell &cell, std::vector<std::vector<Reflection> > &reflections);
diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp
index c3831097..8a3678ff 100644
--- a/tools/jfjoch_scale.cpp
+++ b/tools/jfjoch_scale.cpp
@@ -240,31 +240,19 @@ int main(int argc, char **argv) {
 
     experiment.ImportScalingSettings(scaling_settings);
 
-    logger.Info("Running scaling (mosaicity refinement) ...");
-    auto reflections = reader.ReadReflections(start_image, end_image);
-    size_t nimages = 0;
-    size_t nrefl = 0;
-    for (auto& image : reflections) {
-        if (!image.empty()) {
-            nimages++;
-            nrefl += image.size();
-
-            for (auto &r: image) {
-                r.image_scale_corr = 1.0; // restart these numbers, since we don't know any better
-            }
-        }
-    }
-
-    logger.Info("Read {} reflections from {} images", nrefl, nimages);
-
-    if (experiment.GetUnitCell()) {
-        UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections);
-        logger.Info("Reflection resolution updated based on experiment unit cell");
-    } else {
+    if (!experiment.GetUnitCell()) {
         logger.Error("Experiment unit cell not found, cannot update reflection resolution");
         exit(EXIT_FAILURE);
     }
 
+    logger.Info("Loading reflections");
+
+    auto reflections = reader.ReadReflections(start_image, end_image);
+
+    auto refl_stats = UpdateReflectionResolution(experiment.GetUnitCell().value(), reflections);
+
+    logger.Info("Read {} reflections from {} images", refl_stats.n_reflections, refl_stats.n_images);
+    
     std::vector<float> mosaicity(end_image - start_image + 1);
     for (int i = 0; i < end_image - start_image + 1; i++) {
         mosaicity[i] = dataset->mosaicity_deg[start_image + i];