In this tutorial you will learn how to handle models from Python. This is done using the GammaLib classes GModels, GModel, GModelSpatial, GModelSpectral, and GModelTemporal. You can find all the information on the GammaLib classes in the Doxygen documentation.
In GammaLib there are some conventions about units. If not otherwise specified:
energies are in MeV
photon fluxes are in photons/cm2/s
differential photon fluxes are in photons/cm2/s/MeV
energy fluxes are in erg/cm2/s
To start we import the gammalib Python module.
In [1]:
import gammalib
In [2]:
container = gammalib.GModels('$CTOOLS/share/models/crab.xml')
The easiest way to inspect a model container is by using the function print that will display in the terminal its entire content.
In [3]:
print(container)
The model container contains two models, for a total of ten parameters. The first model is a sky model named Crab, the second a background model. Each model has a spatial, a spectral, and a temporal component. Each model component has parameters. For each parameter you have a value, an error (in our example it is always zero because no fit to the data was done yet), an allowed value range for fitting the parameter to data, and, whenever relevant, units. In the parenthesis you can se if the parameter is free or fixed, and some more technical information, i.e, if the model has an analytical gradient w.r.t. this parameter implemented (can be fit to data quickly under certain conditions), and if internally the parameter is scaled (parameters should be scaled so that value/scale is of order unity for a more efficient handling by the optimiser during a fit to the data). Furthemore each model can be restricted to a specific instrument or observation. In the example the only restriction is that the CTA background model is used only for CTA data (Instruments ...............: CTA).
Model containers can be parsed using the model and parameter names if they are know to you. Below for example we directly access the properties of the Prefactor for the spectral model of the Crab source.
In [4]:
print('value', container['Crab']['Prefactor'].value())
print('error', container['Crab']['Prefactor'].error())
print('min value', container['Crab']['Prefactor'].min())
print('max value', container['Crab']['Prefactor'].max())
print('is free?', container['Crab']['Prefactor'].is_free())
You can access specific model components through the spatial, spectral, and temporal methods, for example:
In [5]:
print(container['Crab'].spectral())
You can also blindly parse models and model parameters. Below we get names and values for the spectral parameters of all the models in the container.
In [6]:
for s in range(container.size()):
model = container[s]
print('###### Name:', model.name())
for i in range(model.spectral().size()):
param = model.spectral()[i]
print(param.name(),'=',param.value())
Finally there are convenience methods to compute for a model component the flux and energy flux over an energy interval (spectral components), and the flux within a circular region (spatial models).
In [7]:
# calculate Crab spectral model flux in the 1-10 TeV energy range
# define energy bounds
emin = gammalib.GEnergy(1.,'TeV')
emax = gammalib.GEnergy(10.,'TeV')
# extract model from container'
crab = container['Crab']
# photon flux (cm-2 s-1)
flux = crab.spectral().flux(emin,emax)
print('Crab photon flux:', flux, 'cm-2 s-1')
# energy flux (cm-2 s-1)
eflux = crab.spectral().eflux(emin,emax)
print('Crab energy flux', eflux, 'erg cm-2 s-1')
In [8]:
# calculate Crab spatial model flux in 0.2 deg region centred on the source
# region centre = source position
centre = crab.spatial().dir()
# circular region
reg = gammalib.GSkyRegionCircle(centre,0.2)
# flux in the circular region (relative value)
flux = crab.spatial().flux(reg)
print('Crab flux in region:',flux,'(relative value)')
Spatial models can be normalised to unity. This is always the case for analytical models including point sources. Therefore we do expect to have a spatial model flux of 1 for a region encompassing the source.
Once you have a model container you can update its parameters, change model components and attributes, delete existing models or append new ones.
In the first example we change the position of the Crab source by 0.1 deg toward positive R.A. and we free the spatial parameters, but we restrict the allowed range around the known position.
In [9]:
# change R.A.
crab['RA'].value(crab['RA'].value() + .1)
# free spatial parameters
crab['RA'].free()
crab['DEC'].free()
# restric range
crab['RA'].min(79.)
crab['RA'].max(89.)
crab['DEC'].min(21.5)
crab['DEC'].max(22.5)
# inspect model to verify changes
print(crab)
In the second example we change the spectral model from a power law to an exponentially-cutoff power law, keeping the power law parameters to the original values, and adding a cutoff at 50 TeV.
In [10]:
expplaw = gammalib.GModelSpectralExpPlaw()
expplaw['Prefactor'].value(crab['Prefactor'].value())
expplaw['Index'].value(crab['Index'].value())
expplaw['PivotEnergy'].value(crab['PivotEnergy'].value())
expplaw['CutoffEnergy'].value(50.e6) # value in MeV
crab.spectral(expplaw)
# inspect model to verify changes
print(crab)
In the last example we remove the IRF background model, and append a cube background model from another container.
In [11]:
# remove existing background model
container.remove('CTABackgroundModel')
# open example background cube model
bkgcube_container = gammalib.GModels('$CTOOLS/share/models/bkg_cube.xml')
# append first model in new container to old container
container.append(bkgcube_container[0])
# inspect model to verify changes
print(container)
In [12]:
container.save('my_crab.xml')
New models can be created directly in Python. First, create a new empty model container.
In [13]:
my_container = gammalib.GModels()
To create a new source we need to define at least spatial and spectral components (if not specified the temporal component is taken to be constants). Refer to the Doxygen documentation for finding out all the models available in GammaLib and how to use them.
We first create a Gaussian spatial component.
In [14]:
# define source direction
srcdir = gammalib.GSkyDir()
# set R.A. and Dec
srcdir.radec_deg(54.,-19.)
# Gaussian spatial component
spatial = gammalib.GModelSpatialRadialGauss(srcdir,0.5) # centre and radius in deg
We create a power-law spectral component.
In [15]:
# define pivot energy
pivot = gammalib.GEnergy(1,'TeV')
# power law
spectral = gammalib.GModelSpectralPlaw(1.e-18,-2.5,pivot) # differential photon flux, index, pivot energy
We create the source and append it to the model container.
In [16]:
# create source
source = gammalib.GModelSky(spatial,spectral)
source.name('my_source')
# append to container
my_container.append(source)
Out[16]:
We also add a CTA IRF background model with power-law spectral correction.
In [17]:
# spectral correction
spectral = gammalib.GModelSpectralPlaw(1, 0, gammalib.GEnergy(1, 'TeV'))
# create background model
bkgmodel = gammalib.GCTAModelIrfBackground(spectral)
bkgmodel.name('Background')
bkgmodel.instruments('CTA')
# append to container
my_container.append(bkgmodel)
Out[17]:
We inspect the model we just created.
In [18]:
print(my_container)
The container can be used in Python or written to disk in an XML file.
In [19]:
my_container.save('mymodel.xml')