From bf117650f6fdcce1dfaa27f1166e6c5d5f55df2d Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Wed, 13 Nov 2024 13:55:55 +0100
Subject: [PATCH 01/10] Start gammapy interface module

---
 src/iact_estimator/gammapy_interface.py | 208 ++++++++++++++++++++++++
 1 file changed, 208 insertions(+)
 create mode 100644 src/iact_estimator/gammapy_interface.py

diff --git a/src/iact_estimator/gammapy_interface.py b/src/iact_estimator/gammapy_interface.py
new file mode 100644
index 0000000..a611d39
--- /dev/null
+++ b/src/iact_estimator/gammapy_interface.py
@@ -0,0 +1,208 @@
+"""Interface to gammapy related to computation."""
+
+import logging
+
+import astropy.units as u
+from astropy.coordinates import Angle
+from astropy.table import Table
+import numpy as np
+
+from gammapy.data import DataStore, FixedPointingInfo, Observation, Observations
+from gammapy.datasets import Datasets, SpectrumDatasetOnOff
+
+from gammapy.maps import MapAxis, RegionGeom
+
+from gammapy.estimators import FluxPointsEstimator
+from gammapy.stats import wstat
+
+from regions import PointSkyRegion, CircleSkyRegion
+
+from .core import log_and_raise
+
+logger = logging.getLogger(__name__)
+
+
+def load_real_observations(data_store_path, required_irfs):
+    data_store = DataStore.from_dir(data_store_path)
+    observations = data_store.get_observations(required_irf=required_irfs)
+
+    return observations
+
+
+def create_simulated_observation(
+    pointing_coordinates,
+    observer,
+    observation_livetime,
+    irfs,
+):
+    simulated_observation = Observation.create(
+        pointing=FixedPointingInfo(
+            fixed_icrs=pointing_coordinates.icrs,
+        ),
+        livetime=observation_livetime,
+        irfs=irfs,
+        location=observer.location,
+    )
+    observations = Observations([simulated_observation])
+    return observations
+
+
+def has_radmax_2d(observation):
+    rad_max = observation.rad_max
+
+    return True if (rad_max is not None and rad_max.has_offset_axis) else False
+
+
+def load_energy_axis_from_config(config, which):
+    axis_config = config["datasets"]["geom"]["axes"][which]
+    axis = MapAxis.from_energy_bounds(
+        u.Quantity(axis_config["min"]),
+        u.Quantity(axis_config["max"]),
+        u.Quantity(axis_config["nbins"]),
+        per_decade=axis_config["per_decade"],
+        name=which,
+    )
+
+    return axis
+
+
+def define_on_region_geometry(
+    target_source_coordinates,
+    observation,
+    offset,
+    reconstructed_energy_axis,
+    on_region_radius=None,
+):
+    if has_radmax_2d(observation):
+        on_region = PointSkyRegion(target_source_coordinates)
+    else:
+        if not on_region_radius:
+            log_and_raise(
+                logger,
+                "Your IRFs do not have angular energy dependence and you did not define a fixed ON region radius.",
+                exc_type=ValueError,
+            )
+
+        on_region_radius = Angle(on_region_radius)
+
+        center = target_source_coordinates.directional_offset_by(
+            position_angle=0 * u.deg, separation=Angle(offset)
+        )
+
+        on_region = CircleSkyRegion(center=center, radius=on_region_radius)
+
+    on_region_geometry = RegionGeom.create(
+        region=on_region, axes=[reconstructed_energy_axis]
+    )
+
+    return on_region_geometry
+
+
+def fake_onoff_from_real_observations(
+    observations,
+    empty_spectrum_dataset,
+    observation_livetime,
+    spectrum_dataset_maker,
+    background_maker,
+    sky_models,
+):
+    """Produce ON/OFF spectrum datasets from real observations.
+
+    This is the use case for experiments like MAGIC which provide DL3
+    files with run-wise IRFs and real OFF counts in place of a background model."""
+
+    logger.debug("Faking ON / OFF spectrum dataset from real observations.")
+
+    fake_spectrum_datasets_on_off = Datasets()
+
+    for observation in observations:
+        spectrum_dataset = spectrum_dataset_maker.run(
+            empty_spectrum_dataset.copy(name=str(observation.obs_id)), observation
+        )
+
+        spectrum_dataset.models = sky_models
+
+        spectrum_dataset.exposure *= (
+            observation_livetime / observation.observation_live_time_duration.to("h")
+        )
+        spectrum_dataset.fake()
+
+        # Make a real ON/OFF dataset because we need the real OFF counts
+        # since these IRFs do not have simulated background
+        fake_spectrum_dataset_on_off = background_maker.run(
+            spectrum_dataset, observation
+        )
+        fake_spectrum_dataset_on_off.models = sky_models
+        # then fake it using its OFF counts
+        fake_spectrum_dataset_on_off.fake(
+            npred_background=fake_spectrum_dataset_on_off.npred_background()
+        )
+        fake_spectrum_datasets_on_off.append(fake_spectrum_dataset_on_off)
+
+    return fake_spectrum_datasets_on_off
+
+
+def fake_onoff_from_fake_observation(
+    simulated_observation,
+    empty_spectrum_dataset,
+    spectrum_dataset_maker,
+    sky_models,
+    n_off_regions,
+    target_source,
+):
+    """Produce ON/OFF spectrum datasets from a simulated observation.
+
+    This is the use case for CTA which provide IRFs with a background model
+    but no real events (for now!).
+    """
+
+    spectrum_dataset = spectrum_dataset_maker.run(
+        empty_spectrum_dataset, simulated_observation
+    )
+
+    sky_models[target_source.name].spatial_model = None
+    spectrum_dataset.models = sky_models
+    spectrum_dataset.fake()
+
+    spectrum_dataset_on_off = SpectrumDatasetOnOff.from_spectrum_dataset(
+        dataset=spectrum_dataset, acceptance=1, acceptance_off=n_off_regions
+    )
+    spectrum_dataset_on_off.fake(npred_background=spectrum_dataset.npred_background())
+
+    return Datasets(spectrum_dataset_on_off)
+
+
+def estimate_sed(target_source, energy_axis, spectrum_datasets_on_off, **kwargs):
+    flux_points_estimator = FluxPointsEstimator(
+        source=target_source.name,
+        energy_edges=energy_axis.edges,
+        selection_optional="all",
+        **kwargs,
+    )
+
+    flux_points = flux_points_estimator.run(spectrum_datasets_on_off)
+
+    return flux_points
+
+
+def get_wstat_table(spectrum_dataset_on_off):
+    table = Table()
+    table["mu_sig"] = spectrum_dataset_on_off.npred_signal().data.flatten()
+    table["n_on"] = spectrum_dataset_on_off.counts.data.flatten()
+    table["n_off"] = spectrum_dataset_on_off.counts_off.data.flatten()
+    table["alpha"] = spectrum_dataset_on_off.alpha.data.flatten()
+    table["wstat"] = wstat(
+        n_on=table["n_on"],
+        n_off=table["n_off"],
+        alpha=table["alpha"],
+        mu_sig=table["mu_sig"],
+    )
+    table["li_ma_sigma"] = np.sqrt(
+        wstat(
+            n_on=table["n_on"],
+            n_off=table["n_off"],
+            alpha=table["alpha"],
+            mu_sig=table["mu_sig"],
+        )
+    )
+    return table

