How to handle models from Python?

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.

import gammalib

Reading and updating an existing XML model

Opening an XML model and parsing its content

An existing XML model can be read into a model container (GModels class).

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.

print(container)

=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 10
=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: PointSource
 Model components ..........: "PointSource" * "PowerLaw" * "Constant"
 Number of parameters ......: 6
 Number of spatial par's ...: 2
  RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
  DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.48 +/- 0 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0
=== GCTAModelIrfBackground ===
 Name ......................: CTABackgroundModel
 Instruments ...............: CTA
 Observation identifiers ...: all
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 1 +/- 0 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: 0 +/- 0 [-5,5]  (free,scale=1,gradient)
  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

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.

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

value 5.7e-16
error 0.0
min value 1e-23
max value 1e-13
is free? True

You can access specific model components through the spatial, spectral, and temporal methods, for example:

print(container['Crab'].spectral())

=== GModelSpectralPlaw ===
 Number of parameters ......: 3
  Prefactor ................: 5.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.48 +/- 0 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)

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.

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

###### Name: Crab
Prefactor = 5.7e-16
Index = -2.48
PivotEnergy = 300000.0
###### Name: CTABackgroundModel
Prefactor = 1.0
Index = 0.0
PivotEnergy = 1000000.0

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

# 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')

Crab photon flux: 1.8803968833520812e-11 cm-2 s-1
Crab energy flux: 6.426073324618545e-11 erg cm-2 s-1

# 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)')

Crab flux in region: 1.0 (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.

Updating a model

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.

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

=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: PointSource
 Model components ..........: "PointSource" * "PowerLaw" * "Constant"
 Number of parameters ......: 6
 Number of spatial par's ...: 2
  RA .......................: 83.7331 +/- 0 [79,89] deg (free,scale=1)
  DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 5.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.48 +/- 0 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0

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.

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)

=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: PointSource
 Model components ..........: "PointSource" * "ExponentialCutoffPowerLaw" * "Constant"
 Number of parameters ......: 7
 Number of spatial par's ...: 2
  RA .......................: 83.7331 +/- 0 [79,89] deg (free,scale=1)
  DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)
 Number of spectral par's ..: 4
  Prefactor ................: 5.7e-16 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: -2.48 +/- 0 [-10,10]  (free,scale=1,gradient)
  CutoffEnergy .............: 50000000 +/- 0 [0.1,infty[ MeV (free,scale=1,gradient)
  PivotEnergy ..............: 300000 MeV (fixed,scale=1,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0

In the last example we remove the IRF background model, and append a cube background model from another container.

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

=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 11
=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: PointSource
 Model components ..........: "PointSource" * "ExponentialCutoffPowerLaw" * "Constant"
 Number of parameters ......: 7
 Number of spatial par's ...: 2
  RA .......................: 83.7331 +/- 0 [79,89] deg (free,scale=1)
  DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)
 Number of spectral par's ..: 4
  Prefactor ................: 5.7e-16 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: -2.48 +/- 0 [-10,10]  (free,scale=1,gradient)
  CutoffEnergy .............: 50000000 +/- 0 [0.1,infty[ MeV (free,scale=1,gradient)
  PivotEnergy ..............: 300000 MeV (fixed,scale=1,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0
=== GCTAModelCubeBackground ===
 Name ......................: Background
 Instruments ...............: CTA
 Observation identifiers ...: all
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 1 +/- 0 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: 0 +/- 0 [-5,5]  (free,scale=1,gradient)
  PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

Save model changes to an XML file

You can use directly the model container in Python to pass it to ctools or cscripts, but you can also write it to an XML file.

container.save('my_crab.xml')

Creating a new model

New models can be created directly in Python. First, create a new empty model container.

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.

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

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

# create source
source = gammalib.GModelSky(spatial,spectral)
source.name('my_source')

# append to container
my_container.append(source)

<gammalib.model.GModelSky; proxy of <Swig Object of type 'GModelSky *' at 0x11252dae0> >

We also add a CTA IRF background model with power-law spectral correction.

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

<gammalib.cta.GCTAModelIrfBackground; proxy of <Swig Object of type 'GCTAModelIrfBackground *' at 0x11252d960> >

We inspect the model we just created.

print(my_container)

=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 11
=== GModelSky ===
 Name ......................: my_source
 Instruments ...............: all
 Observation identifiers ...: all
 Model type ................: ExtendedSource
 Model components ..........: "RadialGaussian" * "PowerLaw" * "Constant"
 Number of parameters ......: 7
 Number of spatial par's ...: 3
  RA .......................: 54 deg (fixed,scale=1)
  DEC ......................: -19 deg (fixed,scale=1)
  Sigma ....................: 0.5 +/- 0 [0.0002778,infty[ deg (free,scale=1)
 Number of spectral par's ..: 3
  Prefactor ................: 1e-18 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1e-18,gradient)
  Index ....................: -2.5 +/- 0 [10,-10]  (free,scale=-2.5,gradient)
  PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
 Number of scale par's .....: 0
=== GCTAModelIrfBackground ===
 Name ......................: Background
 Instruments ...............: CTA
 Observation identifiers ...: all
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 1 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: 0 +/- 0 [-10,10]  (free,scale=1,gradient)
  PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

The container can be used in Python or written to disk in an XML file.

my_container.save('mymodel.xml')