From 9bb91ee84900ca981dce6ecac7abcf8f1fe0743f Mon Sep 17 00:00:00 2001
From: leonarski_f <filip.leonarski@psi.ch>
Date: Sat, 28 Feb 2026 13:06:23 +0100
Subject: [PATCH] ScaleAndMerge: Minor fixes

---
 image_analysis/scale_merge/ScaleAndMerge.cpp | 21 ++++++++++++--------
 1 file changed, 13 insertions(+), 8 deletions(-)

diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp
index bc9b931b..74f339f7 100644
--- a/image_analysis/scale_merge/ScaleAndMerge.cpp
+++ b/image_analysis/scale_merge/ScaleAndMerge.cpp
@@ -158,7 +158,7 @@ namespace {
             : Iobs_(static_cast<double>(r.I)),
               weight_(SafeInv(sigma_obs, 1.0)),
               lp_(SafeInv(r.rlp, 1.0)),
-              dist_ewald_sq_(r.dist_ewald) {
+              dist_ewald_sq_(r.dist_ewald * r.dist_ewald) {
         }
 
         template<typename T>
@@ -246,6 +246,7 @@ namespace {
     void scale(const ScaleMergeOptions &opt,
                std::vector<double> &g,
                std::vector<double> &mosaicity,
+               std::vector<double> &R_sq,
                const std::vector<uint8_t> &image_slot_used,
                bool rotation_crystallography,
                size_t nhkl,
@@ -293,7 +294,7 @@ namespace {
                         problem.AddResidualBlock(cost,
                                                  nullptr,
                                                  &g[o.img_id],
-                                                 &mosaicity[o.img_id],
+                                                 &R_sq[o.img_id],
                                                  &Itrue[o.hkl_slot]);
                     }
                 case ScaleMergeOptions::PartialityModel::Unity: {
@@ -350,10 +351,10 @@ namespace {
         }
 
         if (opt.partiality_model == ScaleMergeOptions::PartialityModel::Still) {
-            for (int i = 0; i < mosaicity.size(); ++i) {
+            for (int i = 0; i < R_sq.size(); ++i) {
                 if (image_slot_used[i]) {
-                    problem.SetParameterLowerBound(&mosaicity[i], 0, 1e-9);
-                    problem.SetParameterUpperBound(&mosaicity[i], 0, 1.0);
+                    problem.SetParameterLowerBound(&R_sq[i], 0, 1e-9);
+                    problem.SetParameterUpperBound(&R_sq[i], 0, 1.0);
                 }
             }
         }
@@ -590,6 +591,7 @@ namespace {
     void calc_obs(const ScaleMergeOptions &opt,
                std::vector<double> &g,
                std::vector<double> &mosaicity,
+               std::vector<double> &R_sq,
                const std::vector<ObsRef> &obs,
                std::vector<CorrectedObs> &corr_obs) {
 
@@ -616,7 +618,7 @@ namespace {
                 }
                     break;
                 case ScaleMergeOptions::PartialityModel::Still:
-                    partiality = std::exp(-r.dist_ewald * r.dist_ewald / mosaicity[o.img_id]);
+                    partiality = std::exp(-r.dist_ewald * r.dist_ewald / R_sq[o.img_id]);
                     break;
                 case ScaleMergeOptions::PartialityModel::Unity:
                     break;
@@ -721,16 +723,19 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
 
     std::vector<double> g(n_image_slots, 1.0);
     std::vector<double> mosaicity(n_image_slots, opt.mosaicity_init_deg);
+    std::vector<double> R_sq(n_image_slots, 0.001 * 0.001);
+
     for (int i = 0; i < n_image_slots; i++) {
         if (!image_slot_used[i]) {
             mosaicity[i] = NAN;
             g[i] = NAN;
+            R_sq[i] = NAN;
         } else if (opt.mosaicity_init_deg_vec.size() > i && std::isfinite(opt.mosaicity_init_deg_vec[i])) {
             mosaicity[i] = opt.mosaicity_init_deg_vec[i];
         }
     }
 
-    scale(opt, g, mosaicity, image_slot_used, rotation_crystallography, nhkl, obs);
+    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs);
 
     ScaleMergeResult out;
 
@@ -775,7 +780,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         }
     }
 
-    calc_obs(opt, g, mosaicity, obs, corr_obs);
+    calc_obs(opt, g, mosaicity, R_sq, obs, corr_obs);
     merge(nhkl, out, corr_obs);
     stats(opt, nhkl, out, corr_obs);