From 7ba8a5b2e48b50dba74dd79844fa6fcc795226fa Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Wed, 13 Nov 2024 13:56:25 +0100
Subject: [PATCH 02/10] Add gammapy interface configuration

---
 src/iact_estimator/resources/config.yml | 28 ++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/src/iact_estimator/resources/config.yml b/src/iact_estimator/resources/config.yml
index 45190c4..684a501 100644
--- a/src/iact_estimator/resources/config.yml
+++ b/src/iact_estimator/resources/config.yml
@@ -1,4 +1,6 @@
 assumed_model:
+  # NOTE: this applies to the old version scripts
+  # for gammapy use models.yml
   # see https://docs.gammapy.org/1.1/user-guide/model-gallery/index.html#spectral-models
   name: gammapy.modeling.models.LogParabolaSpectralModel
   parameters:
@@ -44,9 +46,33 @@ plotting_options:
   file_format: "pdf"
   merge_horizon_profiles: True
 
+gammapy_interface:
+  # this section follows loosely gammapy's HLI but does NOT use it
+  # see https://github.com/gammapy/gammapy/issues/4169
+  # add any key that you want to overwrite based on gammapy low-level interface
+  input: # path to IRF file or path to real DL3 data compatible with GADF / `gammapy.data.DataStore`
+  datasets:
+    type: 1d
+    geom:
+        axes:
+            energy: {min: 70 GeV, max: 20 TeV, nbins: 5, per_decade: True}
+            energy_true: {min: 10 GeV, max: 100 TeV, nbins: 10, per_decade: True}
+    on_region:
+      frame: icrs,
+      radius: null # null for energy dependence
+      offset: "0.4 deg"
+    containment_correction: false
+  estimators:
+    flux_points:
+  plotting_options:
+    sed:
+      ts_profile: True
+      energy_flux_unit: "TeV cm-2 s-1"
+
+
 use_seaborn: True
 seaborn_options:
-  context: talk
+  context: paper
   style: whitegrid
   palette: viridis
   font: sans-serif

From 19286c5bf6077fb366b2330b1a3d0350e923122b Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Wed, 13 Nov 2024 13:56:46 +0100
Subject: [PATCH 03/10] Add models configuration file

---
 src/iact_estimator/resources/models.yml | 29 +++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 src/iact_estimator/resources/models.yml

diff --git a/src/iact_estimator/resources/models.yml b/src/iact_estimator/resources/models.yml
new file mode 100644
index 0000000..6b2e0c8
--- /dev/null
+++ b/src/iact_estimator/resources/models.yml
@@ -0,0 +1,29 @@
+# NOTE: This file is only for use with the gammapy interface
+components:
+- name: Crab Nebula
+  type: SkyModel
+  spatial:
+    type: PointSpatialModel
+    frame: icrs
+    parameters:
+    - name: lon_0
+      value: 83.63
+      unit: deg
+    - name: lat_0
+      value: 22.014
+      unit: deg
+  spectral:
+    type: LogParabolaSpectralModel
+    parameters:
+    - name: amplitude
+      value: 3.39e-11
+      unit: cm-2 s-1 TeV-1
+    - name: alpha
+      value: 2.51
+      unit: ''
+    - name: beta
+      value: 0.09
+      unit: ''
+    - name: reference
+      value: 1.0
+      unit: TeV

From 7ef10dfc651a0095986078deb1fb1f40bed78709 Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Wed, 13 Nov 2024 13:57:16 +0100
Subject: [PATCH 04/10] Add core.log_and_raise() function

---
 src/iact_estimator/core.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/iact_estimator/core.py b/src/iact_estimator/core.py
index c42b5ff..a170a3a 100644
--- a/src/iact_estimator/core.py
+++ b/src/iact_estimator/core.py
@@ -85,6 +85,11 @@ def setup_logging(log_level, source_name):
     return logger
 
 
