diff --git a/benchmarks/insert_muon/analysis/muon_plots.py b/benchmarks/insert_muon/analysis/muon_plots.py
index 5c81b52..ac2aa10 100644
--- a/benchmarks/insert_muon/analysis/muon_plots.py
+++ b/benchmarks/insert_muon/analysis/muon_plots.py
@@ -55,7 +55,7 @@ def Landau(x, normalization,location,stdev):
     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 b9d2002..424295c 100644
--- a/benchmarks/insert_neutron/analysis/neutron_plots.py
+++ b/benchmarks/insert_neutron/analysis/neutron_plots.py
@@ -79,7 +79,7 @@ def gauss(x, A,mu, sigma):
 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 @@ 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(np.abs(coeff[2]))
             dsigmas.append(np.sqrt(var_matrix[2][2]))
             xvals.append(p)
@@ -151,7 +151,7 @@ def gauss(x, A,mu, sigma):
         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 @@ def gauss(x, A,mu, sigma):
         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: