Reflected regions background#

Overview#

This method is used in classical Cherenkov astronomy to estimate background for spectral analysis. As illustrated in Fig. 1, a region on the sky, the ON region, is chosen to select events around the studied source position. In the absence of a solid template of the residual hadronic background, a classical method to estimate it is the so-called Reflected Region Background. The underlying assumption is that the background is approximately purely radial in the field-of-view. A set of OFF counts is found in the observation, by rotating the ON region selected around the pointing position. To avoid that the reflected regions contain actual gamma-ray signal from other objects, one has to remove the gamma-ray bright parts of the field-of-view with a exclusion mask. More details on the reflected regions method can be found in [Berge2007].

../../_images/hgps_spectrum_background_estimation.png — Fig.1, Illustration of the reflected regions background estimation method, taken from [Abdalla2018].#

The extraction of the OFF events from the EventList of a set of observations is performed by the ReflectedRegionsBackgroundMaker. The latter uses the ReflectedRegionsFinder to create reflected regions for a given circular on region and exclusion mask.

from gammapy.makers import SpectrumDatasetMaker, ReflectedRegionsBackgroundMaker, SafeMaskMaker
from gammapy.datasets import SpectrumDataset, Datasets
from gammapy.data import DataStore
from gammapy.maps import MapAxis, Map, WcsGeom, RegionGeom

data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1")
observations = data_store.get_observations([23592, 23559])

energy_axis = MapAxis.from_energy_bounds("0.5 TeV", "10 TeV", nbin=15)
energy_axis_true = MapAxis.from_energy_bounds("0.3 TeV", "20 TeV", nbin=40, name="energy_true")

geom = RegionGeom.create(region="icrs;circle(83.63,22.01,0.11)", axes=[energy_axis])
dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true)

wcsgeom = WcsGeom.create(skydir=geom.center_skydir, width=5, binsz=0.02)
exclusion_mask = wcsgeom.region_mask(geom.region, inside=False)

maker = SpectrumDatasetMaker()
safe_mask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)

bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask)

datasets = Datasets()
for obs in observations:
    dataset = maker.run(dataset_empty.copy(), obs)
    dataset = safe_mask_maker.run(dataset, obs)
    dataset_on_off = bkg_maker.run(dataset, obs)
    datasets.append(dataset)

Using regions#

The on region is a SkyRegion. It is typically a circle (CircleSkyRegion) for pointlike source analysis, but can be a more complex region such as a CircleAnnulusSkyRegion a EllipseSkyRegion, a RectangleSkyRegion etc.

The following example shows how to create such regions:

"""Example how to compute and plot reflected regions."""
import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import CircleSkyRegion, EllipseAnnulusSkyRegion, RectangleSkyRegion
import matplotlib.pyplot as plt
from gammapy.maps import WcsNDMap

position = SkyCoord(83.63, 22.01, unit="deg", frame="icrs")

on_circle = CircleSkyRegion(position, 0.3 * u.deg)

on_ellipse_annulus = EllipseAnnulusSkyRegion(
    center=position,
    inner_width=1.5 * u.deg,
    outer_width=2.5 * u.deg,
    inner_height=3 * u.deg,
    outer_height=4 * u.deg,
    angle=130 * u.deg,
)

another_position = SkyCoord(80.3, 22.0, unit="deg")
on_rectangle = RectangleSkyRegion(
    center=another_position, width=2.0 * u.deg, height=4.0 * u.deg, angle=50 * u.deg
)

# Now we plot those regions. We first create an empty map
empty_map = WcsNDMap.create(
    skydir=position, width=10 * u.deg, binsz=0.1 * u.deg, proj="TAN"
)
empty_map.data += 1.0
empty_map.plot(cmap="gray", vmin=0, vmax=1)

# To plot the regions, we convert them to PixelRegion with the map wcs
on_circle.to_pixel(empty_map.geom.wcs).plot()
on_rectangle.to_pixel(empty_map.geom.wcs).plot()
on_ellipse_annulus.to_pixel(empty_map.geom.wcs).plot()

plt.show()

(png, hires.png, pdf)

The reflected region finder#

The following example illustrates how to create reflected regions for a given circular on region and exclusion mask using the ReflectedRegionsFinder. In particular, it shows how to change the minimal distance between the ON region and the reflected regions. This is useful to limit contamination by events leaking out the ON region. It also shows how to change the minimum distance between adjacent regions as well as the maximum number of reflected regions.

from astropy.coordinates import Angle, SkyCoord
from regions import CircleSkyRegion
import matplotlib.pyplot as plt
from gammapy.makers import ReflectedRegionsFinder
from gammapy.maps import RegionGeom, WcsNDMap

# Exclude a rectangular region
exclusion_mask = WcsNDMap.create(npix=(801, 701), binsz=0.01, skydir=(83.6, 23.0))

coords = exclusion_mask.geom.get_coord().skycoord
data = (Angle("23 deg") < coords.dec) & (coords.dec < Angle("24 deg"))
exclusion_mask.data = ~data

pos = SkyCoord(83.633, 22.014, unit="deg")
radius = Angle(0.3, "deg")
on_region = CircleSkyRegion(pos, radius)
center = SkyCoord(83.633, 24, unit="deg")

# One can impose a minimal distance between ON region and first reflected regions
finder = ReflectedRegionsFinder(
    min_distance_input="0.2 rad",
)
regions, _ = finder.run(
    region=on_region,
    center=center,
    exclusion_mask=exclusion_mask,
)

fig, axes = plt.subplots(
    ncols=3,
    subplot_kw={"projection": exclusion_mask.geom.wcs},
    figsize=(12, 3),
)


def plot_regions(ax, regions, on_region, exclusion_mask):
    """Little helper function to plot off regions"""
    exclusion_mask.plot_mask(ax=ax, colors="gray")
    on_region.to_pixel(ax.wcs).plot(ax=ax, edgecolor="tab:orange")
    geom = RegionGeom.from_regions(regions)
    geom.plot_region(ax=ax, color="tab:blue")


ax = axes[0]
ax.set_title("Min. distance first region")
plot_regions(ax=ax, regions=regions, on_region=on_region, exclusion_mask=exclusion_mask)


# One can impose a minimal distance between two adjacent regions
finder = ReflectedRegionsFinder(
    min_distance="0.1 rad",
)
regions, _ = finder.run(
    region=on_region,
    center=center,
    exclusion_mask=exclusion_mask,
)

ax = axes[1]
ax.set_title("Min. distance all regions")
plot_regions(ax=ax, regions=regions, on_region=on_region, exclusion_mask=exclusion_mask)


# One can impose a maximal number of regions to be extracted
finder = ReflectedRegionsFinder(
    max_region_number=5,
    min_distance="0.1 rad",
)
regions, _ = finder.run(
    region=on_region,
    center=center,
    exclusion_mask=exclusion_mask,
)

ax = axes[2]
ax.set_title("Max. number of regions")
plot_regions(ax=ax, regions=regions, on_region=on_region, exclusion_mask=exclusion_mask)
plt.show()

(png, hires.png, pdf)

../../_images/make_reflected_regions.png

Using the reflected background estimator#

In practice, the user does not usually need to directly interact with the ReflectedRegionsFinder. This actually is done via the ReflectedRegionsBackgroundMaker, which extracts the ON and OFF events for an Observations object. The last example shows how to run it on a few observations with a rectangular region.

"""Example how to compute and plot reflected regions."""
import astropy.units as u
from astropy.coordinates import SkyCoord
from regions import RectangleSkyRegion
import matplotlib.pyplot as plt
from gammapy.data import DataStore
from gammapy.datasets import SpectrumDataset
from gammapy.makers import ReflectedRegionsBackgroundMaker, SpectrumDatasetMaker
from gammapy.maps import Map, MapAxis, RegionGeom
from gammapy.visualization import plot_spectrum_datasets_off_regions

data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1/")
mask = data_store.obs_table["TARGET_NAME"] == "Crab"
obs_ids = data_store.obs_table["OBS_ID"][mask].data
observations = data_store.get_observations(obs_ids)

crab_position = SkyCoord(83.63, 22.01, unit="deg", frame="icrs")

# The ON region center is defined in the icrs frame. The angle is defined w.r.t. to its axis.
rectangle = RectangleSkyRegion(
    center=crab_position, width=0.5 * u.deg, height=0.4 * u.deg, angle=0 * u.deg
)

bkg_maker = ReflectedRegionsBackgroundMaker(min_distance=0.1 * u.rad)
dataset_maker = SpectrumDatasetMaker(selection=["counts"])

energy_axis = MapAxis.from_energy_bounds(0.1, 100, 30, unit="TeV")
geom = RegionGeom.create(region=rectangle, axes=[energy_axis])
dataset_empty = SpectrumDataset.create(geom=geom)

datasets = []

for obs in observations:

    dataset = dataset_maker.run(dataset_empty.copy(name=f"obs-{obs.obs_id}"), obs)
    dataset_on_off = bkg_maker.run(observation=obs, dataset=dataset)
    datasets.append(dataset_on_off)

m = Map.create(skydir=crab_position, width=(8, 8), proj="TAN")

ax = m.plot(vmin=-1, vmax=0)

rectangle.to_pixel(ax.wcs).plot(ax=ax, color="black")

plot_spectrum_datasets_off_regions(datasets=datasets, ax=ax)
plt.show()