diff --git a/src/aspire/abinitio/commonline_sync3n.py b/src/aspire/abinitio/commonline_sync3n.py
index db1dc69c12..6a7e5390d3 100644
--- a/src/aspire/abinitio/commonline_sync3n.py
+++ b/src/aspire/abinitio/commonline_sync3n.py
@@ -714,13 +714,26 @@ def _triangle_scores(
             * (a + 1)
         )
         # normalization of 2nd component: B = P*N_delta/sum(f), where f is the component formula
-        B0 = P ** (self.n_img * (self.n_img - 1) * (self.n_img - 2) / 2) / np.sum(
-            ((1 - hist_x) ** b) * np.exp(-b / (1 - x0) * (1 - hist_x))
-        )
-        start_values = np.array([B0, P, b, x0], dtype=np.float64)
+        # B0 = (
+        #     P
+        #     * (self.n_img * (self.n_img - 1) * (self.n_img - 2) / 2)
+        #     / np.sum(((1 - hist_x) ** b) * np.exp(-b / (1 - x0) * (1 - hist_x)))
+        # )
+        # P must be in lower and upper bounds or `curve_fit` will error
+        # This was not the case for MATLAB...
+        # P0 = np.clip(P, Pmin**3, Pmax**3)
+        # Note, MATLAB suggests the following, but I feel it is a bug.
+        # Will discuss with Yoel about the original code's intent.
+        #     np.array([B0, P0, b, x0], dtype=np.float64)
+        start_values = None
         lower_bounds = np.array([0, Pmin**3, 2, 0], dtype=np.float64)
         upper_bounds = np.array([np.inf, Pmax**3, np.inf, 1], dtype=np.float64)
 
+        with np.printoptions(precision=2):
+            logger.info(f"curve_fit lower_bounds:{lower_bounds}")
+            logger.info(f"curve_fit start_values:{start_values}")
+            logger.info(f"curve_fit upper_bounds:{upper_bounds}")
+
         # Fit distribution
         def fun(x, B, P, b, x0, A=A, a=a):
             """Function to fit. x is data vector."""
@@ -742,6 +755,8 @@ def fun(x, B, P, b, x0, A=A, a=a):
         P = P ** (1 / 3)
         sigma = (1 - x0) / peak2sigma
 
+        logger.info(f"Estimated CL Errors P,STD:\t{100*P:.2f}%\t{sigma:.2f}")
+
         # Initialize probability computations
         # Local histograms analysis
         A = a + 1  # distribution 1st component normalization factor
@@ -768,6 +783,10 @@ def fun(x, B, P, b, x0, A=A, a=a):
             )
             Pij = np.nan_to_num(Pij)
 
+        logger.info(
+            f"Common lines probabilities to be indicative Pij={100*np.mean(Pij):.2f}%"
+        )
+
         return P, sigma, Pij, scores_hist
 
     ###########################################