Improve handling of background IRFs with AltAz alignment by slicing t…

…he background evaluation based on the error on the FoV rotation.

Signed-off-by: Gabriel Emery <[email protected]>

gammapy/makers/map.py

@@ -39,6 +39,8 @@ class MapDatasetMaker(Maker):
         Pad one bin in offset for 2d background map.
         This avoids extrapolation at edges and use the nearest value.
         Default is True.
+    fov_rotation_error_limit : `~astropy.units.Quantity`
+        Maximum error on the rotation angle between AltAz and RaDec frames during background evaluation.
 
     Examples
     --------
@@ -107,10 +109,12 @@ def __init__(
         background_oversampling=None,
         background_interp_missing_data=True,
         background_pad_offset=True,
+        fov_rotation_error_limit=0.5 * u.deg,
     ):
         self.background_oversampling = background_oversampling
         self.background_interp_missing_data = background_interp_missing_data
         self.background_pad_offset = background_pad_offset
+        self.fov_rotation_error_limit = fov_rotation_error_limit
         if selection is None:
             selection = self.available_selection
 
@@ -238,6 +242,7 @@ def make_background(self, geom, observation):
             ontime=observation.observation_time_duration,
             bkg=bkg,
             geom=geom,
+            fov_rotation_error_limit=self.fov_rotation_error_limit,
             oversampling=self.background_oversampling,
             use_region_center=use_region_center,
             obstime=observation.tmid,

gammapy/makers/tests/test_utils.py

@@ -190,6 +190,7 @@ def test_map_background_2d(bkg_2d, fixed_pointing_info):
         ontime="42 s",
         bkg=bkg_2d,
         geom=geom,
+        fov_rotation_error_limit=0.5 * u.deg,
     )
 
     assert_allclose(bkg.data[:, 1, 1], [1.869025, 0.186903], rtol=1e-5)
@@ -200,6 +201,7 @@ def test_map_background_2d(bkg_2d, fixed_pointing_info):
         ontime="42 s",
         bkg=bkg_2d,
         geom=geom,
+        fov_rotation_error_limit=0.5 * u.deg,
         obstime=obstime,
     )
     assert_allclose(bkg.data, bkg_fpi.data, rtol=1e-5)
@@ -213,6 +215,7 @@ def make_map_background_irf_with_symmetry(fpi, symmetry="constant"):
         ontime="42 s",
         bkg=bkg_3d_custom(symmetry),
         geom=WcsGeom.create(npix=(3, 3), binsz=4, axes=[axis], skydir=fpi.fixed_icrs),
+        fov_rotation_error_limit=0.5 * u.deg,
         obstime=obstime,
     )
 
@@ -262,6 +265,7 @@ def test_make_map_background_irf(bkg_3d, pars, fixed_pointing_info):
             ebounds=pars["ebounds"],
             skydir=fixed_pointing_info.fixed_icrs,
         ),
+        fov_rotation_error_limit=0.5 * u.deg,
         oversampling=10,
         obstime=Time("2020-01-01T20:00"),
     )
@@ -322,9 +326,29 @@ def test_make_map_background_irf_skycoord(fixed_pointing_info_aligned):
             ontime="42 s",
             bkg=bkg_3d_custom("asymmetric", "ALTAZ"),
             geom=WcsGeom.create(npix=(3, 3), binsz=4, axes=[axis], skydir=position),
+            fov_rotation_error_limit=0.5 * u.deg,
         )
 
 
+@requires_data()
+def test_make_map_background_irf_AltAz_align(fixed_pointing_info):
+    axis = MapAxis.from_edges([0.1, 1, 10], name="energy", unit="TeV", interp="log")
+    obstime = Time("2020-01-01T20:00:00")
+    make_map_background_irf(
+        pointing=fixed_pointing_info,
+        ontime="42 s",
+        bkg=bkg_3d_custom("asymmetric", "ALTAZ"),
+        geom=WcsGeom.create(
+            npix=(3, 3),
+            binsz=4,
+            axes=[axis],
+            skydir=fixed_pointing_info.get_icrs(obstime),
+        ),
+        fov_rotation_error_limit=0.5 * u.deg,
+        obstime=obstime,
+    )
+
+
 def test_make_edisp_kernel_map():
     migra = MapAxis.from_edges(np.linspace(0.5, 1.5, 50), unit="", name="migra")
     etrue = MapAxis.from_energy_bounds(0.5, 2, 6, unit="TeV", name="energy_true")

gammapy/makers/utils.py

@@ -98,6 +98,38 @@ def _get_fov_coords(pointing, irf, geom, use_region_center=True, obstime=None):
     return coords
 
 
