# 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].

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
from gammapy.utils.regions import make_region

data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1") observations = data_store.get_observations([23592, 23559]) on_region = make_region("icrs; circle(83.63,22.01,0.11)") wcsgeom = WcsGeom.create(skydir=on_region.center, width=5, binsz=0.02) exclusion_mask = wcsgeom.region_mask([on_region], inside=False) 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=on_region, axes=[energy_axis]) dataset_empty = SpectrumDataset.create(geom=geom, energy_axis_true=energy_axis_true) 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. import matplotlib.pyplot as plt from astropy.coordinates import Angle, SkyCoord from regions import CircleSkyRegion from gammapy.makers import ReflectedRegionsFinder from gammapy.maps import WcsNDMap, RegionGeom # 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( region=on_region, center=center, exclusion_mask=exclusion_mask, min_distance_input="0.2 rad", ) regions = finder.run() 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( region=on_region, center=center, exclusion_mask=exclusion_mask, min_distance="0.1 rad", ) regions = finder.run() 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( region=on_region, center=center, exclusion_mask=exclusion_mask, max_region_number=5, min_distance="0.1 rad", ) regions = finder.run() 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) ## 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/")
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
)

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()


(png, hires.png, pdf)

The following notebook shows an example using ReflectedRegionsBackgroundMaker to perform a spectral extraction and fitting: