We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Models are lost in FluxProfileEstimator
**Gammapy version**
gammapy 1.0.1

**Bug description**
All models attached to my datasets are discarded when I calculate a flux profile with `FluxProfileEstimator`.

**Expected behavior**
The predicted counts from the models should become part of the background when `to_spectrum_datasets` is called in `FluxProfileEstimator.run()`




--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gammapy/estimators/points/profile.py`
Content:
```
1 # Licensed under a 3-clause BSD style license - see LICENSE.rst
2 """Tools to create profiles (i.e. 1D "slices" from 2D images)."""
3 from astropy import units as u
4 from regions import CircleAnnulusSkyRegion
5 from gammapy.datasets import Datasets
6 from gammapy.maps import MapAxis
7 from gammapy.modeling.models import PowerLawSpectralModel, SkyModel
8 from .core import FluxPoints
9 from .sed import FluxPointsEstimator
10 
11 __all__ = ["FluxProfileEstimator"]
12 
13 
14 class FluxProfileEstimator(FluxPointsEstimator):
15     """Estimate flux profiles
16 
17     Parameters
18     ----------
19     regions : list of `~regions.SkyRegion`
20         regions to use
21     spectrum : `~gammapy.modeling.models.SpectralModel` (optional)
22         Spectral model to compute the fluxes or brightness.
23         Default is power-law with spectral index of 2.
24     **kwargs : dict
25         Keywords forwarded to the `FluxPointsEstimator` (see documentation
26         there for further description of valid keywords)
27 
28     Examples
29     --------
30     This example shows how to compute a counts profile for the Fermi galactic
31     center region::
32 
33     >>> from astropy import units as u
34     >>> from astropy.coordinates import SkyCoord
35     >>> from gammapy.data import GTI
36     >>> from gammapy.estimators import FluxProfileEstimator
37     >>> from gammapy.utils.regions import make_orthogonal_rectangle_sky_regions
38     >>> from gammapy.datasets import MapDataset
39     >>> from gammapy.maps import RegionGeom
40 
41     >>> # load example data
42     >>> filename = "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz"
43     >>> dataset = MapDataset.read(filename, name="fermi-dataset")
44 
45     >>> # configuration
46     >>> dataset.gti = GTI.create("0s", "1e7s", "2010-01-01")
47 
48     >>> # creation of the boxes and axis
49     >>> start_pos = SkyCoord("-1d", "0d", frame='galactic')
50     >>> end_pos = SkyCoord("1d", "0d", frame='galactic')
51 
52     >>> regions = make_orthogonal_rectangle_sky_regions(
53                 start_pos=start_pos,
54                 end_pos=end_pos,
55                 wcs=dataset.counts.geom.wcs,
56                 height=2 * u.deg,
57                 nbin=21
58             )
59 
60     >>> # set up profile estimator and run
61     >>> prof_maker = FluxProfileEstimator(regions=regions, energy_edges=[10, 2000] * u.GeV)
62     >>> fermi_prof = prof_maker.run(dataset)
63     >>> print(fermi_prof)
64     FluxPoints
65     ----------
66     <BLANKLINE>
67       geom                   : RegionGeom
68       axes                   : ['lon', 'lat', 'energy', 'projected-distance']
69       shape                  : (1, 1, 1, 21)
70       quantities             : ['norm', 'norm_err', 'ts', 'npred', 'npred_excess', 'stat', 'counts', 'success']  # noqa: E501
71       ref. model             : pl
72       n_sigma                : 1
73       n_sigma_ul             : 2
74       sqrt_ts_threshold_ul   : 2
75       sed type init          : likelihood
76 
77     """
78 
79     tag = "FluxProfileEstimator"
80 
81     def __init__(self, regions, spectrum=None, **kwargs):
82         if len(regions) <= 1:
83             raise ValueError(
84                 "Please provide at least two regions for flux profile estimation."
85             )
86 
87         self.regions = regions
88 
89         if spectrum is None:
90             spectrum = PowerLawSpectralModel()
91 
92         self.spectrum = spectrum
93         super().__init__(**kwargs)
94 
95     @property
96     def projected_distance_axis(self):
97         """Get projected distance from the first region.
98 
99         For normal region this is defined as the distance from the
100         center of the region. For annulus shaped regions it is the
101         mean between the inner and outer radius.
102 
103         Returns
104         -------
105         axis : `MapAxis`
106             Projected distance axis
107         """
108         distances = []
109         center = self.regions[0].center
110 
111         for idx, region in enumerate(self.regions):
112             if isinstance(region, CircleAnnulusSkyRegion):
113                 distance = (region.inner_radius + region.outer_radius) / 2.0
114             else:
115                 distance = center.separation(region.center)
116 
117             distances.append(distance)
118 
119         return MapAxis.from_nodes(
120             u.Quantity(distances, "deg"), name="projected-distance"
121         )
122 
123     def run(self, datasets):
124         """Run flux profile estimation
125 
126         Parameters
127         ----------
128         datasets : list of `~gammapy.datasets.MapDataset`
129             Map datasets.
130 
131         Returns
132         -------
133         profile : `~gammapy.estimators.FluxPoints`
134             Profile flux points.
135         """
136         datasets = Datasets(datasets=datasets)
137 
138         maps = []
139 
140         for region in self.regions:
141             datasets_to_fit = datasets.to_spectrum_datasets(region=region)
142             datasets_to_fit.models = SkyModel(self.spectrum, name="test-source")
143             fp = super().run(datasets_to_fit)
144             maps.append(fp)
145 
146         return FluxPoints.from_stack(
147             maps=maps,
148             axis=self.projected_distance_axis,
149         )
150 
151     @property
152     def config_parameters(self):
153         """Config parameters"""
154         pars = self.__dict__.copy()
155         pars = {key.strip("_"): value for key, value in pars.items()}
156         pars.pop("regions")
157         return pars
158
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gammapy/estimators/points/profile.py b/gammapy/estimators/points/profile.py
--- a/gammapy/estimators/points/profile.py
+++ b/gammapy/estimators/points/profile.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """Tools to create profiles (i.e. 1D "slices" from 2D images)."""
+import numpy as np
 from astropy import units as u
 from regions import CircleAnnulusSkyRegion
 from gammapy.datasets import Datasets
