This is a fixed-text formatted version of a Jupyter notebook.

# Template background model production with Gammapy¶

## Introduction¶

In this tutorial, we will create a template background model in the bkg_2d format, i.e. with offset and energy axes (see spec).

We will be working with only 4 H.E.S.S. runs on the Crab nebula here, just as an example.

To build a coherent background model you normally use 100s of runs of AGN observations or intentional “off” runs that point at parts of the sky containing no known gamma-ray sources.

We will mainly be using the following classes:

## Setup¶

As always, we start the notebook with some setup and imports.

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt

In [2]:

import shutil
import numpy as np
import astropy.units as u
from astropy.table import Table
from astropy.coordinates import SkyCoord, Angle

In [3]:

from gammapy.extern.pathlib import Path
from gammapy.utils.energy import EnergyBounds
from gammapy.utils.nddata import sqrt_space
from gammapy.data import DataStore, ObservationGroupAxis, ObservationGroups
from gammapy.background import EnergyOffsetBackgroundModel
from gammapy.background import OffDataBackgroundMaker
from gammapy.catalog import SourceCatalogGammaCat


## Compute background model¶

Computing a set of template background model has two major steps: 1. Define group of runs for each background model 2. Run the OffDataBackgroundMaker

We also need a scratch directory, and a table of known gamma-ray sources to exclude.

### Make a scratch directory¶

Background model production is a little pipeline that needs a “scratch” directory to put some files while running. Let’s make ourselves a fresh empty scratch sub-directory called background in the current working directory.

In [4]:

def make_fresh_dir(path):
"""Make a fresh directory. Delete first if exists"""
path = Path(path)
if path.is_dir():
shutil.rmtree(str(path))
path.mkdir()
return path

In [5]:

scratch_dir = make_fresh_dir('background')
scratch_dir

Out[5]:

PosixPath('background')


### Make an observation table defining the run grouping¶

Prepare a scheme to group observations with similar observing conditions and create a new ObservationTable with the grouping ID for each run

In [6]:

# Create a background model from the 4 Crab runs for the counts ouside the exclusion region so here outside the Crab
data_store = DataStore.from_dir("$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2") # Define the grouping you want to use to group the obervations to make the acceptance curves # Here we use 2 Zenith angle bins only, you can also add efficiency bins for example etc... axes = [ObservationGroupAxis('ZEN_PNT', [0, 49, 90], fmt='edges')] # Create the ObservationGroups object obs_groups = ObservationGroups(axes) # write it to file filename = str(scratch_dir / 'group-def.fits') obs_groups.obs_groups_table.write(filename, overwrite=True) # Create a new ObservationTable with the column group_id # You give the runs list you want to use to produce the background model that are in your obs table. # Here very simple only the 4 Crab runs... list_ids = [23523, 23526, 23559, 23592] obs_table_with_group_id = obs_groups.apply(data_store.obs_table.select_obs_id(list_ids))  ### Make table of known gamma-ray sources to exclude¶ We need a mask to remove known sources from the observation. We use gamma-cat and exclude a circular region of at least 0.3° radius. Here since we use Crab runs, we will remove the Crab events from the FOV to select only the OFF events to build the acceptance curves. Of cource normally you use thousand of AGN runs to build coherent acceptance curves. In [7]:  cat = SourceCatalogGammaCat() exclusion_table = cat.table.copy() exclusion_table.rename_column('ra', 'RA') exclusion_table.rename_column('dec', 'DEC') radius = exclusion_table['morph_sigma'].data radius[np.isnan(radius)] = 0.3 exclusion_table['Radius'] = radius * u.deg exclusion_table = Table(exclusion_table)  ### Run the OffDataBackgroundMaker¶ Make the acceptance curves in the different group of observation conditions you defined above using the obs_table containaing the group id for each observation used to compute the bkg model In [8]:  bgmaker = OffDataBackgroundMaker( data_store=data_store, outdir=str(scratch_dir), run_list=None, obs_table=obs_table_with_group_id, ntot_group=obs_groups.n_groups, excluded_sources=exclusion_table, ) # Define the energy and offset binning to use ebounds = EnergyBounds.equal_log_spacing(0.1, 100, 15, 'TeV') offset = sqrt_space(start=0, stop=2.5, num=100) * u.deg # Make the model (i.e. stack counts and livetime) bgmaker.make_model("2D", ebounds=ebounds, offset=offset) # Smooth the model bgmaker.smooth_models("2D") # Write the model to disk bgmaker.save_models("2D") bgmaker.save_models(modeltype="2D", smooth=True)  WARNING: AstropyDeprecationWarning: The truth value of a Quantity is ambiguous. In the future this will raise a ValueError. [astropy.units.quantity]  Congratulations, you have produced a background model. The following files were generated in our scratch directory: In [9]:  [path.name for path in scratch_dir.glob('*')]  Out[9]:  ['smooth_background_2D_group_001_table.fits.gz', 'background_2D_group_001_table.fits.gz', 'group-def.fits', 'background_2D_group_000_table.fits.gz', 'smooth_background_2D_group_000_table.fits.gz']  ## Inspect the background model¶ Our template background model has two axes: offset and energy. Let’s make a few plots to see what it looks like: 1. Acceptance curve (background rate as a function of field of view offset for a given energy) 1. Rate spectrum (background rate as a function of energy for a given offset) 1. Rate image (background rate as a function of energy and offset) ### Acceptance curve¶ In [10]:  # Read one of the background models from file filename = scratch_dir / 'smooth_background_2D_group_000_table.fits.gz' model = EnergyOffsetBackgroundModel.read(str(filename))  In [11]:  offset = model.bg_rate.offset_bin_center energies = model.bg_rate.energy iE = 6 x = offset y = model.bg_rate.data[iE,:] plt.plot(x, y, label="bkg model smooth") title = "energy band: "+str("%.2f"%energies[iE].value)+"-"+str("%.2f"%energies[iE+1].value)+" TeV" plt.title(title) plt.xlabel("Offset (degree)") plt.ylabel("Bkg rate (MeV-1 s-1 sr-1)") plt.legend()  Out[11]:  <matplotlib.legend.Legend at 0x1130ab710>  ### Background rate spectrum¶ In [12]:  x = energies.log_centers y = model.bg_rate.data[:,10] plt.loglog(x, y, label="bkg model smooth") plt.title("offset: "+str("%.2f"%offset[10].value)+" deg") plt.xlabel("Energy (TeV)") plt.ylabel("Bkg rate (MeV-1 s-1 sr-1)")  Out[12]:  Text(0,0.5,'Bkg rate (MeV-1 s-1 sr-1)')  ### Background rate image with energy and offset axes¶ It doesn’t look good in this case. To do this well, you need to use more off or AGN runs to build the background model! In [13]:  model.bg_rate.plot()  Out[13]:  (<matplotlib.axes._subplots.AxesSubplot at 0x1139aebe0>, <matplotlib.colorbar.Colorbar at 0x113a5df60>)  ## Make new index table¶ Here we first copy the dataset of the 4 crab runs from gammapy-extra in a new directory containing the data you will use for the analysis. We use the same dataset to produce the bkg or for the analysis. Of course normally you produce the bkg model using thousands of AGN runs not the 4 Crab test runs. In [14]:  # Make a new hdu table in your dataset directory that contains the link to the acceptance curve to use to build the bkg model in your cube analysis data_dir = make_fresh_dir('data')  In [15]:  ds = DataStore.from_dir("$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2")
ds.copy_obs(ds.obs_table, data_dir)

WARNING: Attribute BASE_DIR of type <class 'gammapy.extern.pathlib.PosixPath'> cannot be added to FITS Header - skipping [astropy.io.fits.convenience]


The hdu_table in this directory contains no link to a bkg model for each observation.

In [16]:

data_store = DataStore.from_dir(data_dir)
data_store.hdu_table

Out[16]:

HDUIndexTable length=28
OBS_IDHDU_TYPEHDU_CLASSFILE_DIRFILE_NAMEHDU_NAMESIZEMTIMEMD5
int64bytes6bytes10bytes26bytes30bytes12int64float64bytes32
23523gtigtirun023400-023599/run023523hess_events_023523.fits.gzGTI6209751455089616.34075029e402094c3a3e05ae4199b7cc9a01215
23523eventseventsrun023400-023599/run023523hess_events_023523.fits.gzEVENTS6209751455089616.34075029e402094c3a3e05ae4199b7cc9a01215
23523aeffaeff_2drun023400-023599/run023523hess_aeff_2d_023523.fits.gzAEFF_2D37271455089616.34075026430c082176f092e0aed0f2bf9840915
23523edispedisp_2drun023400-023599/run023523hess_edisp_2d_023523.fits.gzEDISP_2D289631455089616.3407502f580ea6cb104e4d6735b8d2940ac6774
23523psfpsf_3gaussrun023400-023599/run023523hess_psf_3gauss_023523.fits.gzPSF_2D_GAUSS30271455089616.344750287f2d5c5ca56575a4a083b33e9700312
23523psfpsf_kingrun023400-023599/run023523hess_psf_king_023523.fits.gzPSF_2D_KING18231455089616.34475027760e349a40883345406c7e3ea1cbd54
23526psfpsf_kingrun023400-023599/run023526hess_psf_king_023526.fits.gzPSF_2D_KING18341455089616.5647545873e4ec0771bdfcc6d82608199a4067
23526psfpsf_3gaussrun023400-023599/run023526hess_psf_3gauss_023526.fits.gzPSF_2D_GAUSS30021455089616.564754f27364b40bbf8e35c747828e60224b28
23526edispedisp_2drun023400-023599/run023526hess_edisp_2d_023526.fits.gzEDISP_2D288821455089616.560754c7e99f4d282a55c7fdd4024167bdbef5
...........................
23559aeffaeff_2drun023400-023599/run023559hess_aeff_2d_023559.fits.gzAEFF_2D37151455089616.77675850eb94ec65ac146ff60845860d84a1a2
23592psfpsf_kingrun023400-023599/run023592hess_psf_king_023592.fits.gzPSF_2D_KING18181455089616.9967618514f94397507d1160341387ed56e80ba
23592gtigtirun023400-023599/run023592hess_events_023592.fits.gzGTI5981071455089616.98876173c94433be8e9f29aa198065403570dbb
23592eventseventsrun023400-023599/run023592hess_events_023592.fits.gzEVENTS5981071455089616.98876173c94433be8e9f29aa198065403570dbb
23592aeffaeff_2drun023400-023599/run023592hess_aeff_2d_023592.fits.gzAEFF_2D37211455089616.99276165e2f5dc70fecb5060b6e1662d59d9507
23592edispedisp_2drun023400-023599/run023592hess_edisp_2d_023592.fits.gzEDISP_2D289311455089616.992761637db40409b6c37a9060aef8817a14ffe
23592psfpsf_tablerun023400-023599/run023592hess_psf_table_023592.fits.gzPSF_2D_TABLE2217091455089617.000761738be05a13aa0bd492c9979c5b7d3d47f

In order to produce a background image or background cube we have to create a hdu table that contains for each observation a link to the bkg model to use depending of the observation conditions of the run.

In [17]:

#Copy the background directory in the one where is located the hdu table, here data
shutil.move(str(scratch_dir), str(data_dir))

# Create the new hdu table with a link to the background model
group_filename = data_dir / 'background/group-def.fits'

#relat_path= (scratch_dir.absolute()).relative_to(data_dir.absolute())
hdu_index_table = bgmaker.make_total_index_table(
data_store=data_store,
modeltype='2D',
out_dir_background_model=scratch_dir,
filename_obs_group_table=str(group_filename),
smooth=False,
)

# Write the new hdu table
filename = data_dir / 'hdu-index.fits.gz'
hdu_index_table.write(str(filename), overwrite=True)

WARNING: Attribute BASE_DIR of type <class 'gammapy.extern.pathlib.PosixPath'> cannot be added to FITS Header - skipping [astropy.io.fits.convenience]

In [18]:

hdu_index_table

Out[18]:

OBS_IDHDU_TYPEHDU_CLASSFILE_DIRFILE_NAMEHDU_NAMESIZEMTIMEMD5
int64str6str10str26str37str12int64float64bytes32
23523gtigtirun023400-023599/run023523hess_events_023523.fits.gzGTI6209751455089616.34075029e402094c3a3e05ae4199b7cc9a01215
23523eventseventsrun023400-023599/run023523hess_events_023523.fits.gzEVENTS6209751455089616.34075029e402094c3a3e05ae4199b7cc9a01215
23523aeffaeff_2drun023400-023599/run023523hess_aeff_2d_023523.fits.gzAEFF_2D37271455089616.34075026430c082176f092e0aed0f2bf9840915
23523edispedisp_2drun023400-023599/run023523hess_edisp_2d_023523.fits.gzEDISP_2D289631455089616.3407502f580ea6cb104e4d6735b8d2940ac6774
23523psfpsf_3gaussrun023400-023599/run023523hess_psf_3gauss_023523.fits.gzPSF_2D_GAUSS30271455089616.344750287f2d5c5ca56575a4a083b33e9700312
23523psfpsf_kingrun023400-023599/run023523hess_psf_king_023523.fits.gzPSF_2D_KING18231455089616.34475027760e349a40883345406c7e3ea1cbd54
23526psfpsf_kingrun023400-023599/run023526hess_psf_king_023526.fits.gzPSF_2D_KING18341455089616.5647545873e4ec0771bdfcc6d82608199a4067
23526psfpsf_3gaussrun023400-023599/run023526hess_psf_3gauss_023526.fits.gzPSF_2D_GAUSS30021455089616.564754f27364b40bbf8e35c747828e60224b28
23526edispedisp_2drun023400-023599/run023526hess_edisp_2d_023526.fits.gzEDISP_2D288821455089616.560754c7e99f4d282a55c7fdd4024167bdbef5
...........................
23592gtigtirun023400-023599/run023592hess_events_023592.fits.gzGTI5981071455089616.98876173c94433be8e9f29aa198065403570dbb
23592eventseventsrun023400-023599/run023592hess_events_023592.fits.gzEVENTS5981071455089616.98876173c94433be8e9f29aa198065403570dbb
23592aeffaeff_2drun023400-023599/run023592hess_aeff_2d_023592.fits.gzAEFF_2D37211455089616.99276165e2f5dc70fecb5060b6e1662d59d9507
23592edispedisp_2drun023400-023599/run023592hess_edisp_2d_023592.fits.gzEDISP_2D289311455089616.992761637db40409b6c37a9060aef8817a14ffe
23592psfpsf_tablerun023400-023599/run023592hess_psf_table_023592.fits.gzPSF_2D_TABLE2217091455089617.000761738be05a13aa0bd492c9979c5b7d3d47f
23526bkgbkg_2dbackgroundbackground_2D_group_000_table.fits.gzbkg_2d------
23559bkgbkg_2dbackgroundbackground_2D_group_000_table.fits.gzbkg_2d------
23592bkgbkg_2dbackgroundbackground_2D_group_000_table.fits.gzbkg_2d------
23523bkgbkg_2dbackgroundbackground_2D_group_001_table.fits.gzbkg_2d------
In [19]:

print(hdu_index_table[-4]["FILE_DIR"], " ", hdu_index_table[-4]["FILE_NAME"])

background   background_2D_group_000_table.fits.gz

In [20]:

print(hdu_index_table[-1]["FILE_DIR"], " ", hdu_index_table[-1]["FILE_NAME"])

background   background_2D_group_001_table.fits.gz


## Exercises¶

• Use real AGN run
• Change the binning for the grouping: thinner zenithal bin, add efficiency binning ….
• Change the energy binning (ebounds) and the offset (offset) used to compute the acceptance curve

## What next?¶

In this tutorial we have created a template background model in the bkg_2d format, i.e. with offset and energy axes (see spec).

In future tutorials, we will use this background model as one of the model components for source analysis.