+def log_and_raise(logger, msg, exc_type=ValueError):
+    logger.error(msg)
+    raise exc_type(msg)
+
+
 def get_horizon_stereo_profile(M1_data, M2_data):
     """Merge the horizon profiles into a stereo one.
 

From 8bd0b7038ebecfd5e27dff6c8de8ef0d087f7fce Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Wed, 13 Nov 2024 13:58:56 +0100
Subject: [PATCH 05/10] Add utility functions for seaborn figures

---
 src/iact_estimator/plots.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)

diff --git a/src/iact_estimator/plots.py b/src/iact_estimator/plots.py
index 8af48a9..36f5523 100644
--- a/src/iact_estimator/plots.py
+++ b/src/iact_estimator/plots.py
@@ -9,6 +9,7 @@
 from astroplan.plots import plot_sky_24hr, plot_altitude
 import matplotlib.pyplot as plt
 import numpy as np
+import seaborn as sns
 
 from .core import observed_flux, get_horizon_stereo_profile
 from .spectral import crab_nebula_spectrum
@@ -26,6 +27,34 @@
 logger = logging.getLogger(__name__)
 
 
+def get_default_seaborn_font_scaling(context):
+    context_params = sns.plotting_context(context)
+    font_scale = (
+        context_params["font.size"] / sns.plotting_context("paper")["font.size"]
+    )
+    return font_scale
+
+
+def set_context_with_scaled_figsize(
+    context="notebook", base_figure_size=plt.rcParams["figure.figsize"]
+):
+    """
+    Sets Seaborn context and adjusts figure size based on context's default scaling factor.
+
+    Parameters:
+        context (str): Seaborn context ('paper', 'notebook', 'talk', 'poster')
+    """
+
+    context_params = sns.plotting_context(context)
+    font_scale = (
+        context_params["font.size"] / sns.plotting_context("paper")["font.size"]
+    )
+
+    figure_size = (base_figure_size[0] * font_scale, base_figure_size[1] * font_scale)
+
+    return figure_size
+
+
 def plot_spectrum(
     config,
     energy_bounds,

From 61c378d838ae80665b2c06b90878e79b398e6feb Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Thu, 14 Nov 2024 10:41:48 +0100
Subject: [PATCH 06/10] Update plots module

- add plot_gammapy_sed
add plot_irf_source_pointing

Fix quotes formatting
---
 src/iact_estimator/plots.py | 139 +++++++++++++++++++++++++++++++++++-
 1 file changed, 138 insertions(+), 1 deletion(-)

diff --git a/src/iact_estimator/plots.py b/src/iact_estimator/plots.py
index 36f5523..0f35568 100644
--- a/src/iact_estimator/plots.py
+++ b/src/iact_estimator/plots.py
@@ -4,13 +4,16 @@
 from pathlib import Path
 
 import astropy.units as u
+from astropy.coordinates import AltAz
 from astropy.visualization import quantity_support
-from astroplan import FixedTarget
+from astroplan import FixedTarget, time_grid_from_range
 from astroplan.plots import plot_sky_24hr, plot_altitude
 import matplotlib.pyplot as plt
 import numpy as np
 import seaborn as sns
 
+from gammapy.maps import MapAxis
+
 from .core import observed_flux, get_horizon_stereo_profile
 from .spectral import crab_nebula_spectrum
 from iact_estimator import HORIZON_PROFILE_M1, HORIZON_PROFILE_M2
@@ -22,6 +25,8 @@
     "plot_altitude_airmass",
     "plot_exposure",
     "plot_rates",
+    "plot_gammapy_sed",
+    "plot_irf_source_pointing",
 ]
 
 logger = logging.getLogger(__name__)
@@ -404,3 +409,135 @@ def plot_rates(performance_data, title=None, ax=None):
         ax.set_xscale("log")
         ax.set_yscale("log")
     return ax
+
+
+def plot_gammapy_sed(
+    observation_livetime,
+    flux_points,
+    assumed_spectral_model,
+    plot_ts_profiles=True,
+    energy_flux_unit="TeV cm-2 s-1",
+    source_name=None,
+    save=True,
+    output_path=None,
+    **kwargs,
+):
+    """Plot an SED with gammapy shoowing the significance per bin.
+
+    Optionally plot TS profiles.
+
+    Parameters
+    ==========
+    observation_livetime : `astropy.units.Quantity`
+    flux_points : `~gammapy.estimators.FluxPoints`
+    assumed_spectral_model : `~gammapy.modeling.models.SpectralModel`
+    plot_ts_profiles=True,
+    energy_flux_unit : str, default="TeV cm-2 s-1"
+    source_name : str, default=None
+    savefig : bool, default=True
+    output_path : str or `pathlib.Path`, default=None
+    kwargs :
+        Keyword arguments for `~matplotlib.axis.Axis.grid()`
+        and `~matplotlib.figure.Figure.savefig()`.
+    """
+
+    grid_kwargs = kwargs.get("grid", {})
+    savefig_kwargs = kwargs.get("savefig", {})
+    seaborn_kwargs = kwargs.get("seaborn", {})
+
+    figsize = set_context_with_scaled_figsize(
+        context=seaborn_kwargs["context"],
+        base_figure_size=plt.rcParams["figure.figsize"],
+    )
+    fig, ax = plt.subplots(figsize=figsize, layout="constrained")
+
+    sed_type = "e2dnde"
+    energy_flux_unit = u.Unit(energy_flux_unit)
+
+    ax.yaxis.units = energy_flux_unit
+
+    flux_points_table = flux_points.to_table(
+        sed_type="e2dnde", format="gadf-sed", formatted=True
+    )
+
+    energy_axis = MapAxis.from_energy_edges(
+        np.concatenate(
+            (
+                flux_points_table["e_min"].data,
+                np.array([flux_points_table["e_max"][-1]]),
+            )
+        )
+        * flux_points_table["e_min"].unit
+    )
+
+    assumed_spectral_model.plot(
+        [energy_axis.edges.min().value, energy_axis.edges.max().value]
+        * energy_axis.edges.unit,
+        ax=ax,
+        sed_type=sed_type,
+        label="Assumed spectral model",
+        color="green",
+    )
+
+    flux_points.plot(
+        ax=ax,
+        sed_type=sed_type,
+        color="darkorange",
+        label=f"Estimated observation ({observation_livetime})",
+    )
+
+    if plot_ts_profiles:
+        flux_points.plot_ts_profiles(ax=ax, sed_type=sed_type, cmap="Blues")
+
+    for flux_point in flux_points_table:
+        ax.annotate(
+            rf"{flux_point['sqrt_ts']:.1f} $\sigma$",
+            xy=(
+                flux_point["e_min"],
+                (
+                    flux_point["e2dnde"] * flux_point.table["e2dnde"].unit
+                    + 3 * flux_point["e2dnde_err"] * flux_point.table["e2dnde_err"].unit
+                ).to_value(energy_flux_unit),
+            ),
+            rotation=45,
+        )
+
+    ax.grid(**grid_kwargs)
+
+    ax.set_ylabel(
+        rf"$E^{2}\dfrac{{dN}}{{dE}}\quad$({energy_flux_unit.to_string('latex',fraction=False)})"
+    )
+    ax.legend()
+
+    if save:
+        output_path = output_path if output_path is not None else Path.cwd()
+        fig.savefig(
+            output_path / f"{source_name or ''}_gammapy_sed.{savefig_kwargs['format']}",
+            **savefig_kwargs,
+        )
+
+
+def plot_irf_source_pointing(target_source, observation, ax=None, **kwargs):
+    ax = plt.gca() if ax else None
+
+    observation_times = time_grid_from_range(
+        [observation.tstart, observation.tstop], time_resolution=1 * u.h
+    )
+
+    altaz_pointings = observation.get_pointing_altaz(observation_times)
+    target_source_altaz = target_source.coord.transform_to(
+        AltAz(
+            obstime=observation_times, location=observation.observatory_earth_location
+        )
+    )
+
+    ax.plot(altaz_pointings.az, altaz_pointings.alt, label="IRFs pointing")
+    ax.plot(target_source_altaz.az, target_source_altaz.alt, label="Target source")
+
+    ax.set_ylabel("Zenith angle (deg)")
+    ax.set_xlabel("Azimuth angle (deg)")
+    ax.grid()
+
+    ax.legend()
+
+    return ax

From 0906e05e8e7a2ca7e1fc906640dabd1560f0bb5b Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Thu, 14 Nov 2024 10:42:28 +0100
Subject: [PATCH 07/10] Add gammapy interface execution to main

---
 src/iact_estimator/scripts/main.py | 184 ++++++++++++++++++++++++++++-
 1 file changed, 179 insertions(+), 5 deletions(-)

diff --git a/src/iact_estimator/scripts/main.py b/src/iact_estimator/scripts/main.py
index e9cb4a5..d2d65df 100644
--- a/src/iact_estimator/scripts/main.py
+++ b/src/iact_estimator/scripts/main.py
@@ -9,6 +9,15 @@
 from astropy.time import Time
 import astropy.units as u
 from astropy.visualization import quantity_support
+from gammapy.irf import load_irf_dict_from_file
+
+from gammapy.datasets import SpectrumDataset
+from gammapy.makers import (
+    ReflectedRegionsBackgroundMaker,
+    SpectrumDatasetMaker,
+    WobbleRegionsFinder,
+)
+from gammapy.modeling.models import Models
 import matplotlib.pyplot as plt
 
 from .. import __version__
@@ -21,7 +30,24 @@
     source_detection,
     calculate,
 )
-from ..plots import plot_spectrum, plot_sed, plot_transit, plot_altitude_airmass
+from ..gammapy_interface import (
+    load_real_observations,
+    create_simulated_observation,
+    has_radmax_2d,
+    define_on_region_geometry,
+    estimate_sed,
+    fake_onoff_from_fake_observation,
+    fake_onoff_from_real_observations,
+    load_energy_axis_from_config,
+    get_wstat_table,
+)
+from ..plots import (
+    plot_spectrum,
+    plot_sed,
+    plot_transit,
+    plot_altitude_airmass,
+    plot_gammapy_sed,
+)
 from .. import RESOURCES_PATH
 
 parser = argparse.ArgumentParser()
@@ -32,7 +58,7 @@
 
 parser_config = subparsers.add_parser(
     "config",
-    help="Copy the default configuration file to your new project directory.",
+    help="Copy the configuration files to your new project directory.",
 )
 
 parser_config.add_argument(
@@ -47,9 +73,19 @@
     help="Launch the estimation of observability and physical properties of the target source.",
 )
 
+parser_run.add_argument(
+    "--use-gammapy", action="store_true", help="Use also the new gammapy interface."
+)
+
 parser_run.add_argument(
     "--config", required=True, type=str, help="Path to configuration file."
 )
+parser_run.add_argument(
+    "--models",
+    required=False,
+    type=str,
+    help="Path to gammapy model file (models.yml).",
+)
 parser_run.add_argument(
     "--source-name",
     default="test_source",
@@ -88,6 +124,7 @@ def main():
             Path.cwd() if args_config.to is None else Path(args_config.to)
         )
         shutil.copy(RESOURCES_PATH / "config.yml", config_output_path)
+        shutil.copy(RESOURCES_PATH / "models.yml", config_output_path)
 
         parser.exit(0, message="You're halfway there! Have fun!\n")
 
@@ -154,9 +191,9 @@ def main():
             energy_bins, gamma_rate, background_rate, config, assumed_spectrum
         )
 
-        combined_significance = source_detection(
-            sigmas, u.Quantity(config["observation_time"])
-        )
+        simulated_observation_livetime = u.Quantity(config["observation_time"])
+
+        combined_significance = source_detection(sigmas, simulated_observation_livetime)
 
         with quantity_support():
             plot_sed(
@@ -181,6 +218,143 @@ def main():
 
         crab = FixedTarget.from_name("Crab")
 
+        target_source_coordinates = target_source.coord
+
+        if args.use_gammapy:
+            gammapy_config = config["gammapy_interface"]
+
+            gammapy_analysis_type = gammapy_config["datasets"]["type"]
+
+            if gammapy_analysis_type == "1d":
+                required_irfs = "point-like"
+            else:
+                raise NotImplementedError(
+                    "Only 1D analyses from point-like simulations are available at the moment with the gammapy interface."
+                )
+
+            gammapy_input = Path(gammapy_config["input"])
+
+            sky_models = Models.from_dict(read_yaml(args.models))
+
+            if gammapy_input.is_dir():
+                logger.info(
+                    "Creating data store from real observations at %s", gammapy_input
+                )
+                observations = load_real_observations(gammapy_input, required_irfs)
+            else:
+                logger.info(
+                    "Loading instrument response functions from %s", gammapy_input
+                )
+                irfs = load_irf_dict_from_file(gammapy_input)
+
+                observations = create_simulated_observation(
+                    target_source_coordinates,
+                    observer,
+                    simulated_observation_livetime,
+                    irfs,
+                )
+
+            reco_energy_axis = load_energy_axis_from_config(gammapy_config, "energy")
+            true_energy_axis = load_energy_axis_from_config(
+                gammapy_config, "energy_true"
+            )
+
+            on_region_config = gammapy_config["datasets"]["on_region"]
+            offset = on_region_config["offset"]
+            on_region_radius = on_region_config["radius"] or None
+            on_region_geometry = define_on_region_geometry(
+                target_source_coordinates,
+                observations[0],
+                offset,
+                reco_energy_axis,
+                on_region_radius=on_region_radius,
+            )
+
+            empty_spectrum_dataset = SpectrumDataset.create(
+                geom=on_region_geometry,
+                energy_axis_true=true_energy_axis,
+            )
+
+            logger.info("Angular cut depends on reconstructed energy.")
+
+            is_point_like = has_radmax_2d(observations[0])
+
+            containment_correction = gammapy_config["datasets"][
+                "containment_correction"
+            ]
+            if is_point_like:
+                logger.warning(
+                    "Input observation is point-like, so containment correction will be disabled."
+                )
+                containment_correction = False
+
+            n_off_regions = config["n_off_regions"]
+
+            if len(observations) >= 1 and has_radmax_2d(observations[0]):
+                region_finder = WobbleRegionsFinder(n_off_regions=n_off_regions)
+                background_maker = ReflectedRegionsBackgroundMaker(
+                    region_finder=region_finder
+                )
+
+                spectrum_dataset_maker = SpectrumDatasetMaker(
+                    containment_correction=containment_correction,
+                    selection=["counts", "exposure", "edisp"],
+                )
+
+                empty_spectrum_dataset = SpectrumDataset.create(
+                    geom=on_region_geometry,
+                    energy_axis_true=true_energy_axis,
+                )
+
+                spectrum_datasets_on_off = fake_onoff_from_real_observations(
+                    observations,
+                    empty_spectrum_dataset,
+                    simulated_observation_livetime,
+                    spectrum_dataset_maker,
+                    background_maker,
+                    sky_models,
+                )
+
+            else:
+                spectrum_dataset_maker = SpectrumDatasetMaker(
+                    containment_correction=containment_correction,
+                    selection=["exposure", "edisp", "background"],
+                )
+
+                spectrum_datasets_on_off = fake_onoff_from_fake_observation(
+                    observations[0],
+                    empty_spectrum_dataset,
+                    spectrum_dataset_maker,
+                    sky_models,
+                    n_off_regions,
+                    target_source,
+                )
+
+            spectrum_datasets_on_off_stacked = spectrum_datasets_on_off.stack_reduce()
+            spectrum_datasets_on_off_stacked.models = sky_models
+            wstat_table = get_wstat_table(spectrum_datasets_on_off_stacked)
+            logging.info("Computing WSTAT on stacked spectrum ON/OFF dataset")
+            print(wstat_table)
+
+            flux_points = estimate_sed(
+                target_source,
+                reco_energy_axis,
+                spectrum_datasets_on_off,
+                **gammapy_config["estimators"]["flux_points"] or {},
+            )
+
+            plot_gammapy_sed(
+                simulated_observation_livetime,
+                flux_points,
+                assumed_spectrum,
+                plot_ts_profiles=True,
+                energy_flux_unit="TeV cm-2 s-1",
+                source_name=target_source.name,
+                grid={"which": "both", "axis": "both", "alpha": 0.5, "ls": "dashed"},
+                savefig={"format": "pdf"},
+                seaborn={"context": config["seaborn_options"]["context"]},
+            )
+
         with quantity_support():
             plot_transit(
                 config,

From 569f91bddc5de154fdd500114bda8054b07245eb Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Thu, 14 Nov 2024 13:52:45 +0100
Subject: [PATCH 08/10] Update documentation for initial gammapy interface

---
 docs/source/changes/33.added.rst |  1 +
 docs/source/userguide.rst        | 34 +++++++++++++++++++++++++++++---
 2 files changed, 32 insertions(+), 3 deletions(-)
 create mode 100644 docs/source/changes/33.added.rst

diff --git a/docs/source/changes/33.added.rst b/docs/source/changes/33.added.rst
new file mode 100644
index 0000000..68e841f
--- /dev/null
+++ b/docs/source/changes/33.added.rst
@@ -0,0 +1 @@
+Added gammapy interface with support for 1D spectral analysis (no EBL absorption).
diff --git a/docs/source/userguide.rst b/docs/source/userguide.rst
index 5a189f9..94209ea 100644
--- a/docs/source/userguide.rst
+++ b/docs/source/userguide.rst
@@ -11,13 +11,20 @@ You can use the :ref:`iact-estimator-cfg` command-line tool
 
 .. code-block::
 
-    iact-estimator-cfg --output-path /where/to/save/config/file
+    iact-estimator config --to /where/to/save/config/file
 
-to get the following example configuration file,
+to get the example configuration files:
+
+- a global configuration file
 
 .. literalinclude:: /../../src/iact_estimator/resources/config.yml
     :language: yaml
 
+- a configuration files to define models to be used with the :ref:`gammapy_interface`,`
+
+.. literalinclude:: /../../src/iact_estimator/resources/models.yml
+    :language: yaml
+
 Launch the simulation
 =====================
 
@@ -32,6 +39,27 @@ for a source named "my_source",
 
     iact-estimator --config config.yml --source-name my_source
 
+Gammapy interface
+-----------------
+
+In addition to the legacy scripts (which can still be used
+in case all you have are gamma and background rates as a function of reconstructed energy)
+_iact-estimator_ provides also an interface based on
+`gammapy <https://gammapy.org/>`_.
+
+In order to use it you only need to add the `--use-gammapy` flag
+to the `run` subcommand.
+
+This interface defines all sky models via the `models.yml` file,
+ignoring the `assumed_model` section in the global configuration file.
+
+.. note::
+
+    The configuration schema will soon change,
+    unifying the sky model specifications.
+
+At the moment only the 1D analysis case without EBL absorption is supported.
+
 Output
 ======
 
@@ -45,7 +73,7 @@ bin satisfies the conditions for the detection.
 Plots
 -----
 
-Tha package comes with a small plotting library
+The package comes with a plotting library
 that allows to plot information about the observability of the source
 and its spectral properties as seen by the telescopes.
 

From e0920b37f9ed94f31c1edf1871dd625d098563ab Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Thu, 14 Nov 2024 14:32:08 +0100
Subject: [PATCH 09/10] Fix Point spatial model for 1D analysis to on region
 center

---
 src/iact_estimator/gammapy_interface.py |  2 --
 src/iact_estimator/resources/models.yml | 10 ----------
 src/iact_estimator/scripts/main.py      | 10 ++++++++--
 3 files changed, 8 insertions(+), 14 deletions(-)

diff --git a/src/iact_estimator/gammapy_interface.py b/src/iact_estimator/gammapy_interface.py
index a611d39..cad7e12 100644
--- a/src/iact_estimator/gammapy_interface.py
+++ b/src/iact_estimator/gammapy_interface.py
@@ -148,7 +148,6 @@ def fake_onoff_from_fake_observation(
     spectrum_dataset_maker,
     sky_models,
     n_off_regions,
-    target_source,
 ):
     """Produce ON/OFF spectrum datasets from a simulated observation.
 
