def METHOD_NAME(p, X, Y, Z):
    """Pair of 3D Gaussian model functions with linear background
     - parameter vector [A, x0, y0, z0, x1, y1, z1, wxy, wz, background]
     Note: Assumes sames sigma for both Gaussian (dumbell). """
    A, x0, y0, z0, x1, y1, z1, sxy, sz, b = p
    r0 = A*np.exp(-((X[:,None, None]-x0)**2 + (Y[None,:,None]-y0)**2)/(2*sxy**2) - ((Z[None,None, :]-z0)**2)/(2*sz**2))
    r1 = A*np.exp(-((X[:,None, None]-x1)**2 + (Y[None,:,None]-y1)**2)/(2*sxy**2) - ((Z[None,None, :]-z1)**2)/(2*sz**2))
    r = r0+r1+b
    return r