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

CTA simulation tools

Introduction

In this tutorial we will simulate the expected counts of a Fermi/LAT source in the CTA energy range.

We will go through the following topics:

Setup

In order to deal with plots we will begin with matplotlib import:

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt

PKS 2155-304 selection from the 3FHL Fermi/LAT catalogue

We will start by selecting the source PKS 2155-304 in the 3FHL Fermi/LAT catalogue for further use.

In [2]:
from gammapy.catalog import SourceCatalog3FHL

# load catalogs
fermi_3fhl = SourceCatalog3FHL()
name = 'PKS 2155-304'
source = fermi_3fhl[name]

We can then access the caracteristics of the source via the data attribut and select its spectral model for further use.

In [3]:
redshift = source.data['Redshift']
src_spectral_model = source.spectral_model

Here is an example on how to plot the source spectra

In [4]:
# plot the Fermi/LAT model
import astropy.units as u
src_spectral_model.plot(energy_range=[10 * u.GeV, 2 *u.TeV])
Out[4]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f1641669828>
../_images/notebooks_cta_simulation_8_1.png

Select a model for EBL absorption

We will need to modelise EBL (extragalactic background light) attenuation to have get a ‘realistic’ simulation. Different models are available in GammaPy. Here is an example on how to deal with the absorption coefficients.

In [5]:
from gammapy.spectrum.models import Absorption

# Load models for PKS 2155-304 redshift
dominguez = Absorption.read_builtin('dominguez').table_model(redshift)
franceschini = Absorption.read_builtin('franceschini').table_model(redshift)
finke = Absorption.read_builtin('finke').table_model(redshift)

From here you can have access to the absorption coefficient for a given energy.

In [6]:
energy = 1 * u.TeV
abs_value = dominguez.evaluate(energy=energy, scale=1)
print('absorption({} {}) = {}'.format(energy.value, energy.unit, abs_value))
absorption(1.0 TeV) = 0.2877190690093156

Below is an example to plot EBL absorption for different models

In [7]:
# start customised plot
energy_range = [0.08, 3] * u.TeV
ax = plt.gca()
opts = dict(energy_range=energy_range, energy_unit='TeV', ax=ax)
franceschini.plot(label='Franceschini 2008', **opts)
finke.plot(label='Finke 2010', **opts)
dominguez.plot(label='Dominguez 2011', **opts)

# tune plot
ax.set_ylabel(r'Absorption coefficient [$\exp{(-        au(E))}$]')
ax.set_xlim(energy_range.value)  # we get ride of units
ax.set_ylim([1.e-1, 2.])
ax.set_yscale('log')
ax.set_title('EBL models (z=' + str(redshift) + ')')
plt.grid(which='both')
plt.legend(loc='best') # legend

# show plot
plt.show()

../_images/notebooks_cta_simulation_14_0.png

CTA instrument response functions

Here we are going to deal with CTA point-like instrument response functions (public version, production 2). Data format for point-like IRF is still missing. For now, a lot of efforts is made to define full-containment IRFs (https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/index.html). In the meantime a temporary format is used in gammapy. It will evolved.

To simulate one observation we need the following IRFs: - effective area as a function of true energy (energy-dependent theta square cute) - background rate as a function of reconstructed energy (energy-dependent theta square cute) - migration matrix, e_reco/e_true as a function of true energy

To handle CTA’s responses we will use the CTAPerf class

In [8]:
from gammapy.scripts import CTAPerf
# South array optimisation for faint source
filename = '$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/South_50h.fits.gz'
cta_perf = CTAPerf.read(filename)
cta_perf.peek()
/home/kingj/Software/miniconda3/envs/headversions/lib/python3.5/site-packages/astropy-3.0.dev20201-py3.5-linux-x86_64.egg/astropy/units/quantity.py:1053: RuntimeWarning: invalid value encountered in true_divide
  return super().__truediv__(other)
../_images/notebooks_cta_simulation_16_1.png

Different optimisations are available for different type of source (bright, 0.5h; medium, 5h; faint, 50h). Here is an example to have a quick look to the different optimisation

In [9]:
prod_dir = '$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/'
opti = ['0.5h', '5h', '50h']
site = ['South', 'North']
cta_perf_list = []  # will be filled with different performance
labels = []  # will be filled with different performance labels for the legend
for isite in site:
    for iopti in opti:
        filename = prod_dir + '/' + isite + '_' + iopti + '.fits.gz'
        cta_perf = CTAPerf.read(filename)
        cta_perf_list.append(cta_perf)
        labels.append(isite + ' (' + iopti + ')')

CTAPerf.superpose_perf(cta_perf_list, labels)
/home/kingj/Software/miniconda3/envs/headversions/lib/python3.5/site-packages/astropy-3.0.dev20201-py3.5-linux-x86_64.egg/astropy/units/quantity.py:1053: RuntimeWarning: invalid value encountered in true_divide
  return super().__truediv__(other)
../_images/notebooks_cta_simulation_18_1.png

CTA simulation of an observation

Now we are going to simulate the expected counts in the CTA energy range. To do so we will need to specify a target (caracteristics of the source) and the parameters of the observation (such as time, ON/OFF normalisation, etc.)

Target definition

