How can I use Gaussian Mixture to fit 2 gaussians in my data?

Question

I'm trying to fit 2 gaussians in my red/blue population data. But it does not fit. What am I doing wrong?

The code is:

mag1 = 'z'
mag2 = 'y'
mask_rich = (wazp_membership['NGALS']>25)&(wazp_membership['NGALS']<35)
mask_z = (wazp_membership['ZP_CL']>0.6)&(wazp_membership['ZP_CL']<0.7)
mask = (mask_rich)&(mask_z)

wazp_membership_color = wazp_membership[f'mag_{mag1}'] - wazp_membership[f'mag_{mag2}']
color = np.array(wazp_membership_color[mask])
pmem = np.array(wazp_membership['PMEM'][mask])

plt.figure(figsize=[8,5])
y,x,_=plt.hist(color, bins=30, weights=pmem, fc="moccasin", ec='sandybrown', label='color distribution')
x=(x[1:]+x[:-1])/2
#plt.title(f'Color histogram ({round(z1,1)}<$z$<{round(z2,1)} and {ngals_min}<$Ngals$<{ngals_max})', fontsize=12)

x = x.reshape(-1, 1)

gmm = GaussianMixture(n_components=2, covariance_type="full", tol=0.001)


mean = gmm.fit(x).means_  
covs  = gmm.fit(x).covariances_
weights = gmm.fit(x).weights_

x_axis = np.linspace(x.min(), x.max(), 500)
y_axis0 = norm.pdf(x_axis, float(mean[0][0]), np.sqrt(float(covs[0][0][0])))*weights[0] # 1st gaussian
y_axis1 = norm.pdf(x_axis, float(mean[1][0]), np.sqrt(float(covs[1][0][0])))*weights[1] # 2nd gaussian

plt.plot(x_axis, y_axis0, lw=3, c='C0')
plt.plot(x_axis, y_axis1, lw=3, c='C1')
plt.plot(x_axis, y_axis0+y_axis1, lw=3, c='C2', ls='dashed')

And I get this: enter image description here