diff --git a/cqlib/benchmark/irb.py b/cqlib/benchmark/irb.py
index 36d86b938044b0f71f8505b6fb9aad3e404c9fc0..b5aaa33c3b29d144512f573540c5f09afdc4ca98 100644
--- a/cqlib/benchmark/irb.py
+++ b/cqlib/benchmark/irb.py
@@ -251,25 +251,43 @@ def _fit_decay(
     p_init: float,
     sigma_floor: float,
 ) -> dict[str, float | None]:
-    """Fit y(m)=A*p^m+B for one RB branch (ref or int)."""
+    """Fit y(m)=A*p^m+B for one RB branch (ref or int).
+
+    The fit is intentionally *unweighted*. Weighting by the per-length SEM
+    (Standard Error of the Mean), as in ``curve_fit(sigma=sem,
+    absolute_sigma=True)``, over-weights the short-sequence points (small SEM)
+    and suppresses the long-sequence tail (whose SEM grows with m). That tail
+    is precisely what constrains the asymptote B, so SEM-weighting leaves B
+    under-determined and a single-start fit lands in inconsistent basins,
+    destabilizing p (and hence the gate fidelity p_int/p_ref). We fit
+    unweighted with a small multi-start over the offset B and keep the lowest
+    sum-of-squared-errors solution. ``sem``/``sigma_floor`` are retained for
+    API compatibility.
+    """
 
     def model(m: np.ndarray, a: float, p: float, b: float) -> np.ndarray:
         return a * (p**m) + b
 
-    sigma = np.where(sem > 1e-6, sem, float(sigma_floor))
     a0 = float(max(mean[0] - mean[-1], 1e-3))
-    b0 = float(np.clip(mean[-1], 0.0, 1.0))
-
-    params, cov = curve_fit(
-        model,
-        m_data,
-        mean,
-        p0=[a0, p_init, b0],
-        sigma=sigma,
-        absolute_sigma=True,
-        bounds=([-1.0, 0.0, 0.0], [1.0, 1.0, 1.0]),
-        maxfev=20000,
-    )
+    best: tuple[float, np.ndarray, np.ndarray] | None = None
+    for b0 in (0.0, 0.25, 0.5):
+        try:
+            params, cov = curve_fit(
+                model,
+                m_data,
+                mean,
+                p0=[a0, p_init, b0],
+                bounds=([-1.0, 0.0, 0.0], [1.0, 1.0, 1.0]),
+                maxfev=20000,
+            )
+        except RuntimeError:
+            continue
+        sse = float(np.sum((model(m_data, *params) - mean) ** 2))
+        if best is None or sse < best[0]:
+            best = (sse, params, cov)
+    if best is None:
+        raise RuntimeError("RB decay fit failed to converge from any start")
+    _, params, cov = best
     a, p, b = [float(v) for v in params]
     p_std = float(np.sqrt(max(cov[1, 1], 0.0))) if cov.shape == (3, 3) else None
     y_hat = model(m_data, a, p, b)