diff --git a/common/BraggIntegrationSettings.h b/common/BraggIntegrationSettings.h
index 073a5d30..b8292916 100644
--- a/common/BraggIntegrationSettings.h
+++ b/common/BraggIntegrationSettings.h
@@ -16,7 +16,7 @@ class BraggIntegrationSettings {
     IntegratorMode integrator_mode = IntegratorMode::ProfileGaussian;
     float r_1 = 4;
     float r_2 = 6;
-    float r_3 = 8;
+    float r_3 = 10;
     float d_min_limit_A = 1.0;
     std::optional<float> fixed_profile_radius;
     float minimum_sigma_in_regards_to_i = 0.02;
diff --git a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp
index 29accbee..39a347bb 100644
--- a/image_analysis/bragg_integration/ProfileIntegrate2D.cpp
+++ b/image_analysis/bragg_integration/ProfileIntegrate2D.cpp
@@ -17,6 +17,28 @@ constexpr int N_SHELL = 6;
 constexpr double STRONG_I_OVER_SIGMA = 5.0;
 constexpr int MIN_STRONG_PER_SHELL = 30;
 
+// Mark the r2 signal disk of every predicted reflection. A neighbour whose disk falls inside this
+// reflection's r2..r3 background ring would otherwise bias the background mean high and over-subtract;
+// excluding the marked pixels from the ring removes that contamination (same as BraggIntegrate2D).
+std::vector<uint8_t> BuildReflectionMask(const std::vector<Reflection> &predicted, size_t npredicted,
+                                         size_t xpixel, size_t ypixel, float r2) {
+    std::vector<uint8_t> mask(xpixel * ypixel, 0);
+    const float r2_sq = r2 * r2;
+    for (size_t i = 0; i < npredicted; ++i) {
+        const auto &r = predicted[i];
+        const int64_t x0 = std::max<int64_t>(0, std::floor(r.predicted_x - r2 - 1.0f));
+        const int64_t x1 = std::min<int64_t>(xpixel - 1, std::ceil(r.predicted_x + r2 + 1.0f));
+        const int64_t y0 = std::max<int64_t>(0, std::floor(r.predicted_y - r2 - 1.0f));
+        const int64_t y1 = std::min<int64_t>(ypixel - 1, std::ceil(r.predicted_y + r2 + 1.0f));
+        for (int64_t y = y0; y <= y1; ++y)
+            for (int64_t x = x0; x <= x1; ++x) {
+                const double d2 = (x - r.predicted_x) * (x - r.predicted_x) + (y - r.predicted_y) * (y - r.predicted_y);
+                if (d2 < r2_sq) mask[y * xpixel + x] = 1;
+            }
+    }
+    return mask;
+}
+
 // Rough box-sum of one reflection: total intensity over the inner disk (r1) minus the mean of the
 // background ring (r2..r3), like BraggIntegrate2D. Used to seed the fit, pick strong spots for
 // profile learning, and provide the per-reflection background. ok=false when the inner disk is
@@ -29,7 +51,8 @@ struct Rough {
 
 template<class T>
 Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int64_t saturation,
-             const Reflection &r, float r1_sq, float r2_sq, float r3_sq, float r3, float min_sigma_ratio,
+             const Reflection &r, const std::vector<uint8_t> &refl_mask,
+             float r1_sq, float r2_sq, float r3_sq, float r3, float min_sigma_ratio,
              bool apply_bkg_clip) {
     Rough out;
     const int64_t x0 = std::max<int64_t>(0, std::floor(r.predicted_x - r3 - 1.0));
@@ -52,6 +75,7 @@ Rough BoxSum(const T *image, size_t xpixel, size_t ypixel, int64_t special, int6
                 I_sum += px;
                 ++n_inner_valid;
             } else if (d2 >= r2_sq && d2 < r3_sq) {
+                if (refl_mask[y * xpixel + x]) continue;
                 if (px == special || px == special + 1 || px == saturation || px == saturation - 1) continue;
                 bkg_sum += static_cast<double>(px);
                 ++n_bkg;
@@ -110,10 +134,11 @@ std::vector<Reflection> ProfileIntegrateInternal(const DiffractionExperiment &ex
     const size_t xpixel = image.GetWidth(), ypixel = image.GetHeight();
 
     // --- Pass A: box-sum every reflection (rough I, background, centre, strong flag). ---
+    const auto refl_mask = BuildReflectionMask(predicted, npredicted, xpixel, ypixel, r2);
     std::vector<Rough> rough(npredicted);
     double inv_d2_min = std::numeric_limits<double>::max(), inv_d2_max = 0.0;
     for (size_t i = 0; i < npredicted; ++i) {
-        rough[i] = BoxSum<T>(ptr, xpixel, ypixel, special, saturation, predicted[i],
+        rough[i] = BoxSum<T>(ptr, xpixel, ypixel, special, saturation, predicted[i], refl_mask,
                              r1_sq, r2_sq, r3_sq, r3, min_sigma_ratio, apply_bkg_clip);
         if (rough[i].ok && predicted[i].d > 0.0f) {
             const double inv_d2 = 1.0 / (predicted[i].d * predicted[i].d);