gammapy: Covariance matrix is not properly set on Datasets

Gammapy version from at least v0.17 to 1.0rc1

Bug description The covariance matrix is not properly set on Datasets when Models contain more than one model. Only the sub-covariances of individual models are properly set. No positive covariance terms between independent model components are set. Thanks @AtreyeeS for noting this.

Expected behavior Correlations between independent model parameters should be non-zero.

To Reproduce

import numpy as np
from gammapy.datasets import Datasets
from gammapy.modeling import Covariance

ds = Datasets.read("$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml",
              filename_models= "$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_models.yaml")

factor_matrix = np.ones((9,9))
ds.models.covariance = Covariance.from_factor_matrix(
        parameters=ds.parameters, matrix=factor_matrix
    )
ds.models.covariance.plot_correlation()

yields:

while:

ds[0].models.covariance = Covariance.from_factor_matrix(
        parameters=ds[0].models.parameters, matrix=factor_matrix
    )
ds[0].models.covariance.plot_correlation()

yields:

But ds.models.covariance.plot_correlation() still plots the incorrect matrix.