diff --git a/benchmarks/zdc_lambda/analysis/lambda_plots.py b/benchmarks/zdc_lambda/analysis/lambda_plots.py
index bf855d9..931a44b 100644
--- a/benchmarks/zdc_lambda/analysis/lambda_plots.py
+++ b/benchmarks/zdc_lambda/analysis/lambda_plots.py
@@ -191,7 +191,7 @@ def gauss(x, A,mu, sigma):
 fnc=gauss
 p0=[100, 0, 0.05]
 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}_{\\Lambda}-\\theta^{*\\rm truth}_{\\Lambda}$ [mrad]")
@@ -214,7 +214,7 @@ def gauss(x, A,mu, sigma):
     #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(coeff[2])
         dsigmas.append(np.sqrt(var_matrix[2][2]))
         xvals.append(p)
@@ -259,7 +259,7 @@ def gauss(x, A,mu, sigma):
 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')
 print(coeff[2], np.sqrt(var_matrix[2][2]))
@@ -284,7 +284,7 @@ def gauss(x, A,mu, sigma):
     #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(coeff[2])
         dsigmas.append(np.sqrt(var_matrix[2][2]))
         xvals.append(p)
@@ -327,7 +327,7 @@ def gauss(x, A,mu, sigma):
 fnc=gauss
 p0=[100, lambda_mass, 0.04]
 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]))
@@ -350,7 +350,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, lambda_mass, 0.05]
     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]))
diff --git a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
index 94a9310..a2cf533 100644
--- a/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
+++ b/benchmarks/zdc_photon/analysis/zdc_photon_plots.py
@@ -57,7 +57,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, p, 10]
     #print(list(y), list(x))
     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)
     if p==100:
         xx=np.linspace(p*0.75,p*1.25, 100)
         plt.plot(xx, fnc(xx,*coeff))
@@ -78,7 +78,7 @@ def gauss(x, A,mu, sigma):
 fnc=lambda E,a: a/np.sqrt(E)
 #pvals, resvals, dresvals
 coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
-                                 sigma=dresvals)
+                                 sigma=dresvals, maxfev=10000)
 
 xx=np.linspace(15, 275, 100)
 plt.plot(xx, fnc(xx, *coeff), label=f'fit:  $\\frac{{{coeff[0]*100:.0f}\\%}}{{\\sqrt{{E}}}}$')
@@ -129,7 +129,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, 0, 0.1]
     #print(list(y), list(x))
     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)
     if p==100:
         xx=np.linspace(-0.5,0.5, 100)
         plt.plot(xx, fnc(xx,*coeff))
@@ -143,7 +143,7 @@ def gauss(x, A,mu, sigma):
 fnc=lambda E,a, b: np.hypot(a/np.sqrt(E), b)
 #pvals, resvals, dresvals
 coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,.1),
-                                 sigma=dresvals)
+                                 sigma=dresvals, maxfev=10000)
 
 xx=np.linspace(15, 275, 100)
 
diff --git a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
index 1b54430..6d0cddd 100644
--- a/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
+++ b/benchmarks/zdc_pi0/analysis/zdc_pi0_plots.py
@@ -57,7 +57,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, p, 10]
     #print(list(y), list(x))
     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)
     if p==100:
         xx=np.linspace(p*0.5,p*1.5, 100)
         plt.plot(xx, fnc(xx,*coeff))
@@ -76,7 +76,7 @@ def gauss(x, A,mu, sigma):
 fnc=lambda E,a: a/np.sqrt(E)
 #pvals, resvals, dresvals
 coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
-                                 sigma=dresvals)
+                                 sigma=dresvals, maxfev=10000)
 xx=np.linspace(55, 200, 100)
 plt.plot(xx, fnc(xx, *coeff), label=f'fit:  $\\frac{{{coeff[0]:.2f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -133,7 +133,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, 0, 0.1]
     #print(list(y), list(x))
     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)
     if p==100:
         xx=np.linspace(-0.5,0.5, 100)
         plt.plot(xx, fnc(xx,*coeff))
@@ -148,7 +148,7 @@ def gauss(x, A,mu, sigma):
 fnc=lambda E,a: a/np.sqrt(E)
 #pvals, resvals, dresvals
 coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,),
-                                 sigma=dresvals)
+                                 sigma=dresvals, maxfev=10000)
 
 xx=np.linspace(55, 200, 100)
 
@@ -201,7 +201,7 @@ def gauss(x, A,mu, sigma):
     p0=[100, .135, 0.2]
     #print(list(y), list(x))
     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)
     if p==100:
         xx=np.linspace(0,0.2)
         plt.plot(xx, fnc(xx,*coeff))
@@ -218,7 +218,7 @@ def gauss(x, A,mu, sigma):
 fnc=lambda E,a,b: a+b*E
 #pvals, resvals, dresvals
 coeff, var_matrix = curve_fit(fnc, pvals, resvals, p0=(1,1),
-                                 sigma=dresvals)
+                                 sigma=dresvals, maxfev=10000)
 xx=np.linspace(55, 200, 100)
 #plt.plot(xx, fnc(xx, *coeff), label=f'fit:  ${coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times E$ MeV')
 plt.plot(xx, fnc(xx, *coeff), label=f'fit:  $({coeff[0]*1000:.1f}+{coeff[1]*1000:.4f}\\times [E\,in\,GeV])$ MeV')