+def time_limit_rotation(maximum_rotation, pointing_altaz):
+    """
+    Compute the maximum time such that the FoV associated to a fixed RaDec position rotates by less than
+    'maximum_rotation' in AltAz frame. It assumes that the rotation rate at the provided pointing is a good estimate
+    of the rotation rate over the full output duration (first order approximation).
+
+
+    Parameters
+    ----------
+    maximum_rotation : `~astropy.units.Quantity`
+        Maximum rotation angle.
+    pointing_altaz : `~astropy.coordinates.SkyCoord`
+        Pointing direction.
+
+    Returns
+    -------
+    duration : `~astropy.units.Quantity`
+        Time associated with the requested rotation.
+    """
+    earth_rotation_deg_per_second = 0.004167
+    denom = (
+        earth_rotation_deg_per_second
+        * np.cos(pointing_altaz.location.lat.rad)
+        * np.cos(pointing_altaz.az.rad)
+    )
+    if denom == 0:
+        return 3600 * u.s
+    return (
+        maximum_rotation.to_value(u.deg) * np.cos(pointing_altaz.alt.rad) / denom
+    ) * u.s
+
+
 def make_map_exposure_true_energy(
     pointing, livetime, aeff, geom, use_region_center=True
 ):
@@ -185,11 +217,26 @@ def _map_spectrum_weight(map, spectrum=None):
     return map * weights.reshape(shape.astype(int))
 
 
+def evaluate_bkg(pointing, ontime, bkg, geom, use_region_center, obstime, d_omega):
+    coords = _get_fov_coords(
+        pointing=pointing,
+        irf=bkg,
+        geom=geom,
+        use_region_center=use_region_center,
+        obstime=obstime,
+    )
+    coords["energy"] = broadcast_axis_values_to_geom(geom, "energy", False)
+
+    bkg_de = bkg.integrate_log_log(**coords, axis_name="energy")
+    return (bkg_de * d_omega * ontime).to_value("")
+
+
 def make_map_background_irf(
     pointing,
     ontime,
     bkg,
     geom,
+    fov_rotation_error_limit,
     oversampling=None,
     use_region_center=True,
     obstime=None,
@@ -213,6 +260,8 @@ def make_map_background_irf(
         Background rate model.
     geom : `~gammapy.maps.WcsGeom`
         Reference geometry.
+    fov_rotation_error_limit : `~astropy.units.Quantity`
+        Maximum rotation error to create sub-time bins if the irf is 3D and in AltAz.
     oversampling : int
         Oversampling factor in energy, used for the background model evaluation.
     use_region_center : bool, optional
@@ -228,11 +277,8 @@ def make_map_background_irf(
         Background predicted counts sky cube in reconstructed energy.
     """
     # TODO:
-    #  This implementation can be improved in two ways:
-    #  1. Create equal time intervals between TSTART and TSTOP and sum up the
-    #  background IRF for each interval. This is instead of multiplying by
-    #  the total ontime. This then handles the rotation of the FoV.
-    #  2. Use the pointing table (does not currently exist in CTA files) to
+    #  This implementation can be improved in one ways:
+    #  Use the pointing table (does not currently exist in CTA files) to
     #  obtain the RA DEC and time for each interval. This then considers that
     #  the pointing might change slightly over the observation duration
 
@@ -248,18 +294,31 @@ def make_map_background_irf(
     else:
         image_geom = geom.to_image()
         d_omega = image_geom.solid_angle()
-
-    coords = _get_fov_coords(
-        pointing=pointing,
-        irf=bkg,
-        geom=geom,
-        use_region_center=use_region_center,
-        obstime=obstime,
-    )
-    coords["energy"] = broadcast_axis_values_to_geom(geom, "energy", False)
-
-    bkg_de = bkg.integrate_log_log(**coords, axis_name="energy")
-    data = (bkg_de * d_omega * ontime).to_value("")
+    if bkg.has_offset_axis or bkg.fov_alignment == FoVAlignment.RADEC:
+        data = evaluate_bkg(
+            pointing, ontime, bkg, geom, use_region_center, obstime, d_omega
+        )
+    if bkg.fov_alignment == FoVAlignment.ALTAZ:
+        endtime = obstime + ontime
+        data = np.zeros(geom.data_shape)
+        while obstime < endtime:
+            # Evaluate the step time needed to have FoV rotation of less than fov_rotation_error_limit
+            # Minimum step of 1 second to avoid infinite computation very close to zenith
+            steptime = max(
+                time_limit_rotation(
+                    fov_rotation_error_limit, pointing.get_altaz(obstime)
+                ),
+                1 * u.s,
+            )
+            steptime = (
+                steptime
+                if steptime < endtime - obstime - 1 * u.s
+                else endtime - obstime
+            )
+            data += evaluate_bkg(
+                pointing, steptime, bkg, geom, use_region_center, obstime, d_omega
+            )
+            obstime = obstime + steptime
 
     if not use_region_center:
         region_coord, weights = geom.get_wcs_coord_and_weights()