@@ -133,12 +134,20 @@
         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().run(datasets_to_fit)
             maps.append(fp)

{"golden_diff": "diff --git a/gammapy/estimators/points/profile.py b/gammapy/estimators/points/profile.py\n--- a/gammapy/estimators/points/profile.py\n+++ b/gammapy/estimators/points/profile.py\n@@ -1,5 +1,6 @@\n # Licensed under a 3-clause BSD style license - see LICENSE.rst\n \"\"\"Tools to create profiles (i.e. 1D \"slices\" from 2D images).\"\"\"\n+import numpy as np\n from astropy import units as u\n from regions import CircleAnnulusSkyRegion\n from gammapy.datasets import Datasets\n@@ -133,12 +134,20 @@\n         profile : `~gammapy.estimators.FluxPoints`\n             Profile flux points.\n         \"\"\"\n+\n         datasets = Datasets(datasets=datasets)\n \n         maps = []\n-\n         for region in self.regions:\n             datasets_to_fit = datasets.to_spectrum_datasets(region=region)\n+            for dataset_spec, dataset_map in zip(datasets_to_fit, datasets):\n+                dataset_spec.background.data = (\n+                    dataset_map.npred()\n+                    .to_region_nd_map(\n+                        region, func=np.sum, weights=dataset_map.mask_safe\n+                    )\n+                    .data\n+                )\n             datasets_to_fit.models = SkyModel(self.spectrum, name=\"test-source\")\n             fp = super().run(datasets_to_fit)\n             maps.append(fp)\n", "issue": "Models are lost in FluxProfileEstimator\n**Gammapy version**\r\ngammapy 1.0.1\r\n\r\n**Bug description**\r\nAll models attached to my datasets are discarded when I calculate a flux profile with `FluxProfileEstimator`.\r\n\r\n**Expected behavior**\r\nThe predicted counts from the models should become part of the background when `to_spectrum_datasets` is called in `FluxProfileEstimator.run()`\r\n\r\n\r\n\n", "before_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\n\"\"\"Tools to create profiles (i.e. 1D \"slices\" from 2D images).\"\"\"\nfrom astropy import units as u\nfrom regions import CircleAnnulusSkyRegion\nfrom gammapy.datasets import Datasets\nfrom gammapy.maps import MapAxis\nfrom gammapy.modeling.models import PowerLawSpectralModel, SkyModel\nfrom .core import FluxPoints\nfrom .sed import FluxPointsEstimator\n\n__all__ = [\"FluxProfileEstimator\"]\n\n\nclass FluxProfileEstimator(FluxPointsEstimator):\n    \"\"\"Estimate flux profiles\n\n    Parameters\n    ----------\n    regions : list of `~regions.SkyRegion`\n        regions to use\n    spectrum : `~gammapy.modeling.models.SpectralModel` (optional)\n        Spectral model to compute the fluxes or brightness.\n        Default is power-law with spectral index of 2.\n    **kwargs : dict\n        Keywords forwarded to the `FluxPointsEstimator` (see documentation\n        there for further description of valid keywords)\n\n    Examples\n    --------\n    This example shows how to compute a counts profile for the Fermi galactic\n    center region::\n\n    >>> from astropy import units as u\n    >>> from astropy.coordinates import SkyCoord\n    >>> from gammapy.data import GTI\n    >>> from gammapy.estimators import FluxProfileEstimator\n    >>> from gammapy.utils.regions import make_orthogonal_rectangle_sky_regions\n    >>> from gammapy.datasets import MapDataset\n    >>> from gammapy.maps import RegionGeom\n\n    >>> # load example data\n    >>> filename = \"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz\"\n    >>> dataset = MapDataset.read(filename, name=\"fermi-dataset\")\n\n    >>> # configuration\n    >>> dataset.gti = GTI.create(\"0s\", \"1e7s\", \"2010-01-01\")\n\n    >>> # creation of the boxes and axis\n    >>> start_pos = SkyCoord(\"-1d\", \"0d\", frame='galactic')\n    >>> end_pos = SkyCoord(\"1d\", \"0d\", frame='galactic')\n\n    >>> regions = make_orthogonal_rectangle_sky_regions(\n                start_pos=start_pos,\n                end_pos=end_pos,\n                wcs=dataset.counts.geom.wcs,\n                height=2 * u.deg,\n                nbin=21\n            )\n\n    >>> # set up profile estimator and run\n    >>> prof_maker = FluxProfileEstimator(regions=regions, energy_edges=[10, 2000] * u.GeV)\n    >>> fermi_prof = prof_maker.run(dataset)\n    >>> print(fermi_prof)\n    FluxPoints\n    ----------\n    <BLANKLINE>\n      geom                   : RegionGeom\n      axes                   : ['lon', 'lat', 'energy', 'projected-distance']\n      shape                  : (1, 1, 1, 21)\n      quantities             : ['norm', 'norm_err', 'ts', 'npred', 'npred_excess', 'stat', 'counts', 'success']  # noqa: E501\n      ref. model             : pl\n      n_sigma                : 1\n      n_sigma_ul             : 2\n      sqrt_ts_threshold_ul   : 2\n      sed type init          : likelihood\n\n    \"\"\"\n\n    tag = \"FluxProfileEstimator\"\n\n    def __init__(self, regions, spectrum=None, **kwargs):\n        if len(regions) <= 1:\n            raise ValueError(\n                \"Please provide at least two regions for flux profile estimation.\"\n            )\n\n        self.regions = regions\n\n        if spectrum is None:\n            spectrum = PowerLawSpectralModel()\n\n        self.spectrum = spectrum\n        super().__init__(**kwargs)\n\n    @property\n    def projected_distance_axis(self):\n        \"\"\"Get projected distance from the first region.\n\n        For normal region this is defined as the distance from the\n        center of the region. For annulus shaped regions it is the\n        mean between the inner and outer radius.\n\n        Returns\n        -------\n        axis : `MapAxis`\n            Projected distance axis\n        \"\"\"\n        distances = []\n        center = self.regions[0].center\n\n        for idx, region in enumerate(self.regions):\n            if isinstance(region, CircleAnnulusSkyRegion):\n                distance = (region.inner_radius + region.outer_radius) / 2.0\n            else:\n                distance = center.separation(region.center)\n\n            distances.append(distance)\n\n        return MapAxis.from_nodes(\n            u.Quantity(distances, \"deg\"), name=\"projected-distance\"\n        )\n\n    def run(self, datasets):\n        \"\"\"Run flux profile estimation\n\n        Parameters\n        ----------\n        datasets : list of `~gammapy.datasets.MapDataset`\n            Map datasets.\n\n        Returns\n        -------\n        profile : `~gammapy.estimators.FluxPoints`\n            Profile flux points.\n        \"\"\"\n        datasets = Datasets(datasets=datasets)\n\n        maps = []\n\n        for region in self.regions:\n            datasets_to_fit = datasets.to_spectrum_datasets(region=region)\n            datasets_to_fit.models = SkyModel(self.spectrum, name=\"test-source\")\n            fp = super().run(datasets_to_fit)\n            maps.append(fp)\n\n        return FluxPoints.from_stack(\n            maps=maps,\n            axis=self.projected_distance_axis,\n        )\n\n    @property\n    def config_parameters(self):\n        \"\"\"Config parameters\"\"\"\n        pars = self.__dict__.copy()\n        pars = {key.strip(\"_\"): value for key, value in pars.items()}\n        pars.pop(\"regions\")\n        return pars\n", "path": "gammapy/estimators/points/profile.py"}], "after_files": [{"content": "# Licensed under a 3-clause BSD style license - see LICENSE.rst\n\"\"\"Tools to create profiles (i.e. 1D \"slices\" from 2D images).\"\"\"\nimport numpy as np\nfrom astropy import units as u\nfrom regions import CircleAnnulusSkyRegion\nfrom gammapy.datasets import Datasets\nfrom gammapy.maps import MapAxis\nfrom gammapy.modeling.models import PowerLawSpectralModel, SkyModel\nfrom .core import FluxPoints\nfrom .sed import FluxPointsEstimator\n\n__all__ = [\"FluxProfileEstimator\"]\n\n\nclass FluxProfileEstimator(FluxPointsEstimator):\n    \"\"\"Estimate flux profiles\n\n    Parameters\n    ----------\n    regions : list of `~regions.SkyRegion`\n        regions to use\n    spectrum : `~gammapy.modeling.models.SpectralModel` (optional)\n        Spectral model to compute the fluxes or brightness.\n        Default is power-law with spectral index of 2.\n    **kwargs : dict\n        Keywords forwarded to the `FluxPointsEstimator` (see documentation\n        there for further description of valid keywords)\n\n    Examples\n    --------\n    This example shows how to compute a counts profile for the Fermi galactic\n    center region::\n\n    >>> from astropy import units as u\n    >>> from astropy.coordinates import SkyCoord\n    >>> from gammapy.data import GTI\n    >>> from gammapy.estimators import FluxProfileEstimator\n    >>> from gammapy.utils.regions import make_orthogonal_rectangle_sky_regions\n    >>> from gammapy.datasets import MapDataset\n    >>> from gammapy.maps import RegionGeom\n\n    >>> # load example data\n    >>> filename = \"$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz\"\n    >>> dataset = MapDataset.read(filename, name=\"fermi-dataset\")\n\n    >>> # configuration\n    >>> dataset.gti = GTI.create(\"0s\", \"1e7s\", \"2010-01-01\")\n\n    >>> # creation of the boxes and axis\n    >>> start_pos = SkyCoord(\"-1d\", \"0d\", frame='galactic')\n    >>> end_pos = SkyCoord(\"1d\", \"0d\", frame='galactic')\n\n    >>> regions = make_orthogonal_rectangle_sky_regions(\n                start_pos=start_pos,\n                end_pos=end_pos,\n                wcs=dataset.counts.geom.wcs,\n                height=2 * u.deg,\n                nbin=21\n            )\n\n    >>> # set up profile estimator and run\n    >>> prof_maker = FluxProfileEstimator(regions=regions, energy_edges=[10, 2000] * u.GeV)\n    >>> fermi_prof = prof_maker.run(dataset)\n    >>> print(fermi_prof)\n    FluxPoints\n    ----------\n    <BLANKLINE>\n      geom                   : RegionGeom\n      axes                   : ['lon', 'lat', 'energy', 'projected-distance']\n      shape                  : (1, 1, 1, 21)\n      quantities             : ['norm', 'norm_err', 'ts', 'npred', 'npred_excess', 'stat', 'counts', 'success']  # noqa: E501\n      ref. model             : pl\n      n_sigma                : 1\n      n_sigma_ul             : 2\n      sqrt_ts_threshold_ul   : 2\n      sed type init          : likelihood\n\n    \"\"\"\n\n    tag = \"FluxProfileEstimator\"\n\n    def __init__(self, regions, spectrum=None, **kwargs):\n        if len(regions) <= 1:\n            raise ValueError(\n                \"Please provide at least two regions for flux profile estimation.\"\n            )\n\n        self.regions = regions\n\n        if spectrum is None:\n            spectrum = PowerLawSpectralModel()\n\n        self.spectrum = spectrum\n        super().__init__(**kwargs)\n\n    @property\n    def projected_distance_axis(self):\n        \"\"\"Get projected distance from the first region.\n\n        For normal region this is defined as the distance from the\n        center of the region. For annulus shaped regions it is the\n        mean between the inner and outer radius.\n\n        Returns\n        -------\n        axis : `MapAxis`\n            Projected distance axis\n        \"\"\"\n        distances = []\n        center = self.regions[0].center\n\n        for idx, region in enumerate(self.regions):\n            if isinstance(region, CircleAnnulusSkyRegion):\n                distance = (region.inner_radius + region.outer_radius) / 2.0\n            else:\n                distance = center.separation(region.center)\n\n            distances.append(distance)\n\n        return MapAxis.from_nodes(\n            u.Quantity(distances, \"deg\"), name=\"projected-distance\"\n        )\n\n    def run(self, datasets):\n        \"\"\"Run flux profile estimation\n\n        Parameters\n        ----------\n        datasets : list of `~gammapy.datasets.MapDataset`\n            Map datasets.\n\n        Returns\n        -------\n        profile : `~gammapy.estimators.FluxPoints`\n            Profile flux points.\n        \"\"\"\n\n        datasets = Datasets(datasets=datasets)\n\n        maps = []\n        for region in self.regions:\n            datasets_to_fit = datasets.to_spectrum_datasets(region=region)\n            for dataset_spec, dataset_map in zip(datasets_to_fit, datasets):\n                dataset_spec.background.data = (\n                    dataset_map.npred()\n                    .to_region_nd_map(\n                        region, func=np.sum, weights=dataset_map.mask_safe\n                    )\n                    .data\n                )\n            datasets_to_fit.models = SkyModel(self.spectrum, name=\"test-source\")\n            fp = super().run(datasets_to_fit)\n            maps.append(fp)\n\n        return FluxPoints.from_stack(\n            maps=maps,\n            axis=self.projected_distance_axis,\n        )\n\n    @property\n    def config_parameters(self):\n        \"\"\"Config parameters\"\"\"\n        pars = self.__dict__.copy()\n        pars = {key.strip(\"_\"): value for key, value in pars.items()}\n        pars.pop(\"regions\")\n        return pars\n", "path": "gammapy/estimators/points/profile.py"}]}