In [10]:
# define target spectral model absorbed by EBL
from gammapy.spectrum.models import Absorption, AbsorbedSpectralModel

absorption = Absorption.read_builtin('dominguez')
spectral_model = AbsorbedSpectralModel(
    spectral_model=src_spectral_model,
    absorption=absorption,
    parameter=redshift,
)
# define target
from gammapy.scripts.cta_utils import Target
target = Target(
    name=source.data['Source_Name'],  # from the 3FGL catalogue source class
    model=source.spectral_model,  # defined above
)
print(target)
*** Target parameters ***
Name=3FHL J2158.8-3013
amplitude=7.707001703494143e-11 1 / (cm2 GeV s)
reference=18.31732177734375 GeV
alpha=1.8807274103164673
beta=0.14969758689403534

Observation definition

In [11]:
from gammapy.scripts.cta_utils import ObservationParameters
alpha = 0.2 * u.Unit('')  # normalisation between ON and OFF regions
livetime = 5. * u.h
# energy range used for statistics (excess, significance, etc.)
emin, emax = 0.05 * u.TeV, 5 * u.TeV
params = ObservationParameters(
    alpha=alpha, livetime=livetime,
    emin=emin, emax=emax,
)
print(params)
*** Observation parameters summary ***
alpha=0.2 []
livetime=5.0 [h]
emin=0.05 [TeV]
emax=5.0 [TeV]

Performance

In [12]:
from gammapy.scripts import CTAPerf
# PKS 2155-304 is 10 % of Crab at 1 TeV ==> intermediate source
filename = '$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/South_5h.fits.gz'
perf = CTAPerf.read(filename)
/home/kingj/Software/miniconda3/envs/headversions/lib/python3.5/site-packages/astropy-3.0.dev20201-py3.5-linux-x86_64.egg/astropy/units/quantity.py:1053: RuntimeWarning: invalid value encountered in true_divide
  return super().__truediv__(other)

Simulation

Here we are going to simulate what we expect to see with CTA and measure the duration of the simulation

In [13]:
from gammapy.scripts.cta_utils import CTAObservationSimulation

simu = CTAObservationSimulation.simulate_obs(
    perf=perf,
    target=target,
    obs_param=params,
)

# print simulation results
print(simu)
*** Observation summary report ***
Observation Id: 0
Livetime: 5.000 h
On events: 9619
Off events: 10400
Alpha: 0.200
Bkg events in On region: 2080.00
Excess: 7539.00
Excess / Background: 3.62
Gamma rate: 2.51 1 / min
Bkg rate: 0.69 1 / min
Sigma: 102.67
energy range: 0.05 TeV - 5.01 TeV

Now we can take a look at the excess, ON and OFF distributions

In [14]:
CTAObservationSimulation.plot_simu(simu, target)
../_images/notebooks_cta_simulation_30_0.png

We can access simulation parameters via the gammapy.spectrum.SpectrumStats attribute of the gammapy.spectrum.SpectrumObservation class:

In [15]:
stats = simu.total_stats_safe_range
stats_dict = stats.to_dict()
print('excess: {}'.format(stats_dict['excess']))
print('sigma: {:.1f}'.format(stats_dict['sigma']))
excess: 7539.0
sigma: 102.7

Finally, you can get statistics for every reconstructed energy bin with:

In [16]:
table = simu.stats_table()
# Here we only print part of the data from the table
table[['energy_min', 'energy_max', 'excess', 'background', 'sigma']][:10]
/home/kingj/Software/gammapy/gammapy/stats/poisson.py:382: RuntimeWarning: divide by zero encountered in double_scalars
  tt = (alpha + 1) / (n_on + n_off)
/home/kingj/Software/gammapy/gammapy/stats/poisson.py:383: RuntimeWarning: invalid value encountered in multiply
  ll = n_on * np.log(n_on * tt / alpha)
/home/kingj/Software/gammapy/gammapy/stats/poisson.py:384: RuntimeWarning: invalid value encountered in multiply
  mm = n_off * np.log(n_off * tt)
/home/kingj/Software/gammapy/gammapy/stats/poisson.py:384: RuntimeWarning: divide by zero encountered in log
  mm = n_off * np.log(n_off * tt)
Out[16]:
Table length=10
energy_minenergy_maxexcessbackgroundsigma
TeVTeV
float64float64float64float64float64
0.01258925441650.0199526231736124.43701.61.85449984597
0.01995262317360.0316227786243504.03128.07.98444463023
0.03162277862430.05011872574691244.86838.213.2886502977
0.05011872574690.07943282276391569.4818.638.9772009509
0.07943282276390.1258925497531803.2899.842.3713164745
0.1258925497530.1995262354611621.6238.454.9220066511
0.1995262354610.316227763891919.072.046.0016063352
0.3162277638910.501187205315648.227.841.576728621
0.5011872053150.794328212738419.412.634.5498943037
0.7943282127381.25892543793286.24.829.696050389

Exercises

  • do the same thing for the source 1ES 2322-40.9 (faint BL Lac object)
  • repeat the procedure 10 times and average detection results (excess and significance)
  • estimate the time needed to have a 5-sigma detection for Cen A (core)

What next?

In this tutorial we learned how to simulate the expected counts of a Fermi/LAT source in the CTA energy range. Here’s some suggestions where to go next: