res = c1 * np.exp(-(x - mu1) ** 2.0 / (2.0 * sigma1 ** 2.0)) + c2 * np.exp(
    -(x - mu2) ** 2.0 / (2.0 * sigma2 ** 2.0))
return res

def double_gaussian_fit(params):
    fit = double_gaussian(x, params)
return fit - y_proc
y_proc = np.copy(y)
y_proc[y_proc < 5] = 0.0
fit = leastsq(double_gaussian_fit, [13.0, -13.0, 1.0, 60.0, 3.0, 1.0])
plot(x, y, c='b')
plot(x, double_gaussian(x, fit[0]), c='r')
How to fit a double Gaussian distribution in Python?