Repository: gammapy/gammapy
Pull request #3852 (2022-03-15 18:13:55)

Analysis class cleanup to support future extension of the hli
Add the remaining changes proposed in #3788 related to the cleanup of the analysis class :
- the methods .get_xx are replace by specific AnalysisStep classes and can be defined in the configuration file (for example config.general.steps = ["data-reduction", "fit"]), or called interactively by analysis.run(["data-reduction", "fit"])
- in the cli gammapy analysis run is now equivalent to analysis.run() so it can execute all the analysis step available and not only data reduction.

codecov[bot]: # [Codecov](https://codecov.io/gh/gammapy/gammapy/pull/3852?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
> Merging [#3852](https://codecov.io/gh/gammapy/gammapy/pull/3852?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (4b085bc) into [master](https://codecov.io/gh/gammapy/gammapy/commit/ddbc5dc8a4291f125996a9ce0e9b18ecdcc4c6e5?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) (ddbc5dc) will **decrease** coverage by `0.01%`.
> The diff coverage is `94.76%`.

```diff
@@            Coverage Diff             @@
##           master    #3852      +/-   ##
==========================================
- Coverage   93.82%   93.81%   -0.02%     
==========================================
  Files         162      163       +1     
  Lines       19891    19988      +97     
==========================================
+ Hits        18663    18751      +88     
- Misses       1228     1237       +9     
```


| [Impacted Files](https://codecov.io/gh/gammapy/gammapy/pull/3852?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) | Coverage Δ | |
|---|---|---|
| [gammapy/analysis/core.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9hbmFseXNpcy9jb3JlLnB5) | `95.04% <92.85%> (-0.57%)` | :arrow_down: |
| [gammapy/analysis/steps.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9hbmFseXNpcy9zdGVwcy5weQ==) | `94.79% <94.79%> (ø)` | |
| [gammapy/analysis/\_\_init\_\_.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9hbmFseXNpcy9fX2luaXRfXy5weQ==) | `100.00% <100.00%> (ø)` | |
| [gammapy/analysis/config.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9hbmFseXNpcy9jb25maWcucHk=) | `100.00% <100.00%> (ø)` | |
| [gammapy/scripts/analysis.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9zY3JpcHRzL2FuYWx5c2lzLnB5) | `100.00% <100.00%> (ø)` | |
| [gammapy/maps/wcs/ndmap.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9tYXBzL3djcy9uZG1hcC5weQ==) | `91.02% <0.00%> (-0.50%)` | :arrow_down: |
| [gammapy/maps/core.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9tYXBzL2NvcmUucHk=) | `88.21% <0.00%> (ø)` | |
| [gammapy/maps/region/ndmap.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS9tYXBzL3JlZ2lvbi9uZG1hcC5weQ==) | `93.23% <0.00%> (ø)` | |
| [gammapy/utils/interpolation.py](https://codecov.io/gh/gammapy/gammapy/pull/3852/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy#diff-Z2FtbWFweS91dGlscy9pbnRlcnBvbGF0aW9uLnB5) | `94.11% <0.00%> (ø)` | |

:mega: Codecov can now indicate which changes are the most critical in Pull Requests. [Learn more](https://about.codecov.io/product/feature/runtime-insights/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy)

registerrier: This PR implements the basis of a large refactoring of the HLI. A number of design choices will have to be made and a full expanded and improved HLI is an objective for v2.0.
For now, even though the changes proposed in this PR do not break the API, we decided during the co-working week to postpone to v1.1 and implement the HLI in parallel with the existing one.

Patch:
diff --git a/gammapy/scripts/main.py b/gammapy/scripts/main.py
index b16c36aac..69e86062f 100644
--- a/gammapy/scripts/main.py
+++ b/gammapy/scripts/main.py
@@ -80,6 +80,23 @@ def cli_analysis():
     """
 
 
+@cli.group("workflow")
+def cli_workflow():
+    """Automation of configuration driven data reduction process.
+
+    \b
+    Examples
+    --------
+
+    \b
+    $ gammapy workflow config
+    $ gammapy workflow run
+    $ gammapy workflow config --overwrite
+    $ gammapy workflow config --filename myconfig.yaml
+    $ gammapy workflow run --filename myconfig.yaml
+    """
+
+
 @cli.group("download", short_help="Download datasets and notebooks")
 @click.pass_context
 def cli_download(ctx):  # noqa: D301
@@ -130,5 +147,13 @@ def add_subcommands():
 
     cli_analysis.add_command(cli_run_analysis)
 
+    from .workflow import cli_make_config
+
+    cli_workflow.add_command(cli_make_config)
+
+    from .workflow import cli_run_workflow
+
+    cli_workflow.add_command(cli_run_workflow)
+
 
 add_subcommands()
diff --git a/gammapy/scripts/workflow.py b/gammapy/scripts/workflow.py
new file mode 100644
index 000000000..6c3c878cf
--- /dev/null
+++ b/gammapy/scripts/workflow.py
@@ -0,0 +1,37 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import logging
+import click
+from gammapy.workflow import Workflow, WorkflowConfig
+
+log = logging.getLogger(__name__)
+
+
+@click.command(name="config")
+@click.option(
+    "--filename",
+    default="config.yaml",
+    help="Filename to store the default configuration values.",
+    show_default=True,
+)
+@click.option(
+    "--overwrite", default=False, is_flag=True, help="Overwrite existing file."
+)
+def cli_make_config(filename, overwrite):
+    """Writes default configuration file."""
+    config = WorkflowConfig()
+    config.write(filename, overwrite=overwrite)
+    log.info(f"Configuration file produced: {filename}")
+
+
+@click.command(name="run")
+@click.option(
+    "--filename",
+    default="config.yaml",
+    help="Filename with default configuration values.",
+    show_default=True,
+)
+def cli_run_workflow(filename):
+    """Performs automated data reduction process."""
+    config = WorkflowConfig.read(filename)
+    workflow = Workflow(config)
+    workflow.run()
diff --git a/gammapy/workflow/__init__.py b/gammapy/workflow/__init__.py
new file mode 100644
index 000000000..7025add12
--- /dev/null
+++ b/gammapy/workflow/__init__.py
@@ -0,0 +1,31 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Gammapy high level interface (workflow)."""
+from gammapy.utils.registry import Registry
+from .config import WorkflowConfig
+from .core import Workflow
+from .steps import (
+    DataSelectionWorkflowStep,
+    DatasetsWorkflowStep,
+    ExcessMapWorkflowStep,
+    FitWorkflowStep,
+    FluxPointsWorkflowStep,
+    LightCurveWorkflowStep,
+    ObservationsWorkflowStep,
+)
+
+__all__ = [
+    "Workflow",
+    "WorkflowConfig",
+]
+
+WORKFLOW_STEP_REGISTRY = Registry(
+    [
+        DataSelectionWorkflowStep,
+        ObservationsWorkflowStep,
+        DatasetsWorkflowStep,
+        ExcessMapWorkflowStep,
+        FitWorkflowStep,
+        FluxPointsWorkflowStep,
+        LightCurveWorkflowStep,
+    ]
+)
diff --git a/gammapy/workflow/config.py b/gammapy/workflow/config.py
new file mode 100644
index 000000000..e649dd45d
--- /dev/null
+++ b/gammapy/workflow/config.py
@@ -0,0 +1,340 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import html
+import json
+import logging
+from collections import defaultdict
+from collections.abc import Mapping
+from enum import Enum
+from pathlib import Path
+from typing import List, Optional
+from astropy import units as u
+import yaml
+from pydantic import BaseModel, ConfigDict
+from gammapy.makers import MapDatasetMaker
+from gammapy.utils.scripts import read_yaml, to_yaml, write_yaml
+from gammapy.utils.types import AngleType, EnergyType, PathType, TimeType
+
+__all__ = ["WorkflowConfig"]
+
+CONFIG_PATH = Path(__file__).resolve().parent / "config"
+DOCS_FILE = CONFIG_PATH / "docs.yaml"
+
+log = logging.getLogger(__name__)
+
+
+def deep_update(d, u):
+    """Recursively update a nested dictionary.
+
+    Taken from: https://stackoverflow.com/a/3233356/19802442
+    """
+    for k, v in u.items():
+        if isinstance(v, Mapping):
+            d[k] = deep_update(d.get(k, {}), v)
+        else:
+            d[k] = v
+    return d
+
+
+class WorkflowStepEnum(str, Enum):
+    data_reduction = "data-selection"
+    observations = "observations"
+    datasets = "datasets"
+    excess_map = "excess-map"
+    fit = "fit"
+    flux_points = "flux-points"
+    light_curve = "light-curve"
+
+
+class ReductionTypeEnum(str, Enum):
+    spectrum = "1d"
+    cube = "3d"
+
+
+class FrameEnum(str, Enum):
+    icrs = "icrs"
+    galactic = "galactic"
+
+
+class RequiredHDUEnum(str, Enum):
+    events = "events"
+    gti = "gti"
+    aeff = "aeff"
+    bkg = "bkg"
+    edisp = "edisp"
+    psf = "psf"
+    rad_max = "rad_max"
+
+
+class BackgroundMethodEnum(str, Enum):
+    reflected = "reflected"
+    fov = "fov_background"
+    ring = "ring"
+
+
+class SafeMaskMethodsEnum(str, Enum):
+    aeff_default = "aeff-default"
+    aeff_max = "aeff-max"
+    edisp_bias = "edisp-bias"
+    offset_max = "offset-max"
+    bkg_peak = "bkg-peak"
+
+
+class MapSelectionEnum(str, Enum):
+    counts = "counts"
+    exposure = "exposure"
+    background = "background"
+    psf = "psf"
+    edisp = "edisp"
+
+
+class GammapyBaseConfig(BaseModel):
+    model_config = ConfigDict(
+        arbitrary_types_allowed=True,
+        validate_assignment=True,
+        extra="forbid",
+        validate_default=True,
+        use_enum_values=True,
+        json_encoders={u.Quantity: lambda v: f"{v.value} {v.unit}"},
+    )
+
+    def _repr_html_(self):
+        try:
+            return self.to_html()
+        except AttributeError:
+            return f"<pre>{html.escape(str(self))}</pre>"
+
+
+class SkyCoordConfig(GammapyBaseConfig):
+    frame: Optional[FrameEnum] = None
+    lon: Optional[AngleType] = None
+    lat: Optional[AngleType] = None
+
+
+class EnergyAxisConfig(GammapyBaseConfig):
+    min: Optional[EnergyType] = None
+    max: Optional[EnergyType] = None
+    nbins: Optional[int] = None
+
+
+class SpatialCircleConfig(GammapyBaseConfig):
+    frame: Optional[FrameEnum] = None
+    lon: Optional[AngleType] = None
+    lat: Optional[AngleType] = None
+    radius: Optional[AngleType] = None
+
+
+class EnergyRangeConfig(GammapyBaseConfig):
+    min: Optional[EnergyType] = None
+    max: Optional[EnergyType] = None
+
+
+class TimeRangeConfig(GammapyBaseConfig):
+    start: Optional[TimeType] = None
+    stop: Optional[TimeType] = None
+
+
+class FluxPointsConfig(GammapyBaseConfig):
+    energy: EnergyAxisConfig = EnergyAxisConfig()
+    source: str = "source"
+    parameters: dict = {"selection_optional": "all"}
+
+
+class LightCurveConfig(GammapyBaseConfig):
+    time_intervals: TimeRangeConfig = TimeRangeConfig()
+    energy_edges: EnergyAxisConfig = EnergyAxisConfig()
+    source: str = "source"
+    parameters: dict = {"selection_optional": "all"}
+
+
+class FitConfig(GammapyBaseConfig):
+    fit_range: EnergyRangeConfig = EnergyRangeConfig()
+
+
+class ExcessMapConfig(GammapyBaseConfig):
+    correlation_radius: AngleType = "0.1 deg"
+    parameters: dict = {}
+    energy_edges: EnergyAxisConfig = EnergyAxisConfig()
+
+
+class BackgroundConfig(GammapyBaseConfig):
+    method: Optional[BackgroundMethodEnum] = None
+    exclusion: Optional[PathType] = None
+    parameters: dict = {}
+
+
+class SafeMaskConfig(GammapyBaseConfig):
+    methods: List[SafeMaskMethodsEnum] = [SafeMaskMethodsEnum.aeff_default]
+    parameters: dict = {}
+
+
+class EnergyAxesConfig(GammapyBaseConfig):
+    energy: EnergyAxisConfig = EnergyAxisConfig(min="1 TeV", max="10 TeV", nbins=5)
+    energy_true: EnergyAxisConfig = EnergyAxisConfig(
+        min="0.5 TeV", max="20 TeV", nbins=16
+    )
+
+
+class SelectionConfig(GammapyBaseConfig):
+    offset_max: AngleType = "2.5 deg"
+
+
+class WidthConfig(GammapyBaseConfig):
+    width: AngleType = "5 deg"
+    height: AngleType = "5 deg"
+
+
+class WcsConfig(GammapyBaseConfig):
+    skydir: SkyCoordConfig = SkyCoordConfig()
+    binsize: AngleType = "0.02 deg"
+    width: WidthConfig = WidthConfig()
+    binsize_irf: AngleType = "0.2 deg"
+
+
+class GeomConfig(GammapyBaseConfig):
+    wcs: WcsConfig = WcsConfig()
+    selection: SelectionConfig = SelectionConfig()
+    axes: EnergyAxesConfig = EnergyAxesConfig()
+
+
+class DatasetsConfig(GammapyBaseConfig):
+    type: ReductionTypeEnum = ReductionTypeEnum.spectrum
+    stack: bool = True
+    geom: GeomConfig = GeomConfig()
+    map_selection: List[MapSelectionEnum] = MapDatasetMaker.available_selection
+    background: BackgroundConfig = BackgroundConfig()
+    safe_mask: SafeMaskConfig = SafeMaskConfig()
+    on_region: SpatialCircleConfig = SpatialCircleConfig()
+    containment_correction: bool = True
+
+
+class ObservationsConfig(GammapyBaseConfig):
+    datastore: PathType = Path("$GAMMAPY_DATA/hess-dl3-dr1/")
+    obs_ids: List[int] = []
+    obs_file: Optional[PathType] = None
+    obs_cone: SpatialCircleConfig = SpatialCircleConfig()
+    obs_time: TimeRangeConfig = TimeRangeConfig()
+    required_irf: List[RequiredHDUEnum] = ["aeff", "edisp", "psf", "bkg"]
+
+
+class LogConfig(GammapyBaseConfig):
+    level: str = "info"
+    filename: Optional[PathType] = None
+    filemode: Optional[str] = None
+    format: Optional[str] = None
+    datefmt: Optional[str] = None
+
+
+class GeneralConfig(GammapyBaseConfig):
+    log: LogConfig = LogConfig()
+    outdir: str = "."
+    n_jobs: int = 1
+    datasets_file: Optional[PathType] = None
+    models_file: Optional[PathType] = None
+    steps: List[str] = ["data-selection", "fit", "flux-points"]
+    overwrite: List[bool] = [True, True, True]
+
+
+class WorkflowConfig(GammapyBaseConfig):
+    """Gammapy workflow configuration."""
+
+    general: GeneralConfig = GeneralConfig()
+    observations: ObservationsConfig = ObservationsConfig()
+    datasets: DatasetsConfig = DatasetsConfig()
+    fit: FitConfig = FitConfig()
+    flux_points: FluxPointsConfig = FluxPointsConfig()
+    excess_map: ExcessMapConfig = ExcessMapConfig()
+    light_curve: LightCurveConfig = LightCurveConfig()
+
+    def __str__(self):
+        """Display settings in pretty YAML format."""
+        info = self.__class__.__name__ + "\n\n\t"
+        data = self.to_yaml()
+        data = data.replace("\n", "\n\t")
+        info += data
+        return info.expandtabs(tabsize=4)
+
+    @classmethod
+    def read(cls, path):
+        """Read from YAML file.
+
+        Parameters
+        ----------
+        path : str
+            input filepath
+        """
+        config = read_yaml(path)
+        config.pop("metadata", None)
+        return WorkflowConfig(**config)
+
+    @classmethod
+    def from_yaml(cls, config_str):
+        """Create from YAML string.
+
+        Parameters
+        ----------
+        config_str : str
+            yaml str
+
+        """
+        settings = yaml.safe_load(config_str)
+        return WorkflowConfig(**settings)
+
+    def write(self, path, overwrite=False):
+        """Write to YAML file.
+
+        Parameters
+        ----------
+        path : `pathlib.Path` or str
+            Path to write files.
+        overwrite : bool, optional
+            Overwrite existing file. Default is False.
+        """
+        yaml_str = self.to_yaml()
+        write_yaml(yaml_str, path, overwrite=overwrite)
+
+    def to_yaml(self):
+        """Convert to YAML string."""
+        data = json.loads(self.model_dump_json())
+        return to_yaml(data)
+
+    def set_logging(self):
+        """Set logging config.
+
+        Calls ``logging.basicConfig``, i.e. adjusts global logging state.
+        """
+        self.general.log.level = self.general.log.level.upper()
+        logging.basicConfig(**self.general.log.model_dump())
+        log.info("Setting logging config: {!r}".format(self.general.log.model_dump()))
+
+    def update(self, config=None):
+        """Update config with provided settings.
+
+        Parameters
+        ----------
+        config : str or `WorkflowConfig` object, optional
+            Configuration settings provided in dict() syntax. Default is None.
+        """
+        if isinstance(config, str):
+            other = WorkflowConfig.from_yaml(config)
+        elif isinstance(config, WorkflowConfig):
+            other = config
+        else:
+            raise TypeError(f"Invalid type: {config}")
+
+        config_new = deep_update(
+            self.model_dump(exclude_defaults=True),
+            other.model_dump(exclude_defaults=True),
+        )
+        return WorkflowConfig(**config_new)
+
+    @staticmethod
+    def _get_doc_sections():
+        """Return dictionary with commented docs from docs file."""
+        doc = defaultdict(str)
+        with open(DOCS_FILE) as f:
+            for line in filter(lambda line: not line.startswith("---"), f):
+                line = line.strip("\n")
+                if line.startswith("# Section: "):
+                    keyword = line.replace("# Section: ", "")
+                doc[keyword] += line + "\n"
+        return doc
diff --git a/gammapy/workflow/config/docs.yaml b/gammapy/workflow/config/docs.yaml
new file mode 100644
index 000000000..26484401c
--- /dev/null
+++ b/gammapy/workflow/config/docs.yaml
@@ -0,0 +1,68 @@
+# Section: general
+# General settings for the high level interface / optional
+general:
+    # logging settings for the session
+    log:
+        # choose one of the example values for level
+        level: INFO            # also CRITICAL, ERROR, WARNING, DEBUG
+        filename: filename.log
+        filemode: w
+        format: "%(asctime)s - %(message)s"
+        datefmt: "%d-%b-%y %H:%M:%S"
+    # output folder where files will be stored
+    outdir: .
+
+# Section: observations
+# Observations used in the workflow / mandatory
+observations:
+    # path to data store where to fetch observations
+    datastore: $GAMMAPY_DATA/hess-dl3-dr1/
+    obs_ids: [23523, 23526]
+    obs_file:   # csv file with obs_ids
+    # spatial /time filters applied on the obs_ids
+    obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 3 deg}
+    obs_time: {start: '2019-12-01', stop: '2020-03-01'}
+
+# Section: datasets
+# Process of data reduction / mandatory
+datasets:
+    type: 3d   # also 1d
+    stack: false
+    geom:
+        wcs:
+            skydir: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg}
+            binsize: 0.1 deg
+            width: {width: 7 deg, height: 5 deg}
+            binsize_irf: 0.1 deg
+        axes:
+            energy: {min: 0.1 TeV, max: 10 TeV, nbins: 30}
+            energy_true: {min: 0.1 TeV, max: 10 TeV, nbins: 30}
+    map_selection: ['counts', 'exposure', 'background', 'psf', 'edisp']
+    background:
+        method: ring        # also fov_background, reflected for 1d
+        exclusion:          # fits file for exclusion mask
+        parameters: {r_in: 0.7 deg, width: 0.7 deg} # ring
+    safe_mask:
+            methods: ['aeff-default', 'offset-max']
+            parameters: {offset_max: 2.5 deg}
+    on_region: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 3 deg}
+    containment_correction: true
+
+# Section: fit
+# Fitting process / optional
+fit:
+    fit_range: {min: 0.1 TeV, max: 10 TeV}
+
+# Section: flux_points
+# Flux estimation process /optional
+flux_points:
+    energy: {min: 0.1 TeV, max: 10 TeV, nbins: 30}
+    source: "source"
+    parameters: {}
+
+# Section: excess_map
+# Flux/significance map process /optional
+excess_map:
+    correlation_radius: 0.1 deg
+    energy_edges: {min: 0.1 TeV, max: 10 TeV, nbins: 1}
+    parameters: {}
diff --git a/gammapy/workflow/config/example-1d.yaml b/gammapy/workflow/config/example-1d.yaml
new file mode 100644
index 000000000..7ee150865
--- /dev/null
+++ b/gammapy/workflow/config/example-1d.yaml
@@ -0,0 +1,19 @@
+observations:
+    datastore: $GAMMAPY_DATA/hess-dl3-dr1
+    obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 5 deg}
+    obs_ids: [23592, 23559]
+
+datasets:
+    type: 1d
+    stack: false
+    geom:
+        axes:
+            energy: {min: 0.01 TeV, max: 100 TeV, nbins: 73}
+    on_region: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 0.1 deg}
+    containment_correction: true
+
+fit:
+    fit_range: {min: 0.1 TeV, max: 100 TeV}
+
+flux_points:
+    energy: {min: 1 TeV, max: 10 TeV, nbins: 2}
diff --git a/gammapy/workflow/config/example-3d.yaml b/gammapy/workflow/config/example-3d.yaml
new file mode 100644
index 000000000..06caf9f96
--- /dev/null
+++ b/gammapy/workflow/config/example-3d.yaml
@@ -0,0 +1,37 @@
+observations:
+    datastore: $GAMMAPY_DATA/hess-dl3-dr1
+    obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 5 deg}
+    obs_ids: [23592, 23559]
+
+datasets:
+    type: 3d
+    stack: true
+    geom:
+        wcs:
+            skydir: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg}
+            binsize: 0.04 deg
+            width: {width: 5 deg, height: 5 deg}
+            binsize_irf: 0.2 deg
+        selection:
+            offset_max: 2.5 deg
+        axes:
+            energy: {min: 1 TeV, max: 10 TeV, nbins: 4}
+            energy_true: {min: 0.7 TeV, max: 12 TeV, nbins: 10}
+    map_selection: ['counts', 'exposure', 'background', 'psf', 'edisp']
+    safe_mask:
+        methods: ['offset-max']
+        parameters: {offset_max: 2.5 deg}
+    background:
+        method: fov_background
+        parameters: {method: 'scale'}
+
+fit:
+    fit_range: {min: 1 TeV, max: 10.1 TeV}
+
+flux_points:
+    energy: {min: 1 TeV, max: 10 TeV, nbins: 2}
+
+excess_map:
+    correlation_radius: 0.1 deg
+    energy_edges: {min: 1 TeV, max: 10 TeV, nbins: 1}
+    parameters: {}
diff --git a/gammapy/workflow/config/model-1d.yaml b/gammapy/workflow/config/model-1d.yaml
new file mode 100644
index 000000000..91149162b
--- /dev/null
+++ b/gammapy/workflow/config/model-1d.yaml
@@ -0,0 +1,21 @@
+components:
+    - name: source
+      type: SkyModel
+      spectral:
+          type: PowerLawSpectralModel
+          parameters:
+              - frozen: false
+                name: index
+                scale: 1.0
+                unit: ''
+                value: 2.6
+              - frozen: false
+                name: amplitude
+                scale: 1.0
+                unit: cm-2 s-1 TeV-1
+                value: 5.0e-11
+              - frozen: true
+                name: reference
+                scale: 1.0
+                unit: TeV
+                value: 1.0
diff --git a/gammapy/workflow/config/model.yaml b/gammapy/workflow/config/model.yaml
new file mode 100644
index 000000000..15a4074d2
--- /dev/null
+++ b/gammapy/workflow/config/model.yaml
@@ -0,0 +1,38 @@
+components:
+    - name: source
+      spatial:
+          parameters:
+              - frozen: false
+                max: 85
+                min: 80
+                name: lon_0
+                scale: 1.0
+                unit: deg
+                value: 83.633
+              - frozen: false
+                max: 24
+                min: 20
+                name: lat_0
+                scale: 1.0
+                unit: deg
+                value: 22.14
+          type: PointSpatialModel
+      spectral:
+          parameters:
+              - frozen: false
+                name: index
+                scale: 1.0
+                unit: ''
+                value: 2.6
+              - frozen: false
+                name: amplitude
+                scale: 1.0
+                unit: cm-2 s-1 TeV-1
+                value: 5.0e-11
+              - frozen: true
+                name: reference
+                scale: 1.0
+                unit: TeV
+                value: 1.0
+          type: PowerLawSpectralModel
+      type: SkyModel
diff --git a/gammapy/workflow/core.py b/gammapy/workflow/core.py
new file mode 100644
index 000000000..5463484de
--- /dev/null
+++ b/gammapy/workflow/core.py
@@ -0,0 +1,228 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Session class driving the high level interface API."""
+import html
+import logging
+from gammapy.datasets import Datasets
+from gammapy.modeling import Fit
+from gammapy.modeling.models import DatasetModels, FoVBackgroundModel, Models
+from gammapy.utils.scripts import make_path
+from gammapy.workflow.config import WorkflowConfig
+from .steps import WorkflowStep
+
+__all__ = ["Workflow"]
+
+log = logging.getLogger(__name__)
+
+
+class Workflow:
+    """Config-driven high level workflow interface.
+
+    It is initialized by default with a set of configuration parameters and values declared in
+    an internal high level interface model, though the user can also provide configuration
+    parameters passed as a nested dictionary at the moment of instantiation. In that case these
+    parameters will overwrite the default values of those present in the configuration file.
+
+    Parameters
+    ----------
+    config : dict or `~gammapy.workflow.WorkflowConfig`
+        Configuration options following `WorkflowConfig` schema.
+    """
+
+    def __init__(self, config):
+        self.log = log
+        self.config = config
+        self.config.set_logging()
+        self.datastore = None
+        self.observations = None
+        self.datasets = None
+        self.fit = Fit()
+        self.fit_result = None
+        self.flux_points = None
+
+    def _repr_html_(self):
+        try:
+            return self.to_html()
+        except AttributeError:
+            return f"<pre>{html.escape(str(self))}</pre>"
+
+    @property
+    def models(self):
+        if not self.datasets:
+            raise RuntimeError("No datasets defined. Impossible to set models.")
+        return self.datasets.models
+
+    @models.setter
+    def models(self, models):
+        self.set_models(models, extend=False)
+
+    @property
+    def config(self):
+        """Workflow configuration as an `~gammapy.workflow.WorkflowConfig` object."""
+        return self._config
+
+    @config.setter
+    def config(self, value):
+        if isinstance(value, dict):
+            self._config = WorkflowConfig(**value)
+        elif isinstance(value, WorkflowConfig):
+            self._config = value
+        else:
+            raise TypeError("config must be dict or WorkflowConfig.")
+
+    def run(self, steps=None, overwrite=None, **kwargs):
+        if steps is None:
+            steps = self.config.general.steps
+            overwrite = self.config.general.overwrite
+        else:
+            if overwrite is None:
+                overwrite = [True] * len(steps)
+        for k, step in enumerate(steps):
+            if isinstance(overwrite, list):
+                overwrite_step = overwrite[k]
+            else:
+                overwrite_step = overwrite
+            workflow_step = WorkflowStep.create(
+                step, self.config, log=self.log, overwrite=overwrite_step, **kwargs
+            )
+            workflow_step.run(self)
+
+    # keep these methods to be backward compatible
+    def get_observations(self):
+        """Fetch observations from the data store according to criteria defined in the configuration."""
+        self.run(steps=["observations"])
+
+    def get_datasets(self):
+        """Produce reduced datasets."""
+        self.run(steps=["datasets"])
+
+    def run_fit(self):
+        """Fitting reduced datasets to model."""
+        self.run(steps=["fit"])
+
+    def get_flux_points(self):
+        """Calculate flux points for a specific model component."""
+        self.run(steps=["flux-points"])
+
+    def get_excess_map(self):
+        """Calculate excess map with respect to the current model."""
+        self.run(steps=["excess-map"])
+
+    def get_light_curve(self):
+        """Calculate light curve for a specific model component."""
+        self.run(steps=["light-curve"])
+
+    def set_models(self, models, extend=True):
+        """Set models on datasets.
+
+        Adds `FoVBackgroundModel` if not present already
+
+        Parameters
+        ----------
+        models : `~gammapy.modeling.models.Models` or str
+            Models object or YAML models string.
+        extend : bool, optional
+            Extend the exiting models on the datasets or replace them.
+            Default is True.
+        """
+        if not self.datasets or len(self.datasets) == 0:
+            raise RuntimeError("Missing datasets")
+
+        self.log.info("Reading model.")
+        if isinstance(models, str):
+            models = Models.from_yaml(models)
+        elif isinstance(models, Models):
+            pass
+        elif isinstance(models, DatasetModels) or isinstance(models, list):
+            models = Models(models)
+        else:
+            raise TypeError(f"Invalid type: {models!r}")
+
+        if extend:
+            models.extend(self.datasets.models)
+
+        self.datasets.models = models
+
+        bkg_models = []
+        for dataset in self.datasets:
+            if dataset.tag == "MapDataset" and dataset.background_model is None:
+                bkg_models.append(FoVBackgroundModel(dataset_name=dataset.name))
+        if bkg_models:
+            models.extend(bkg_models)
+            self.datasets.models = models
+
+        self.log.info(models)
+
+    def read_models(self, path, extend=True):
+        """Read models from YAML file.
+
+        Parameters
+        ----------
+        path : str
+            Path to the model file.
+        extend : bool, optional
+            Extend the exiting models on the datasets or replace them.
+            Default is True.
+        """
+        path = make_path(path)
+        models = Models.read(path)
+        self.set_models(models, extend=extend)
+        self.log.info(f"Models loaded from {path}.")
+
+    def write_models(self, overwrite=True, write_covariance=True):
+        """Write models to YAML file.
+
+        File name is taken from the configuration file.
+        """
+        filename_models = self.config.general.models_file
+        if filename_models is not None:
+            self.models.write(
+                filename_models, overwrite=overwrite, write_covariance=write_covariance
+            )
+            log.info(f"Models loaded from {filename_models}.")
+        else:
+            raise RuntimeError("Missing models_file in config.general")
+
+    def read_datasets(self):
+        """Read datasets from YAML file.
+
+        File names are taken from the configuration file.
+        """
+        filename = self.config.general.datasets_file
+        filename_models = self.config.general.models_file
+        if filename is not None:
+            self.datasets = Datasets.read(filename)
+            log.info(f"Datasets loaded from {filename}.")
+            if filename_models is not None:
+                self.read_models(filename_models, extend=False)
+        else:
+            raise RuntimeError("Missing datasets_file in config.general")
+
+    def write_datasets(self, overwrite=True, write_covariance=True):
+        """Write datasets to YAML file.
+
+        File names are taken from the configuration file.
+
+        Parameters
+        ----------
+        overwrite : bool, optional
+            Overwrite existing file. Default is True.
+        write_covariance : bool, optional
+            Save covariance or not. Default is True.
+        """
+        filename = self.config.general.datasets_file
+        filename_models = self.config.general.models_file
+        if filename is not None:
+            self.datasets.write(
+                filename,
+                filename_models,
+                overwrite=overwrite,
+                write_covariance=write_covariance,
+            )
+            log.info(f"Datasets stored to {filename}.")
+            log.info(f"Datasets stored to {filename_models}.")
+        else:
+            raise RuntimeError("Missing datasets_file in config.general")
+
+    def update_config(self, config):
+        """Update the configuration."""
+        self.config = self.config.update(config=config)
diff --git a/gammapy/workflow/steps.py b/gammapy/workflow/steps.py
new file mode 100644
index 000000000..2833d8073
--- /dev/null
+++ b/gammapy/workflow/steps.py
@@ -0,0 +1,522 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""classes containing the workflow steps supported by the high level interface"""
+
+import abc
+import logging
+import os
+from astropy.coordinates import SkyCoord
+from astropy.table import Table
+from regions import CircleSkyRegion
+from gammapy.data import DataStore
+from gammapy.datasets import Datasets, FluxPointsDataset, MapDataset, SpectrumDataset
+from gammapy.estimators import (
+    ExcessMapEstimator,
+    FluxPointsEstimator,
+    LightCurveEstimator,
+)
+from gammapy.makers import (
+    DatasetsMaker,
+    FoVBackgroundMaker,
+    MapDatasetMaker,
+    ReflectedRegionsBackgroundMaker,
+    RingBackgroundMaker,
+    SafeMaskMaker,
+    SpectrumDatasetMaker,
+)
+from gammapy.maps import Map, MapAxis, RegionGeom, WcsGeom
+from gammapy.utils.pbar import progress_bar
+from gammapy.utils.scripts import make_name, make_path
+
+
+class WorkflowStepBase(abc.ABC):
+    tag = "workflow-step"
+
+    def __init__(self, config, log=None, name=None, overwrite=True):
+        self.config = config
+        self.overwrite = overwrite
+        self._name = make_name(name)
+
+        if log is None:
+            log = logging.getLogger(__name__)
+        self.log = log
+
+    @property
+    def name(self):
+        return self._name
+
+    def run(self, data):
+        # TODO: for now data is an workflow class
+        # but it will be replaced by a new specific container
+        # and the type and origin of data required will be defined by config
+        self.data = data
+        self.products = self.data
+        # TODO: for now it is the workflow class
+        # but in future each step will declare the products they create
+        # and init a specific container
+        self._run()
+
+    @abc.abstractmethod
+    def _run(self):
+        pass
+
+
+class WorkflowStep:
+    "Create one of the workflow step class listed in the registry"
+
+    @staticmethod
+    def create(tag, config, **kwargs):
+        from . import WORKFLOW_STEP_REGISTRY
+
+        cls = WORKFLOW_STEP_REGISTRY.get_cls(tag)
+        return cls(config, **kwargs)
+
+
+class DataSelectionWorkflowStep(WorkflowStepBase):
+    tag = "data-selection"
+
+    def _run(self):
+        filepath = self.config.general.datasets_file
+        if filepath is not None and os.path.exists(filepath) and not self.overwrite:
+            self.data.read_datasets()
+        else:
+            obs_step = ObservationsWorkflowStep(
+                self.config, log=self.log, overwrite=self.overwrite
+            )
+            datasets_step = DatasetsWorkflowStep(
+                self.config, log=self.log, overwrite=self.overwrite
+            )
+            obs_step.run(self.data)
+            datasets_step.run(self.data)
+            self.data.write_datasets()
+
+
+class ObservationsWorkflowStep(WorkflowStepBase):
+    tag = "observations"
+
+    def _run(self):
+        """Fetch observations from the data store according to criteria defined in the configuration."""
+        observations_settings = self.config.observations
+        self.data.datastore = self._get_data_store()
+
+        self.log.info("Fetching observations.")
+        ids = self._make_obs_table_selection()
+        required_irf = observations_settings.required_irf
+        self.data.observations = self.data.datastore.get_observations(
+            ids, skip_missing=True, required_irf=required_irf
+        )
+
+        if observations_settings.obs_time.start is not None:
+            start = observations_settings.obs_time.start
+            stop = observations_settings.obs_time.stop
+            if len(start.shape) == 0:
+                time_intervals = [(start, stop)]
+            else:
+                time_intervals = [(tstart, tstop) for tstart, tstop in zip(start, stop)]
+            self.data.observations = self.data.observations.select_time(time_intervals)
+
+        self.log.info(f"Number of selected observations: {len(self.data.observations)}")
+
+        for obs in self.data.observations:
+            self.log.debug(obs)
+
+    def _get_data_store(self):
+        """Set the datastore on the Workflow object."""
+        path = make_path(self.config.observations.datastore)
+        if path.is_file():
+            self.log.debug(f"Setting datastore from file: {path}")
+            return DataStore.from_file(path)
+        elif path.is_dir():
+            self.log.debug(f"Setting datastore from directory: {path}")
+            return DataStore.from_dir(path)
+        else:
+            raise FileNotFoundError(f"Datastore not found: {path}")
+
+    def _make_obs_table_selection(self):
+        """Return list of obs_ids after filtering on datastore observation table."""
+        obs_settings = self.config.observations
+
+        # Reject configs with list of obs_ids and obs_file set at the same time
+        if len(obs_settings.obs_ids) and obs_settings.obs_file is not None:
+            raise ValueError(
+                "Values for both parameters obs_ids and obs_file are not accepted."
+            )
+
+        # First select input list of observations from obs_table
+        if len(obs_settings.obs_ids):
+            selected_obs_table = self.data.datastore.obs_table.select_obs_id(
+                obs_settings.obs_ids
+            )
+        elif obs_settings.obs_file is not None:
+            path = make_path(obs_settings.obs_file)
+            ids = list(Table.read(path, format="ascii", data_start=0).columns[0])
+            selected_obs_table = self.data.datastore.obs_table.select_obs_id(ids)
+        else:
+            selected_obs_table = self.data.datastore.obs_table
+
+        # Apply cone selection
+        if obs_settings.obs_cone.lon is not None:
+            cone = dict(
+                type="sky_circle",
+                frame=obs_settings.obs_cone.frame,
+                lon=obs_settings.obs_cone.lon,
+                lat=obs_settings.obs_cone.lat,
+                radius=obs_settings.obs_cone.radius,
+                border="0 deg",
+            )
+            selected_obs_table = selected_obs_table.select_observations(cone)
+
+        return selected_obs_table["OBS_ID"].tolist()
+
+
+class DatasetsWorkflowStep(WorkflowStepBase):
+    tag = "datasets"
+
+    def _run(self):
+        # TODO: check if exits and read else run and write
+        self.get_datasets()
+
+    def get_datasets(self):
+        """Produce reduced datasets."""
+        datasets_settings = self.config.datasets
+        if not self.data.observations or len(self.data.observations) == 0:
+            raise RuntimeError("No observations have been selected.")
+
+        if datasets_settings.type == "1d":
+            self._spectrum_extraction()
+        else:  # 3d
+            self._map_making()
+
+    @staticmethod
+    def _create_wcs_geometry(wcs_geom_settings, axes):
+        """Create the WCS geometry."""
+        geom_params = {}
+        skydir_settings = wcs_geom_settings.skydir
+        if skydir_settings.lon is not None:
+            skydir = SkyCoord(
+                skydir_settings.lon, skydir_settings.lat, frame=skydir_settings.frame
+            )
+            geom_params["skydir"] = skydir
+
+        if skydir_settings.frame in ["icrs", "galactic"]:
+            geom_params["frame"] = skydir_settings.frame
+        else:
+            raise ValueError(
+                f"Incorrect skydir frame: expect 'icrs' or 'galactic'. Got {skydir_settings.frame}"
+            )
+
+        geom_params["axes"] = axes
+        geom_params["binsz"] = wcs_geom_settings.binsize
+        width = wcs_geom_settings.width.width.to("deg").value
+        height = wcs_geom_settings.width.height.to("deg").value
+        geom_params["width"] = (width, height)
+
+        return WcsGeom.create(**geom_params)
+
+    @staticmethod
+    def _create_region_geometry(on_region_settings, axes):
+        """Create the region geometry."""
+        on_lon = on_region_settings.lon
+        on_lat = on_region_settings.lat
+        on_center = SkyCoord(on_lon, on_lat, frame=on_region_settings.frame)
+        on_region = CircleSkyRegion(on_center, on_region_settings.radius)
+
+        return RegionGeom.create(region=on_region, axes=axes)
+
+    def _create_geometry(self):
+        """Create the geometry."""
+        self.log.debug("Creating geometry.")
+        datasets_settings = self.config.datasets
+        geom_settings = datasets_settings.geom
+        axes = [make_energy_axis(geom_settings.axes.energy)]
+        if datasets_settings.type == "3d":
+            geom = self._create_wcs_geometry(geom_settings.wcs, axes)
+        elif datasets_settings.type == "1d":
+            geom = self._create_region_geometry(datasets_settings.on_region, axes)
+        else:
+            raise ValueError(
+                f"Incorrect dataset type. Expect '1d' or '3d'. Got {datasets_settings.type}."
+            )
+        return geom
+
+    def _create_reference_dataset(self, name=None):
+        """Create the reference dataset for the current workflow."""
+        self.log.debug("Creating target Dataset.")
+        geom = self._create_geometry()
+
+        geom_settings = self.config.datasets.geom
+        geom_irf = dict(energy_axis_true=None, binsz_irf=None)
+        if geom_settings.axes.energy_true.min is not None:
+            geom_irf["energy_axis_true"] = make_energy_axis(
+                geom_settings.axes.energy_true, name="energy_true"
+            )
+        if geom_settings.wcs.binsize_irf is not None:
+            geom_irf["binsz_irf"] = geom_settings.wcs.binsize_irf.to("deg").value
+
+        if self.config.datasets.type == "1d":
+            return SpectrumDataset.create(geom, name=name, **geom_irf)
+        else:
+            return MapDataset.create(geom, name=name, **geom_irf)
+
+    def _create_dataset_maker(self):
+        """Create the Dataset Maker."""
+        self.log.debug("Creating the target Dataset Maker.")
+
+        datasets_settings = self.config.datasets
+        if datasets_settings.type == "3d":
+            maker = MapDatasetMaker(selection=datasets_settings.map_selection)
+        elif datasets_settings.type == "1d":
+            maker_config = {}
+            if datasets_settings.containment_correction:
+                maker_config["containment_correction"] = (
+                    datasets_settings.containment_correction
+                )
+
+            maker_config["selection"] = ["counts", "exposure", "edisp"]
+
+            maker = SpectrumDatasetMaker(**maker_config)
+
+        return maker
+
+    def _create_safe_mask_maker(self):
+        """Create the SafeMaskMaker."""
+        self.log.debug("Creating the mask_safe Maker.")
+
+        safe_mask_selection = self.config.datasets.safe_mask.methods
+        safe_mask_settings = self.config.datasets.safe_mask.parameters
+        return SafeMaskMaker(methods=safe_mask_selection, **safe_mask_settings)
+
+    def _create_background_maker(self):
+        """Create the Background maker."""
+        self.log.info("Creating the background Maker.")
+
+        datasets_settings = self.config.datasets
+        bkg_maker_config = {}
+        if datasets_settings.background.exclusion:
+            path = make_path(datasets_settings.background.exclusion)
+            exclusion_mask = Map.read(path)
+            exclusion_mask.data = exclusion_mask.data.astype(bool)
+            bkg_maker_config["exclusion_mask"] = exclusion_mask
+        bkg_maker_config.update(datasets_settings.background.parameters)
+
+        bkg_method = datasets_settings.background.method
+
+        bkg_maker = None
+        if bkg_method == "fov_background":
+            self.log.debug(
+                f"Creating FoVBackgroundMaker with arguments {bkg_maker_config}"
+            )
+            bkg_maker = FoVBackgroundMaker(**bkg_maker_config)
+        elif bkg_method == "ring":
+            bkg_maker = RingBackgroundMaker(**bkg_maker_config)
+            self.log.debug(
+                f"Creating RingBackgroundMaker with arguments {bkg_maker_config}"
+            )
+            if datasets_settings.geom.axes.energy.nbins > 1:
+                raise ValueError(
+                    "You need to define a single-bin energy geometry for your dataset."
+                )
+        elif bkg_method == "reflected":
+            bkg_maker = ReflectedRegionsBackgroundMaker(**bkg_maker_config)
+            self.log.debug(
+                f"Creating ReflectedRegionsBackgroundMaker with arguments {bkg_maker_config}"
+            )
+        else:
+            self.log.warning("No background maker set. Check configuration.")
+        return bkg_maker
+
+    def _map_making(self):
+        """Make maps and datasets for 3d workflow"""
+        datasets_settings = self.config.datasets
+        offset_max = datasets_settings.geom.selection.offset_max
+
+        self.log.info("Creating reference dataset and makers.")
+        stacked = self._create_reference_dataset(name="stacked")
+
+        maker = self._create_dataset_maker()
+        maker_safe_mask = self._create_safe_mask_maker()
+        bkg_maker = self._create_background_maker()
+
+        makers = [maker, maker_safe_mask, bkg_maker]
+        makers = [maker for maker in makers if maker is not None]
+
+        self.log.info("Start the data reduction loop.")
+
+        datasets_maker = DatasetsMaker(
+            makers,
+            stack_datasets=datasets_settings.stack,
+            n_jobs=self.config.general.n_jobs,
+            cutout_mode="trim",
+            cutout_width=2 * offset_max,
+        )
+        self.products.datasets = datasets_maker.run(stacked, self.data.observations)
+
+    def _spectrum_extraction(self):
+        """Run all steps for the spectrum extraction."""
+        self.log.info("Reducing spectrum datasets.")
+        datasets_settings = self.config.datasets
+        dataset_maker = self._create_dataset_maker()
+        safe_mask_maker = self._create_safe_mask_maker()
+        bkg_maker = self._create_background_maker()
+
+        reference = self._create_reference_dataset()
+
+        datasets = []
+        for obs in progress_bar(self.data.observations, desc="Observations"):
+            self.log.debug(f"Processing observation {obs.obs_id}")
+            dataset = dataset_maker.run(reference.copy(), obs)
+            if bkg_maker is not None:
+                dataset = bkg_maker.run(dataset, obs)
+                if dataset.counts_off is None:
+                    self.log.debug(
+                        f"No OFF region found for observation {obs.obs_id}. Discarding."
+                    )
+                    continue
+            dataset = safe_mask_maker.run(dataset, obs)
+            self.log.debug(dataset)
+            datasets.append(dataset)
+        self.products.datasets = Datasets(datasets)
+
+        if datasets_settings.stack:
+            stacked = self.data.datasets.stack_reduce(name="stacked")
+            self.products.datasets = Datasets([stacked])
+
+
+def make_energy_axis(axis, name="energy"):
+    if axis.min is None or axis.max is None:
+        return None
+    elif axis.nbins is None or axis.nbins < 1:
+        return None
+    else:
+        return MapAxis.from_bounds(
+            name=name,
+            lo_bnd=axis.min.value,
+            hi_bnd=axis.max.to_value(axis.min.unit),
+            nbin=axis.nbins,
+            unit=axis.min.unit,
+            interp="log",
+            node_type="edges",
+        )
+
+
+class ExcessMapWorkflowStep(WorkflowStepBase):
+    tag = "excess-map"
+
+    def _run(self):
+        """Calculate excess map with respect to the current model."""
+        excess_settings = self.config.excess_map
+        # TODO: allow a list of kernel sizes
+        self.log.info("Computing excess maps.")
+
+        if self.config.datasets.type == "1d":
+            raise ValueError("Cannot compute excess map for 1D dataset")
+
+        # Here we could possibly stack the datasets if needed.
+        if len(self.data.datasets) > 1:
+            raise ValueError("Datasets must be stacked to compute the excess map")
+
+        energy_edges = make_energy_axis(excess_settings.energy_edges)
+        if energy_edges is not None:
+            energy_edges = energy_edges.edges
+
+        excess_map_estimator = ExcessMapEstimator(
+            correlation_radius=excess_settings.correlation_radius,
+            energy_edges=energy_edges,
+            **excess_settings.parameters,
+        )
+        self.products.excess_map = excess_map_estimator.run(self.data.datasets[0])
+
+
+class FitWorkflowStep(WorkflowStepBase):
+    tag = "fit"
+
+    def _run(self):
+        """Fitting reduced datasets to model."""
+        if not self.data.models:
+            raise RuntimeError("Missing models")
+
+        fit_settings = self.config.fit
+        for dataset in self.data.datasets:
+            if fit_settings.fit_range:
+                energy_min = fit_settings.fit_range.min
+                energy_max = fit_settings.fit_range.max
+                # TODO : add fit range in lon/lat
+                geom = dataset.counts.geom
+                dataset.mask_fit = geom.energy_mask(energy_min, energy_max)
+
+        self.log.info("Fitting datasets.")
+        result = self.data.fit.run(datasets=self.data.datasets)
+        self.products.fit_result = result
+        self.log.info(self.products.fit_result)
+
+
+class FluxPointsWorkflowStep(WorkflowStepBase):
+    tag = "flux-points"
+
+    def _run(self):
+        """Calculate flux points for a specific model component."""
+
+        if not self.data.datasets:
+            raise RuntimeError(
+                "No datasets defined. Impossible to compute flux points."
+            )
+
+        fp_settings = self.config.flux_points
+        self.log.info("Calculating flux points.")
+
+        energy_edges = make_energy_axis(fp_settings.energy).edges
+        flux_point_estimator = FluxPointsEstimator(
+            energy_edges=energy_edges,
+            source=fp_settings.source,
+            fit=self.data.fit,
+            **fp_settings.parameters,
+        )
+
+        fp = flux_point_estimator.run(datasets=self.data.datasets)
+
+        self.products.flux_points = FluxPointsDataset(
+            data=fp, models=self.data.models[fp_settings.source]
+        )
+        cols = ["e_ref", "dnde", "dnde_ul", "dnde_err", "sqrt_ts"]
+        table = self.products.flux_points.data.to_table(sed_type="dnde")
+        self.log.info("\n{}".format(table[cols]))
+
+
+class LightCurveWorkflowStep(WorkflowStepBase):
+    tag = "light-curve"
+
+    def _run(self):
+        """Calculate light curve for a specific model component."""
+        lc_settings = self.config.light_curve
+        self.log.info("Computing light curve.")
+        energy_edges = make_energy_axis(lc_settings.energy_edges).edges
+
+        if (
+            lc_settings.time_intervals.start is None
+            or lc_settings.time_intervals.stop is None
+        ):
+            self.log.info(
+                "Time intervals not defined. Extract light curve on datasets GTIs."
+            )
+            time_intervals = None
+        else:
+            time_intervals = [
+                (t1, t2)
+                for t1, t2 in zip(
+                    lc_settings.time_intervals.start, lc_settings.time_intervals.stop
+                )
+            ]
+
+        light_curve_estimator = LightCurveEstimator(
+            time_intervals=time_intervals,
+            energy_edges=energy_edges,
+            source=lc_settings.source,
+            fit=self.data.fit,
+            **lc_settings.parameters,
+        )
+        lc = light_curve_estimator.run(datasets=self.data.datasets)
+        self.products.light_curve = lc
+        self.log.info(
+            "\n{}".format(
+                self.products.light_curve.to_table(format="lightcurve", sed_type="flux")
+            )
+        )
diff --git a/setup.cfg b/setup.cfg
index cfddf7a93..2d01d0e6e 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -93,6 +93,7 @@ docs =
 
 [options.package_data]
 gammapy.analysis = config/*
+gammapy.workflow = config/*
 gammapy.irf.psf.tests = data/*
 gammapy.modeling.models.tests = data/*
 gammapy.catalog.tests = data/*


Test patch:
diff --git a/gammapy/scripts/tests/test_workflow.py b/gammapy/scripts/tests/test_workflow.py
new file mode 100644
index 000000000..0a22a8477
--- /dev/null
+++ b/gammapy/scripts/tests/test_workflow.py
@@ -0,0 +1,31 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from gammapy.scripts.main import cli
+from gammapy.utils.testing import requires_data, run_cli
+from ...workflow.tests.test_workflow import get_example_config
+
+
+def test_cli_workflow_config(tmp_path):
+    path_config = tmp_path / "config.yaml"
+    args = ["workflow", "config", f"--filename={path_config}"]
+    run_cli(cli, args)
+    assert path_config.exists()
+
+
+@requires_data()
+def test_cli_workflow_run(tmp_path):
+
+    path_config = tmp_path / "config.yaml"
+    path_datasets = tmp_path / "datasets.yaml"
+    config = get_example_config("1d")
+    config.datasets.background.method = "reflected"
+    config.general.datasets_file = str(path_datasets)
+    config.general.steps = ["data-selection"]
+
+    config.write(path_config)
+    args = [
+        "workflow",
+        "run",
+        f"--filename={path_config}",
+    ]
+    run_cli(cli, args)
+    assert path_datasets.exists()
diff --git a/gammapy/workflow/tests/__init__.py b/gammapy/workflow/tests/__init__.py
new file mode 100644
index 000000000..9dce85d06
--- /dev/null
+++ b/gammapy/workflow/tests/__init__.py
@@ -0,0 +1,1 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
diff --git a/gammapy/workflow/tests/test_config.py b/gammapy/workflow/tests/test_config.py
new file mode 100644
index 000000000..0173705ec
--- /dev/null
+++ b/gammapy/workflow/tests/test_config.py
@@ -0,0 +1,141 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from pathlib import Path
+import pytest
+from astropy.coordinates import Angle
+from astropy.time import Time
+from astropy.units import Quantity
+from pydantic import ValidationError
+from gammapy.utils.testing import assert_allclose
+from gammapy.workflow.config import GeneralConfig, WorkflowConfig
+
+CONFIG_PATH = Path(__file__).resolve().parent.parent / "config"
+DOC_FILE = CONFIG_PATH / "docs.yaml"
+
+
+def test_config_default_types():
+    config = WorkflowConfig()
+    assert config.observations.obs_cone.frame is None
+    assert config.observations.obs_cone.lon is None
+    assert config.observations.obs_cone.lat is None
+    assert config.observations.obs_cone.radius is None
+    assert config.observations.obs_time.start is None
+    assert config.observations.obs_time.stop is None
+    assert config.datasets.geom.wcs.skydir.frame is None
+    assert config.datasets.geom.wcs.skydir.lon is None
+    assert config.datasets.geom.wcs.skydir.lat is None
+    assert isinstance(config.datasets.geom.wcs.binsize, Angle)
+    assert isinstance(config.datasets.geom.wcs.binsize_irf, Angle)
+    assert isinstance(config.datasets.geom.axes.energy.min, Quantity)
+    assert isinstance(config.datasets.geom.axes.energy.max, Quantity)
+    assert isinstance(config.datasets.geom.axes.energy_true.min, Quantity)
+    assert isinstance(config.datasets.geom.axes.energy_true.max, Quantity)
+    assert isinstance(config.datasets.geom.selection.offset_max, Angle)
+    assert config.fit.fit_range.min is None
+    assert config.fit.fit_range.max is None
+    assert isinstance(config.excess_map.correlation_radius, Angle)
+    assert config.excess_map.energy_edges.min is None
+    assert config.excess_map.energy_edges.max is None
+    assert config.excess_map.energy_edges.nbins is None
+
+
+def test_config_not_default_types():
+    config = WorkflowConfig()
+    config.observations.obs_cone = {
+        "frame": "galactic",
+        "lon": "83.633 deg",
+        "lat": "22.014 deg",
+        "radius": "1 deg",
+    }
+    config.fit.fit_range = {"min": "0.1 TeV", "max": "10 TeV"}
+    assert config.observations.obs_cone.frame == "galactic"
+    assert isinstance(config.observations.obs_cone.lon, Angle)
+    assert isinstance(config.observations.obs_cone.lat, Angle)
+    assert isinstance(config.observations.obs_cone.radius, Angle)
+    config.observations.obs_time.start = "2019-12-01"
+    assert isinstance(config.observations.obs_time.start, Time)
+    with pytest.raises(ValueError):
+        config.flux_points.energy.min = "1 deg"
+    assert isinstance(config.fit.fit_range.min, Quantity)
+    assert isinstance(config.fit.fit_range.max, Quantity)
+
+
+def test_config_basics():
+    config = WorkflowConfig()
+    assert "WorkflowConfig" in str(config)
+    config = WorkflowConfig.read(DOC_FILE)
+    assert config.general.outdir == "."
+
+
+def test_config_create_from_dict():
+    data = {"general": {"log": {"level": "warning"}}}
+    config = WorkflowConfig(**data)
+    assert config.general.log.level == "warning"
+
+
+def test_config_create_from_yaml():
+    config = WorkflowConfig.read(DOC_FILE)
+    assert isinstance(config.general, GeneralConfig)
+    config_str = Path(DOC_FILE).read_text()
+    config = WorkflowConfig.from_yaml(config_str)
+    assert isinstance(config.general, GeneralConfig)
+
+
+def test_config_to_yaml(tmp_path):
+    config = WorkflowConfig()
+    assert "level: info" in config.to_yaml()
+    config = WorkflowConfig()
+    fpath = Path(tmp_path) / "temp.yaml"
+    config.write(fpath)
+    text = Path(fpath).read_text()
+    assert "stack" in text
+    with pytest.raises(IOError):
+        config.write(fpath)
+
+
+def test_get_doc_sections():
+    config = WorkflowConfig()
+    doc = config._get_doc_sections()
+    assert "general" in doc.keys()
+
+
+def test_safe_mask_config_validation():
+    config = WorkflowConfig()
+    # Check empty list is accepted
+    config.datasets.safe_mask.methods = []
+
+    with pytest.raises(ValidationError):
+        config.datasets.safe_mask.methods = ["bad"]
+
+
+def test_time_range_iso():
+    cfg = """
+    observations:
+        datastore: $GAMMAPY_DATA/hess-dl3-dr1
+        obs_ids: [23523, 23526]
+        obs_time: {
+            start: [2004-12-04 22:04:48.000, 2004-12-04 22:26:24.000, 2004-12-04 22:53:45.600],
+            stop: [2004-12-04 22:26:24.000, 2004-12-04 22:53:45.600, 2004-12-04 23:31:12.000]
+            }
+    """
+    config = WorkflowConfig.from_yaml(cfg)
+
+    assert_allclose(
+        config.observations.obs_time.start.mjd, [53343.92, 53343.935, 53343.954]
+    )
+
+
+def test_time_range_jyear():
+    cfg = """
+    observations:
+        datastore: $GAMMAPY_DATA/hess-dl3-dr1
+        obs_ids: [23523, 23526]
+        obs_time: {
+            start: [J2004.92654346, J2004.92658453, J2004.92663655],
+            stop: [J2004.92658453, J2004.92663655, J2004.92670773]
+            }
+    """
+    config = WorkflowConfig.from_yaml(cfg)
+
+    assert_allclose(
+        config.observations.obs_time.start.mjd, [53343.92, 53343.935, 53343.954]
+    )
diff --git a/gammapy/workflow/tests/test_workflow.py b/gammapy/workflow/tests/test_workflow.py
new file mode 100644
index 000000000..b2ca4165c
--- /dev/null
+++ b/gammapy/workflow/tests/test_workflow.py
@@ -0,0 +1,469 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import logging
+from pathlib import Path
+import pytest
+from numpy.testing import assert_allclose
+import astropy.units as u
+from astropy.coordinates import SkyCoord
+from regions import CircleSkyRegion
+from pydantic import ValidationError
+from gammapy.datasets import MapDataset, SpectrumDatasetOnOff
+from gammapy.maps import WcsGeom, WcsNDMap
+from gammapy.modeling.models import DatasetModels
+from gammapy.utils.testing import requires_data
+from gammapy.workflow import Workflow, WorkflowConfig
+
+CONFIG_PATH = Path(__file__).resolve().parent / ".." / "config"
+MODEL_FILE = CONFIG_PATH / "model.yaml"
+MODEL_FILE_1D = CONFIG_PATH / "model-1d.yaml"
+
+
+def get_example_config(which):
+    """Example config: which can be 1d or 3d."""
+    return WorkflowConfig.read(CONFIG_PATH / f"example-{which}.yaml")
+
+
+def test_init():
+    cfg = {"general": {"outdir": "test"}}
+    workflow = Workflow(cfg)
+    assert workflow.config.general.outdir == "test"
+    with pytest.raises(TypeError):
+        Workflow("spam")
+
+
+def test_update_config():
+    workflow = Workflow(WorkflowConfig())
+    data = {"general": {"outdir": "test"}}
+    config = WorkflowConfig(**data)
+    workflow.update_config(config)
+    assert workflow.config.general.outdir == "test"
+
+    workflow = Workflow(WorkflowConfig())
+    data = """
+    general:
+        outdir: test
+    """
+    workflow.update_config(data)
+    assert workflow.config.general.outdir == "test"
+
+    workflow = Workflow(WorkflowConfig())
+    with pytest.raises(TypeError):
+        workflow.update_config(0)
+
+
+def test_get_observations_no_datastore():
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.datastore = "other"
+    with pytest.raises(FileNotFoundError):
+        workflow.get_observations()
+
+
+@requires_data()
+def test_get_observations_all():
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.datastore = "$GAMMAPY_DATA/cta-1dc/index/gps/"
+    workflow.get_observations()
+    assert len(workflow.observations) == 4
+
+
+@requires_data()
+def test_get_observations_obs_ids():
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.datastore = "$GAMMAPY_DATA/cta-1dc/index/gps/"
+    workflow.config.observations.obs_ids = ["110380"]
+    workflow.get_observations()
+    assert len(workflow.observations) == 1
+
+
+@requires_data()
+def test_get_observations_obs_cone():
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.datastore = "$GAMMAPY_DATA/hess-dl3-dr1"
+    workflow.config.observations.obs_cone = {
+        "frame": "icrs",
+        "lon": "83d",
+        "lat": "22d",
+        "radius": "5d",
+    }
+    workflow.get_observations()
+    assert len(workflow.observations) == 4
+
+
+@requires_data()
+def test_get_observations_obs_file(tmp_path):
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.get_observations()
+    filename = tmp_path / "obs_ids.txt"
+    filename.write_text("20136\n47829\n")
+    workflow.config.observations.obs_file = filename
+    workflow.get_observations()
+    assert len(workflow.observations) == 2
+
+
+@requires_data()
+def test_get_observations_obs_time(tmp_path):
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.obs_time = {
+        "start": "2004-03-26",
+        "stop": "2004-05-26",
+    }
+    workflow.get_observations()
+    assert len(workflow.observations) == 40
+    workflow.config.observations.obs_ids = [0]
+    with pytest.raises(KeyError):
+        workflow.get_observations()
+
+
+@requires_data()
+def test_get_observations_missing_irf():
+    config = WorkflowConfig()
+    workflow = Workflow(config)
+    workflow.config.observations.datastore = "$GAMMAPY_DATA/joint-crab/dl3/magic/"
+    workflow.config.observations.obs_ids = ["05029748"]
+    workflow.config.observations.required_irf = ["aeff", "edisp"]
+    workflow.get_observations()
+    assert len(workflow.observations) == 1
+
+
+@requires_data()
+def test_set_models():
+    config = get_example_config("3d")
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    models_str = Path(MODEL_FILE).read_text()
+    workflow.set_models(models=models_str)
+    assert isinstance(workflow.models, DatasetModels)
+    assert len(workflow.models) == 2
+    assert workflow.models.names == ["source", "stacked-bkg"]
+    with pytest.raises(TypeError):
+        workflow.set_models(0)
+
+    new_source = workflow.models["source"].copy(name="source2")
+    workflow.set_models(models=[new_source], extend=False)
+    assert len(workflow.models) == 2
+    assert workflow.models.names == ["source2", "stacked-bkg"]
+
+
+@requires_data()
+def test_workflow_1d():
+    cfg = """
+    observations:
+        datastore: $GAMMAPY_DATA/hess-dl3-dr1
+        obs_ids: [23523, 23526]
+        obs_time: {
+            start: [J2004.92654346, J2004.92658453, J2004.92663655],
+            stop: [J2004.92658453, J2004.92663655, J2004.92670773]
+        }
+    datasets:
+        type: 1d
+        background:
+            method: reflected
+        geom:
+            axes:
+                energy_true: {min: 0.01 TeV, max: 300 TeV, nbins: 109}
+        on_region: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 0.11 deg}
+        safe_mask:
+            methods: [aeff-default, edisp-bias]
+            parameters: {bias_percent: 10.0}
+        containment_correction: false
+    flux_points:
+        energy: {min: 1 TeV, max: 50 TeV, nbins: 4}
+    light_curve:
+        energy_edges: {min: 1 TeV, max: 50 TeV, nbins: 1}
+        time_intervals: {
+            start: [J2004.92654346, J2004.92658453, J2004.92663655],
+            stop: [J2004.92658453, J2004.92663655, J2004.92670773]
+        }
+    """
+    config = get_example_config("1d")
+    workflow = Workflow(config)
+    workflow.update_config(cfg)
+    workflow.get_observations()
+    workflow.get_datasets()
+    workflow.read_models(MODEL_FILE_1D)
+    workflow.run_fit()
+    workflow.get_flux_points()
+    workflow.get_light_curve()
+
+    assert len(workflow.datasets) == 3
+    table = workflow.flux_points.data.to_table(sed_type="dnde")
+
+    assert len(table) == 4
+    dnde = table["dnde"].quantity
+    assert dnde.unit == "cm-2 s-1 TeV-1"
+
+    assert_allclose(dnde[0].value, 8.116854e-12, rtol=1e-2)
+    assert_allclose(dnde[2].value, 3.444475e-14, rtol=1e-2)
+
+    axis = workflow.light_curve.geom.axes["time"]
+    assert axis.nbin == 3
+    assert_allclose(axis.time_min.mjd, [53343.92, 53343.935, 53343.954])
+
+    flux = workflow.light_curve.flux.data[:, :, 0, 0]
+    assert_allclose(flux, [[1.688954e-11], [2.347870e-11], [1.604152e-11]], rtol=1e-4)
+
+
+@requires_data()
+def test_geom_workflow_1d():
+    cfg = """
+    observations:
+        datastore: $GAMMAPY_DATA/hess-dl3-dr1
+        obs_ids: [23523]
+    datasets:
+        type: 1d
+        background:
+            method: reflected
+        on_region: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 0.11 deg}
+        geom:
+            axes:
+                energy: {min: 0.1 TeV, max: 30 TeV, nbins: 20}
+                energy_true: {min: 0.03 TeV, max: 100 TeV, nbins: 50}
+        containment_correction: false
+    flux_points:
+        energy: {min: 1 TeV, max: 50 TeV, nbins: 4}
+    """
+    config = get_example_config("1d")
+    workflow = Workflow(config)
+    workflow.update_config(cfg)
+    workflow.get_observations()
+    workflow.get_datasets()
+
+    assert len(workflow.datasets) == 1
+
+    axis = workflow.datasets[0].exposure.geom.axes["energy_true"]
+    assert axis.nbin == 50
+    assert_allclose(axis.edges[0].to_value("TeV"), 0.03)
+    assert_allclose(axis.edges[-1].to_value("TeV"), 100)
+
+
+@requires_data()
+def test_exclusion_region(tmp_path):
+    config = get_example_config("1d")
+    workflow = Workflow(config)
+    region = CircleSkyRegion(center=SkyCoord("85d 23d"), radius=1 * u.deg)
+    geom = WcsGeom.create(npix=(150, 150), binsz=0.05, skydir=SkyCoord("83d 22d"))
+    exclusion_mask = ~geom.region_mask([region])
+
+    filename = tmp_path / "exclusion.fits"
+    exclusion_mask.write(filename)
+    config.datasets.background.method = "reflected"
+    config.datasets.background.exclusion = filename
+    workflow.get_observations()
+    workflow.get_datasets()
+    assert len(workflow.datasets) == 2
+
+    config = get_example_config("3d")
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    geom = workflow.datasets[0]._geom
+    exclusion_mask = ~geom.region_mask([region])
+    filename = tmp_path / "exclusion3d.fits"
+    exclusion_mask.write(filename)
+    config.datasets.background.exclusion = filename
+    workflow.get_datasets()
+    assert len(workflow.datasets) == 1
+
+
+@requires_data()
+def test_workflow_1d_stacked_no_fit_range():
+    cfg = """
+    observations:
+        datastore: $GAMMAPY_DATA/hess-dl3-dr1
+        obs_cone: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 5 deg}
+        obs_ids: [23592, 23559]
+
+    datasets:
+        type: 1d
+        stack: false
+        geom:
+            axes:
+                energy: {min: 0.01 TeV, max: 100 TeV, nbins: 73}
+                energy_true: {min: 0.03 TeV, max: 100 TeV, nbins: 50}
+        on_region: {frame: icrs, lon: 83.633 deg, lat: 22.014 deg, radius: 0.1 deg}
+        containment_correction: true
+        background:
+            method: reflected
+    """
+    config = WorkflowConfig.from_yaml(cfg)
+    workflow = Workflow(config)
+    workflow.update_config(cfg)
+    workflow.config.datasets.stack = True
+    workflow.get_observations()
+    workflow.get_datasets()
+    workflow.read_models(MODEL_FILE_1D)
+    workflow.run_fit()
+    with pytest.raises(ValueError):
+        workflow.get_excess_map()
+
+    assert len(workflow.datasets) == 1
+    assert_allclose(workflow.datasets["stacked"].counts.data.sum(), 184)
+    pars = workflow.models.parameters
+    assert_allclose(workflow.datasets[0].mask_fit.data, True)
+
+    assert_allclose(pars["index"].value, 2.76913, rtol=1e-2)
+    assert_allclose(pars["amplitude"].value, 5.479729e-11, rtol=1e-2)
+
+
+@requires_data()
+def test_workflow_ring_background():
+    config = get_example_config("3d")
+    config.datasets.background.method = "ring"
+    config.datasets.background.parameters = {"r_in": "0.7 deg", "width": "0.7 deg"}
+    config.datasets.geom.axes.energy.nbins = 1
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    workflow.get_excess_map()
+    assert isinstance(workflow.datasets[0], MapDataset)
+    assert_allclose(
+        workflow.datasets[0].npred_background().data[0, 10, 10], 0.091799, rtol=1e-2
+    )
+    assert isinstance(workflow.excess_map["sqrt_ts"], WcsNDMap)
+    assert_allclose(workflow.excess_map.npred_excess.data[0, 62, 62], 134.12389)
+
+
+@requires_data()
+def test_workflow_ring_3d():
+    config = get_example_config("3d")
+    config.datasets.background.method = "ring"
+    config.datasets.background.parameters = {"r_in": "0.7 deg", "width": "0.7 deg"}
+    workflow = Workflow(config)
+    workflow.get_observations()
+    with pytest.raises(ValueError):
+        workflow.get_datasets()
+
+
+@requires_data()
+def test_workflow_no_bkg_1d(caplog):
+    config = get_example_config("1d")
+    workflow = Workflow(config)
+    with caplog.at_level(logging.WARNING):
+        workflow.get_observations()
+        workflow.get_datasets()
+        assert not isinstance(workflow.datasets[0], SpectrumDatasetOnOff)
+        assert "No background maker set. Check configuration." in [
+            _.message for _ in caplog.records
+        ]
+
+
+@requires_data()
+def test_workflow_no_bkg_3d(caplog):
+    config = get_example_config("3d")
+    config.datasets.background.method = None
+    workflow = Workflow(config)
+    with caplog.at_level(logging.WARNING):
+        workflow.get_observations()
+        workflow.get_datasets()
+        assert isinstance(workflow.datasets[0], MapDataset)
+        assert "No background maker set. Check configuration." in [
+            _.message for _ in caplog.records
+        ]
+
+
+@requires_data()
+def test_workflow_3d():
+    config = get_example_config("3d")
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    workflow.read_models(MODEL_FILE)
+    workflow.datasets["stacked"].background_model.spectral_model.tilt.frozen = False
+    workflow.run_fit()
+    workflow.get_flux_points()
+
+    assert len(workflow.datasets) == 1
+    assert len(workflow.models.parameters) == 8
+    res = workflow.models.parameters
+    assert res["amplitude"].unit == "cm-2 s-1 TeV-1"
+
+    table = workflow.flux_points.data.to_table(sed_type="dnde")
+    assert len(table) == 2
+    dnde = table["dnde"].quantity
+
+    assert_allclose(dnde[0].value, 1.2722e-11, rtol=1e-2)
+    assert_allclose(dnde[-1].value, 4.054128e-13, rtol=1e-2)
+    assert_allclose(res["index"].value, 2.772814, rtol=1e-2)
+    assert_allclose(res["tilt"].value, -0.133436, rtol=1e-2)
+
+
+@requires_data()
+def test_workflow_3d_joint_datasets():
+    config = get_example_config("3d")
+    config.datasets.stack = False
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    assert len(workflow.datasets) == 2
+
+    assert_allclose(
+        workflow.datasets[0].background_model.spectral_model.norm.value,
+        1.031743694988066,
+        rtol=1e-6,
+    )
+    assert_allclose(
+        workflow.datasets[0].background_model.spectral_model.tilt.value,
+        0.0,
+        rtol=1e-6,
+    )
+    assert_allclose(
+        workflow.datasets[1].background_model.spectral_model.norm.value,
+        0.9776349021876344,
+        rtol=1e-6,
+    )
+
+
+@requires_data()
+def test_usage_errors():
+    config = get_example_config("1d")
+    workflow = Workflow(config)
+    with pytest.raises(RuntimeError):
+        workflow.get_datasets()
+    with pytest.raises(RuntimeError):
+        workflow.read_datasets()
+    with pytest.raises(RuntimeError):
+        workflow.write_datasets()
+    with pytest.raises(TypeError):
+        workflow.read_models()
+    with pytest.raises(RuntimeError):
+        workflow.write_models()
+    with pytest.raises(RuntimeError):
+        workflow.run_fit()
+    with pytest.raises(RuntimeError):
+        workflow.get_flux_points()
+    with pytest.raises(ValidationError):
+        workflow.config.datasets.type = "None"
+
+
+@requires_data()
+def test_datasets_io(tmpdir):
+    config = get_example_config("3d")
+
+    workflow = Workflow(config)
+    workflow.get_observations()
+    workflow.get_datasets()
+    models_str = Path(MODEL_FILE).read_text()
+    workflow.models = models_str
+
+    config.general.datasets_file = tmpdir / "datasets.yaml"
+    config.general.models_file = tmpdir / "models.yaml"
+    workflow.write_datasets()
+    workflow = Workflow(config)
+    workflow.read_datasets()
+    assert len(workflow.datasets.models) == 2
+    assert workflow.models.names == ["source", "stacked-bkg"]
+
+    workflow.models[0].parameters["index"].value = 3
+    workflow.write_models()
+    workflow = Workflow(config)
+    workflow.read_datasets()
+    assert len(workflow.datasets.models) == 2
+    assert workflow.models.names == ["source", "stacked-bkg"]
+    assert workflow.models[0].parameters["index"].value == 3