def METHOD_NAME(hbar):
    """Test that probabilities are correctly updates for a four modes system under loss"""
    # All this block is to generate the correct covariance matrix.
    # It correnponds to num_modes=4 modes that undergo two mode squeezing between modes i and i + (num_modes / 2).
    # Then they undergo displacement.
    # The signal and idlers see and interferometer with unitary matrix u2x2.
    # And then they see loss by amount etas[i].
    num_modes = 4
    theta = 0.45
    phi = 0.7
    u2x2 = np.array(
        [
            [np.cos(theta / 2), np.exp(1j * phi) * np.sin(theta / 2)],
            [-np.exp(-1j * phi) * np.sin(theta / 2), np.cos(theta / 2)],
        ]
    )
    u4x4 = block_diag(u2x2, u2x2)
    cov = np.identity(2 * num_modes) * hbar / 2
    means = 0.5 * np.random.rand(2 * num_modes) * np.sqrt(hbar / 2)
    rs = [0.1, 0.9]
    n_half = num_modes // 2
    for i, r_val in enumerate(rs):
        Sexpanded = expand(two_mode_squeezing(r_val, 0.0), [i, n_half + i], num_modes)
        cov = Sexpanded @ cov @ (Sexpanded.T)
    Su = expand(interferometer(u4x4), range(num_modes), num_modes)
    cov = Su @ cov @ (Su.T)
    cov_lossless = np.copy(cov)
    means_lossless = np.copy(means)
    etas = [0.9, 0.7, 0.9, 0.1]
    for i, eta in enumerate(etas):
        means, cov = loss(means, cov, eta, i, hbar=hbar)
    cutoff = 3
    probs_lossless = probabilities(means_lossless, cov_lossless, 4 * cutoff, hbar=hbar)
    probs = probabilities(means, cov, cutoff, hbar=hbar)
    probs_updated = update_probabilities_with_loss(etas, probs_lossless)
    assert np.allclose(probs, probs_updated[:cutoff, :cutoff, :cutoff, :cutoff], atol=1e-6)