diff --git a/benchmarks/femc_electron/analysis/femc_electron_plots.py b/benchmarks/femc_electron/analysis/femc_electron_plots.py
index cc34633dc30887057bda2e8331f1ed509c6f33b5..8979e0c79efc5b4fad6cc6c5bd6f1268eb3fdbf0 100644
--- a/benchmarks/femc_electron/analysis/femc_electron_plots.py
+++ b/benchmarks/femc_electron/analysis/femc_electron_plots.py
@@ -132,7 +132,7 @@ for p in arrays_sim:
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -151,7 +151,7 @@ plt.sca(axs[1])
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 xx=np.linspace(7, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -190,7 +190,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)
diff --git a/benchmarks/femc_photon/analysis/femc_photon_plots.py b/benchmarks/femc_photon/analysis/femc_photon_plots.py
index 8170f4a1d0ecd006a5f511bab246cdd0124f9cdd..176eb663801e58edcef015d9e2a3c10ef2eef619 100644
--- a/benchmarks/femc_photon/analysis/femc_photon_plots.py
+++ b/benchmarks/femc_photon/analysis/femc_photon_plots.py
@@ -131,7 +131,7 @@ for p in arrays_sim:
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -150,7 +150,7 @@ plt.sca(axs[1])
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 xx=np.linspace(7, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -188,7 +188,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
         p0=(100, p, 3)
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)
diff --git a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
index cdf9d2319c943721bb8fdf281dd34d7bda7fd069..0bdf26e340d50252f5d2c83af0910e14ed512f03 100644
--- a/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
+++ b/benchmarks/femc_pi0/analysis/femc_pi0_plots.py
@@ -167,7 +167,7 @@ for p in arrays_sim:
     sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
     p0=(100, p, 3)
 
-    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+    coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
     #res=np.abs(coeff[2]/coeff[1])
     if p==50:
         xx=np.linspace(15*p/20,22*p/20, 100)
@@ -186,7 +186,7 @@ plt.sca(axs[1])
 plt.errorbar(pvals, 100*np.array(res), 100*np.array(dres), ls='', marker='o')
 fnc = lambda E, a, b: np.hypot(a,b/np.sqrt(E))
 p0=(.05, .12)
-coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0,sigma=dres)
+coeff, var_matrix = curve_fit(fnc, pvals, res, p0=p0, sigma=dres, maxfev=10000)
 xx=np.linspace(15, 85, 100)
 plt.plot(xx, 100*fnc(xx,*coeff), label=f'fit:{100*coeff[0]:.1f}%$\\oplus\\frac{{{100*coeff[1]:.0f}\\%}}{{\\sqrt{{E}}}}$')
 plt.legend()
@@ -225,7 +225,7 @@ for eta_min, eta_max in zip(partitions[:-1], partitions[1:]):
         p0=(100, p, 3)
 
         try:
-            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0,sigma=list(sigma))
+            coeff, var_matrix = curve_fit(fnc, list(bcs[slc]), list(y[slc]), p0=p0, sigma=list(sigma), maxfev=10000)
             if abs(coeff[1])>100 or np.sqrt(var_matrix[1][1])>100:
                 continue
             pvals.append(p)