diff --git a/benchmarks/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py
index 5c81b52dccb918484991e88e35b308623b68601e..ac2aa100971c1595394b0d32ce79ccf2c4edfd83 100644
--- a/benchmarks/insert_muon/analysis/muon_plots.py
+++ b/benchmarks/insert_muon/analysis/muon_plots.py
@@ -55,7 +55,7 @@ for p in 50,:
     p0=[100, .5, .05]
     #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)
     print(coeff)
     xx=np.linspace(0,.7, 100)
     MIP=coeff[1]/1000
diff --git a/benchmarks/insert_neutron/analysis/neutron_plots.py b/benchmarks/insert_neutron/analysis/neutron_plots.py
index b9d2002a38372b9799c7bb5d90f9904c9d22e31d..424295cf211a5c0657a9df4433b7bae758e008be 100644
--- a/benchmarks/insert_neutron/analysis/neutron_plots.py
+++ b/benchmarks/insert_neutron/analysis/neutron_plots.py
@@ -79,7 +79,7 @@ slc=abs(bc)<3
 fnc=gauss
 sigma=np.sqrt(y[slc])+(y[slc]==0)
 p0=(100, 0, 5)
-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)
 xx=np.linspace(-5,5,100)
 plt.plot(xx,fnc(xx,*coeff))
 # except:
@@ -104,7 +104,7 @@ for eta_min, eta_max in zip(r[:-1],r[1:]):
         #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)
@@ -151,7 +151,7 @@ for p in 20, 30,40,50,60,70, 80:
         sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
         p0=(100, np.mean(r), np.std(r))
         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)
             res=np.abs(coeff[2]/coeff[1])
             
             if res<best_res:
@@ -234,7 +234,7 @@ for eta_min, eta_max in zip(partitions[:-1],partitions[1:]):
         sigma=np.sqrt(y[slc])+0.5*(y[slc]==0)
         p0=(100, np.mean(r), np.std(r))
         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)
             res=np.abs(coeff[2]/coeff[1])
 
             if res<best_res: