prompt: 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()
  
  ```
  
golden_diff: 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)
  
verification_info: {"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"}]}