diff --git a/benchmarks/zdc_sigma/analysis/sigma_plots.py b/benchmarks/zdc_sigma/analysis/sigma_plots.py
index fc422ed6edc75261c2d55f30a41c7efa241580d8..7ded93dae39c069fa9974d531e756f4883dc5f61 100644
--- a/benchmarks/zdc_sigma/analysis/sigma_plots.py
+++ b/benchmarks/zdc_sigma/analysis/sigma_plots.py
@@ -236,7 +236,7 @@ slc=abs(bc)<0.6
 fnc=gauss
 p0=[100, 0, 0.5]
 coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
-                                 sigma=np.sqrt(y[slc])+(y[slc]==0))
+                                 sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
 x=np.linspace(-1, 1)
 plt.plot(x, gauss(x, *coeff), color='tab:orange')
 plt.xlabel("$\\theta^{*\\rm recon}_{\\Sigma}-\\theta^{*\\rm truth}_{\\Sigma}$ [mrad]")
@@ -259,7 +259,7 @@ for p in momenta:
     #print(bc[slc],y[slc])
     sigma=np.sqrt(y[slc])+(y[slc]==0)
     try:
-        coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+        coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
         sigmas.append(np.abs(coeff[2]))
         dsigmas.append(np.sqrt(var_matrix[2][2]))
         xvals.append(p)
@@ -307,7 +307,7 @@ slc=abs(bc)<5
 fnc=gauss
 p0=[100, 0, 1]
 coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
-                                 sigma=np.sqrt(y[slc])+(y[slc]==0))
+                                 sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
 x=np.linspace(-5, 5)
 plt.plot(x, gauss(x, *coeff), color='tab:orange')
 plt.xlabel("$z^{*\\rm recon}_{\\rm vtx}-z^{*\\rm truth}_{\\rm vtx}$ [m]")
@@ -331,7 +331,7 @@ for p in momenta:
     #print(bc[slc],y[slc])
     sigma=np.sqrt(y[slc])+(y[slc]==0)
     try:
-        coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+        coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
         sigmas.append(abs(coeff[2]))
         dsigmas.append(np.sqrt(var_matrix[2][2]))
         xvals.append(p)
@@ -373,7 +373,7 @@ slc=abs(bc-lambda_mass)<0.05
 fnc=gauss
 p0=[100, lambda_mass, 0.03]
 coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
-                                 sigma=np.sqrt(y[slc])+(y[slc]==0))
+                                 sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
 x=np.linspace(0.8, 1.3, 200)
 plt.plot(x, gauss(x, *coeff), color='tab:orange')
 print(coeff[2], np.sqrt(var_matrix[2][2]))
@@ -396,7 +396,7 @@ for p in momenta:
     p0=[100, lambda_mass, 0.03]
     try:
         coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
-                                       sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+                                       sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
         x=np.linspace(0.8, 1.3, 200)
         sigmas.append(coeff[2])
         dsigmas.append(np.sqrt(var_matrix[2][2]))
@@ -437,7 +437,7 @@ slc=abs(bc-sigma_mass)<0.02
 fnc=gauss
 p0=[100, sigma_mass, 0.03]
 coeff, var_matrix = curve_fit(fnc, bc[slc], y[slc], p0=p0,
-                                 sigma=np.sqrt(y[slc])+(y[slc]==0))
+                                 sigma=np.sqrt(y[slc])+(y[slc]==0), maxfev=10000)
 x=np.linspace(0.8, 1.3, 200)
 plt.plot(x, gauss(x, *coeff), color='tab:orange')
 print(coeff[2], np.sqrt(var_matrix[2][2]))
@@ -460,7 +460,7 @@ for p in momenta:
     p0=[100, sigma_mass, 0.03]
     try:
         coeff, var_matrix = curve_fit(fnc, list(bc[slc]), list(y[slc]), p0=p0,
-                                       sigma=list(np.sqrt(y[slc])+(y[slc]==0)))
+                                       sigma=list(np.sqrt(y[slc])+(y[slc]==0)), maxfev=10000)
         sigmas.append(abs(coeff[2]))
         dsigmas.append(np.sqrt(var_matrix[2][2]))
         xvals.append(p)