@@ -160,7 +159,6 @@ def fake_onoff_from_fake_observation(
         empty_spectrum_dataset, simulated_observation
     )
 
-    sky_models[target_source.name].spatial_model = None
     spectrum_dataset.models = sky_models
     spectrum_dataset.fake()
 
diff --git a/src/iact_estimator/resources/models.yml b/src/iact_estimator/resources/models.yml
index 6b2e0c8..a69ded7 100644
--- a/src/iact_estimator/resources/models.yml
+++ b/src/iact_estimator/resources/models.yml
@@ -2,16 +2,6 @@
 components:
 - name: Crab Nebula
   type: SkyModel
-  spatial:
-    type: PointSpatialModel
-    frame: icrs
-    parameters:
-    - name: lon_0
-      value: 83.63
-      unit: deg
-    - name: lat_0
-      value: 22.014
-      unit: deg
   spectral:
     type: LogParabolaSpectralModel
     parameters:
diff --git a/src/iact_estimator/scripts/main.py b/src/iact_estimator/scripts/main.py
index d2d65df..f538640 100644
--- a/src/iact_estimator/scripts/main.py
+++ b/src/iact_estimator/scripts/main.py
@@ -17,7 +17,7 @@
     SpectrumDatasetMaker,
     WobbleRegionsFinder,
 )
