def METHOD_NAME(self, datasets):
    """Run flux profile estimation
    Parameters
    ----------
    datasets : list of `~gammapy.datasets.MapDataset`
        Map datasets.
    Returns
    -------
    profile : `~gammapy.estimators.FluxPoints`
        Profile flux points.
    """
    datasets = Datasets(datasets=datasets)
    maps = []
    for region in self.regions:
        datasets_to_fit = datasets.to_spectrum_datasets(region=region)
        for dataset_spec, dataset_map in zip(datasets_to_fit, datasets):
            dataset_spec.background.data = (
                dataset_map.npred()
                .to_region_nd_map(
                    region, func=np.sum, weights=dataset_map.mask_safe
                )
                .data
            )
        datasets_to_fit.models = SkyModel(self.spectrum, name="test-source")
        fp = super().METHOD_NAME(datasets_to_fit)
        maps.append(fp)
    return FluxPoints.from_stack(
        maps=maps,
        axis=self.projected_distance_axis,
    )