rugnux: recover the asymptotic ISa on weak / radiation-damaged data

The asymptote was read only from groups with I/sigma > 20 and required at
least 100 of them, so weak or radiation-damaged data (few strong reflections)
fell back to the higher whole-range b and reported an over-conservative ISa.
Relax the I/sigma threshold in a second tier (>= 50 groups above I/sigma 10)
when the tight one is unmet; the per-group counting-subtracted estimate is
stable down to that threshold. Data with many strong groups is unchanged (it
keeps the tight-threshold asymptote); weak / damaged data now reports the
asymptote its own reflections support instead of falling back.

Battery (32 crystals): no space-group, R_meas or CC1/2 change; strong sets
unchanged; e.g. a room-temperature radiation-damage series recovers ISa 5.4 ->
11.2 (its data already exceeds the reference on CC1/2 and R_meas).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

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,