Repository: gammapy/gammapy (python)

Issue:
FluxPointsEstimator fails if no edisp is set
**Gammapy version**
present dev


**Bug description**

Running `FluxPointsEstimator` on a dataset without an `edisp` gives a very confusing error message `ValueError: Can only stack equivalent maps or cutout of the same map.`


**To Reproduce**

Set `edisp=None` for any dataset (with e_true=e_reco) and proceed with flux points estimation in multiple energy bands (normal use case).

**Expected behaviour**
The analysis fails (as expected) during the fitting stage if `energy_true` has a different binning from the `energy` axis with a more understandable error 
```
File ~/Gammapy-dev/gammapy/gammapy/datasets/evaluator.py:245, in MapEvaluator._edisp_diagonal(self)
    242 @lazyproperty
    243 def _edisp_diagonal(self):
    244     return EDispKernel.from_diagonal_response(
--> 245         energy_axis_true=self.edisp.axes["energy_true"],
    246         energy_axis=self.edisp.axes["energy"],
    247     )

AttributeError: 'NoneType' object has no attribute 'axes'
```
However, if the `energy_axis` and the `energy_axis_true` have same binning, the fit proceeds fine and then the very confusing error is raised during flux points estimation.

1. Either the flux estimator (and also the fit) should fail with a decent error message consistently for no edisp OR
2. the diagonal edisp can be adapted with the `energy_axis` here 
https://github.com/gammapy/gammapy/blob/e9b22cdff573768ac2ab1e2a4a18250071ba3b60/gammapy/datasets/evaluator.py#L243
to use the dataset energy axis.


Pull request description:
Fix evaluation if no edisp is set
Resolves #5274.
Better solution than #5308 :)

Patch:
diff --git a/gammapy/datasets/evaluator.py b/gammapy/datasets/evaluator.py
index 160ad4c43..7961da006 100644
--- a/gammapy/datasets/evaluator.py
+++ b/gammapy/datasets/evaluator.py
@@ -72,6 +72,8 @@ class MapEvaluator:
         self.contributes = True
         self.psf_containment = None
 
+        self._geom_reco_axis = None
+
         if evaluation_mode not in {"local", "global"}:
             raise ValueError(f"Invalid evaluation_mode: {evaluation_mode!r}")
 
@@ -117,10 +119,12 @@ class MapEvaluator:
     @property
     def _geom_reco(self):
         if self.edisp is not None:
-            energy_axis = self.edisp.axes["energy"].copy(name="energy")
+            energy_axis = self.edisp.axes["energy"]
+        elif self._geom_reco_axis is not None:
+            energy_axis = self._geom_reco_axis
         else:
-            energy_axis = self.geom.axes["energy_true"].copy(name="energy")
-        geom = self.geom.to_image().to_cube(axes=[energy_axis])
+            energy_axis = self.geom.axes["energy_true"]
+        geom = self.geom.to_image().to_cube(axes=[energy_axis.copy(name="energy")])
         return geom
 
     @property
@@ -191,13 +195,15 @@ class MapEvaluator:
         del self.position
         del self.cutout_width
 
+        self._geom_reco_axis = geom.axes["energy"]
+
         # lookup edisp
+        del self._edisp_diagonal
         if edisp:
             energy_axis = geom.axes["energy"]
             self.edisp = edisp.get_edisp_kernel(
                 position=self.position, energy_axis=energy_axis
             )
-            del self._edisp_diagonal
 
         # lookup psf
         if psf and self.model.spatial_model:
@@ -242,8 +248,8 @@ class MapEvaluator:
     @lazyproperty
     def _edisp_diagonal(self):
         return EDispKernel.from_diagonal_response(
-            energy_axis_true=self.edisp.axes["energy_true"],
-            energy_axis=self.edisp.axes["energy"],
+            energy_axis_true=self.geom.axes["energy_true"],
+            energy_axis=self._geom_reco.axes["energy"],
         )
 
     def update_spatial_oversampling_factor(self, geom):
@@ -384,7 +390,7 @@ class MapEvaluator:
         npred_reco : `~gammapy.maps.Map`
             Predicted counts in reconstructed energy bins.
         """
-        if self.model.apply_irf["edisp"]:
+        if self.model.apply_irf["edisp"] and self.edisp:
             return apply_edisp(npred, self.edisp)
         else:
             if "energy_true" in npred.geom.axes.names:


Test patch:
diff --git a/gammapy/estimators/points/tests/test_sed.py b/gammapy/estimators/points/tests/test_sed.py
index 056669323..9f4f3280f 100644
--- a/gammapy/estimators/points/tests/test_sed.py
+++ b/gammapy/estimators/points/tests/test_sed.py
@@ -430,6 +430,31 @@ def test_run_map_pwl_reoptimize(fpe_map_pwl_reoptimize):
     assert_allclose(actual, [9.788123, 0.486066, 17.603708], rtol=1e-2)
 
 
+def test_run_no_edip(fpe_pwl, tmpdir):
+    datasets, fpe = fpe_pwl
+
+    datasets = datasets.copy()
+
+    datasets[0].models["source"].apply_irf["edisp"] = False
+    fp = fpe.run(datasets)
+    table = fp.to_table()
+    actual = table["norm"].data
+    assert_allclose(actual, [1.081434, 0.91077, 0.922176], rtol=1e-3)
+
+    datasets[0].edisp = None
+
+    fp = fpe.run(datasets)
+    table = fp.to_table()
+    actual = table["norm"].data
+    assert_allclose(actual, [1.081434, 0.91077, 0.922176], rtol=1e-3)
+
+    datasets[0].models["source"].apply_irf["edisp"] = True
+    fp = fpe.run(datasets)
+    table = fp.to_table()
+    actual = table["norm"].data
+    assert_allclose(actual, [1.081434, 0.91077, 0.922176], rtol=1e-3)
+
+
 @requires_dependency("iminuit")
 @requires_data()
 def test_run_template_npred(fpe_map_pwl, tmpdir):