diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index b2e1fafa..80d1b6d0 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -1253,8 +1253,10 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
     // intensity raises its systematic term `b`, so 1/b understates that limit. Read it instead directly
     // from the strong equivalents: for each well-measured reflection group the counting-subtracted
     // fractional scatter of its symmetry mates estimates the systematic term, and the robust median over
-    // strong groups is the asymptote. Falls back to the whole-range `b` when too few strong groups exist
-    // (weak / low-multiplicity data). Host-side over the merged fulls, so CPU and GPU agree.
+    // strong groups is the asymptote. The I/sigma threshold that selects "strong" is relaxed on weak /
+    // radiation-damaged data that has few strong reflections (rather than fall back to the higher
+    // whole-range b); it falls back only when even the relaxed set is too small. Host-side over the
+    // merged fulls, so CPU and GPU agree.
     double error_model_b_asymptotic = error_model_b;
     if (error_model_active) {
         struct GroupScatter { double sum = 0, sum_sq = 0, sum_var = 0; int n = 0; };
@@ -1266,19 +1268,26 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool
             auto &g = gs[o.group];
             g.sum += I_corr; g.sum_sq += I_corr * I_corr; g.sum_var += sigma_corr * sigma_corr; ++g.n;
         }
-        std::vector<double> b2_group;
+        // Per-group counting-subtracted fractional systematic variance, paired with the group's I/sigma.
+        std::vector<std::pair<double, double>> group_scatter;   // (systematic b^2, I/sigma)
         for (const auto &g : gs) {
             if (g.n < 5) continue;
             const double mean = g.sum / g.n;
-            const double variance = (g.sum_sq - g.sum * g.sum / g.n) / (g.n - 1);   // scatter of the mates
-            const double counting = g.sum_var / g.n;                                // mean counting variance
-            if (mean <= 0.0 || mean * mean < 400.0 * counting) continue;            // strong groups (I/sig > 20)
-            b2_group.push_back(std::max(variance - counting, 0.0) / (mean * mean));
-        }
-        if (b2_group.size() >= 100) {
-            const double b_asy = std::sqrt(median_of(b2_group));
-            if (b_asy > 0.0) error_model_b_asymptotic = b_asy;
+            const double counting = g.sum_var / g.n;
+            if (mean <= 0.0 || counting <= 0.0) continue;
+            const double variance = (g.sum_sq - g.sum * g.sum / g.n) / (g.n - 1);
+            group_scatter.push_back({std::max(variance - counting, 0.0) / (mean * mean),
+                                     mean / std::sqrt(counting)});
         }
+        auto asymptote_above = [&](double snr_min, size_t min_groups) -> double {
+            std::vector<double> b2;
+            for (const auto &[b2_value, snr] : group_scatter)
+                if (snr >= snr_min) b2.push_back(b2_value);
+            return b2.size() >= min_groups ? std::sqrt(median_of(b2)) : 0.0;
+        };
+        double b_asy = asymptote_above(20.0, 100);   // tight threshold on data that supports it
+        if (b_asy <= 0.0) b_asy = asymptote_above(10.0, 50);   // relaxed for weak / damaged data
+        if (b_asy > 0.0) error_model_b_asymptotic = b_asy;
     }
     if (error_model_active)
         logger.Info("Error model: a={:.3f} b={:.3f} ISa={:.1f} chi2={:.2f}", error_model_a, error_model_b,