-from gammapy.modeling.models import Models
+from gammapy.modeling.models import Models, PointSpatialModel
 import matplotlib.pyplot as plt
 
 from .. import __version__
@@ -321,13 +321,19 @@ def main():
                     selection=["exposure", "edisp", "background"],
                 )
 
+                if gammapy_config["datasets"]["type"] == "1d":
+                    sky_models[
+                        target_source.name
+                    ].spatial_model = PointSpatialModel.from_position(
+                        on_region_geometry.region.center
+                    )
+
                 spectrum_datasets_on_off = fake_onoff_from_fake_observation(
                     observations[0],
                     empty_spectrum_dataset,
                     spectrum_dataset_maker,
                     sky_models,
                     n_off_regions,
-                    target_source,
                 )
 
             spectrum_datasets_on_off_stacked = spectrum_datasets_on_off.stack_reduce()

From 6e61758de6d4b0abb3d3f194eee6c9edf96ef1e6 Mon Sep 17 00:00:00 2001
From: Michele Peresano <peresano@mpp.mpg.de>
Date: Thu, 14 Nov 2024 14:45:33 +0100
Subject: [PATCH 10/10] Add _gammapy_interface ref

---
 docs/source/userguide.rst | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/source/userguide.rst b/docs/source/userguide.rst
index 94209ea..19e53fd 100644
--- a/docs/source/userguide.rst
+++ b/docs/source/userguide.rst
@@ -20,7 +20,7 @@ to get the example configuration files:
 .. literalinclude:: /../../src/iact_estimator/resources/config.yml
     :language: yaml
 
-- a configuration files to define models to be used with the :ref:`gammapy_interface`,`
+- a configuration files to define models to be used with the :ref:`gammapy_interface`,
 
 .. literalinclude:: /../../src/iact_estimator/resources/models.yml
     :language: yaml
@@ -39,6 +39,8 @@ for a source named "my_source",
 
     iact-estimator --config config.yml --source-name my_source
 
+.. _gammapy_interface:
+
 Gammapy interface
 -----------------