Table Of Contents

Previous topic

CTA Instrument Response Functions

Next topic

Analysing current IACT data

Using ctools from Python

Getting started

ctools provides a Python module that allows using all tools and scripts as Python classes. To use ctools from Python all you have to do is to import the ctools module into Python. You should also import the GammaLib module, as ctools without GammaLib is generally not very useful.

>>> import gammalib
>>> import ctools

To illustrate how to use ctools from Python, below is a working example of an event list simulation using the ctobssim class. An instance sim of the ctobssim class is generated and user parameters are set for this instance using the [ ] operator. See the Reference Manual for a list of parameters and their types). The execute() method executes the ctobssim class in the same way as if it were executed from the command line.

>>> import ctools
>>> sim = ctools.ctobssim()
>>> sim["inmodel"] = "${CTOOLS}/share/models/crab.xml"
>>> sim["outevents"] = "events.fits"
>>> sim["caldb"] = "prod2"
>>> sim["irf"] = "South_0.5h"
>>> sim["ra"] = 83.63
>>> sim["dec"] = 22.01
>>> sim["rad"] = 5.0
>>> sim["tmin"] = 0.0
>>> sim["tmax"] = 1800.0
>>> sim["emin"] = 0.1
>>> sim["emax"] = 100.0
>>> sim.execute()

Alternatively, you may “run” the ctobssim tool using

>>> sim.run()

The main difference to the execute() method is that the run() will not write the simulated event file to disk. Why is this useful? Well, after having typed sim.run() the ctobssim class still exists as an object in memory, including all the simulated events. The ctobssim class has an obs() method that returns an observation container that holds the simulated CTA observation with its associated events. To visualise this container, type:

>>> print(sim.obs())
=== GObservations ===
 Number of observations ....: 1
 Number of models ..........: 2
 Number of observed events .: 23099
 Number of predicted events : 0

There is one CTA observation in the container and to visualise the events in that observation you may type:

>>> print(sim.obs()[0].events())
=== GCTAEventList ===
 Number of events ..........: 23099 (disposed in "events.fits")
 Time interval .............: 51544.5 - 51544.5 days
=== GEbounds ===
 Number of intervals .......: 1
 Energy range ..............: 100 GeV - 100 TeV
=== GCTARoi ===
 ROI centre ................: RA=83.63, DEC=22.01 [0,0]
 ROI radius ................: 5 deg

The obs()[0] operator returns the first observation in the observation container, the events() operator returns the event list in that observation. To see what kind of object you actually got, use:

>>> type(sim.obs()[0].events())
<class 'gammalib.cta.GCTAEventList'>

The CTA event list is implement by the GCTAEventList class in the cta module of GammaLib. To visualise the individual events you can iterate over the events using a for loop. This will show the simulated celestial coordinates (RA, DEC), the coordinate in the camera system [DETX, DETY], energies and terrestrial times (TT) of all events.

>>> events = sim.obs()[0].events()
>>> for event in events:
...     print(event)
...
Dir=RA=83.6308, DEC=21.8881 [-0.00212759,1.33661e-05] Energy=106.465 GeV Time=-3.15576e+08 s (TT)
Dir=RA=83.7518, DEC=21.8064 [-0.0035525,0.00197398] Energy=117.706 GeV Time=-3.15576e+08 s (TT)
Dir=RA=83.5545, DEC=22.0933 [0.00145377,-0.00122121] Energy=138.624 GeV Time=-3.15576e+08 s (TT)
...

Performing a maximum likelihood analysis in Python

We can now benefit from the fact that we have some simulated events in memory to fit a model to these events using the ctlike class. We will do this in unbinned mode. Here is what you have to do:

>>> like = ctools.ctlike(sim.obs())
>>> like.run()

This is pretty compact. Where are the user parameters? ctlike doesn’t in fact need any as all the relevant information is already contained in the observation container produced by the ctobssim class. And you may have recognised that we constructed the ctlike instance by using the ctobssim observation container as constructor argument.

To check how the fit went you can inspect the optimiser used by ctlike:

>>> print(like.opt())
=== GOptimizerLM ===
 Optimized function value ..: 154553.422
 Absolute precision ........: 0.005
 Acceptable value decrease .: 2
 Optimization status .......: converged
 Number of parameters ......: 10
 Number of free parameters .: 4
 Number of iterations ......: 2
 Lambda ....................: 1e-05

Apparently, the fit converged after 2 iterations. Out of 10 parameters in the model 4 have been fitted (the others were kept fixed). To inspect the fit results you can print the model container that is a member of the observation container:

>>> print(like.obs().models())
=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 10
=== GModelSky ===
 Name ......................: Crab
 Instruments ...............: all
 Instrument scale factors ..: unity
 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.82698e-16 +/- 1.02186e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
  Index ....................: -2.47534 +/- 0.0154764 [-0,-5]  (free,scale=-1,gradient)
  PivotEnergy ..............: 300000 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
=== GCTAModelIrfBackground ===
 Name ......................: CTABackgroundModel
 Instruments ...............: CTA
 Instrument scale factors ..: unity
 Observation identifiers ...: all
 Model type ................: "PowerLaw" * "Constant"
 Number of parameters ......: 4
 Number of spectral par's ..: 3
  Prefactor ................: 1.01266 +/- 0.0119676 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
  Index ....................: 0.00474762 +/- 0.00731725 [-5,5]  (free,scale=1,gradient)
  PivotEnergy ..............: 1e+06 [10000,1e+09] MeV (fixed,scale=1e+06,gradient)
 Number of temporal par's ..: 1
  Normalization ............: 1 (relative value) (fixed,scale=1,gradient)

Suppose you want to repeat the fit by optimising also the position of the point source. This is easy from Python:

>>> like.obs().models()["Crab"]["RA"].free()
>>> like.obs().models()["Crab"]["DEC"].free()
>>> like.run()
>>> print(like.obs().models())
...
  RA .......................: 83.6327 +/- 0.000916983 [-360,360] deg (free,scale=1)
  DEC ......................: 22.0141 +/- 0.00086378 [-90,90] deg (free,scale=1)

The like.obs().models() method provides the model container, using the ["Crab"] operator we access the Crab model in that container and using the ["RA"] and ["DEC"] methods we access the relevant model parameters. The free() method frees a parameter, the opposite would be a call to the fix() method.

Warning

Passing observation containers between ctools classes is a very convenient and powerful way of building in-memory analysis pipelines. However, this implies that you need some computing ressources when dealing with large observation containers (for example if you want to analyse a few 100 hours of data at once).

Warning

You have to be aware about the scope of the objects you’re dealing with. In the above example, the sim.obs() container is allocated by the ctobssim class, hence it disappears (a.k.a. goes out of scope) once the ctobssim class is deleted, as illustrated by the following example:

>>> obs = sim.obs()
>>> del sim
>>> print(obs)
Segmentation fault

Note that

>>> obs = sim.obs()
>>> del sim
>>> print(like.obs())
=== GObservations ===
 Number of observations ....: 1
 Number of models ..........: 2
 Number of observed events .: 23099
 Number of predicted events : 23099

is okay as the ctools.ctlike(sim.obs()) constructor will create a copy of the observation container that lives within the ctlike instance. To preserve an observation container after a ctools object goes out of scope you have to create a local copy of the container using the copy() method:

>>> obs = sim.obs().copy()
>>> del sim
>>> print(obs)
=== GObservations ===
 Number of observations ....: 1
 Number of models ..........: 2
 Number of observed events .: 23099
 Number of predicted events : 0

Using obsutils

ctools provides the Python module obsutils that may further simplify your analysis efforts. obsutils is a Python script that makes use of the GammaLib and ctools modules to create standard analysis steps. As all Python scripts, obsutils is part of the cscripts module that is imported using

>>> import cscripts

Here an example of how to use obsutils:

>>> import gammalib
>>> import ctools
>>> from cscripts import obsutils
>>> pattern = obsutils.set_obs_patterns("four", ra=83.63, dec=22.01, offset=1.0)
>>> obs = obsutils.set_obs_list(pattern, duration=1800, emin=0.1, emax=100.0, rad=5.0, caldb="prod2", irf="South_0.5h")
>>> print(obs)
=== GObservations ===
 Number of observations ....: 4
 Number of models ..........: 0
 Number of observed events .: 0
 Number of predicted events : 0
>>> obs.models(gammalib.GModels("${CTOOLS}/share/models/crab.xml"))
>>> obs = obsutils.sim(obs)
>>> like = obsutils.fit(obs)
>>> print(like.obs().models())
=== GModels ===
 Number of models ..........: 2
 Number of parameters ......: 10
...

The module is imported using the from cscripts import obsutils directive. The obsutils.set_obs_patterns() function will create a pointing pattern of four observations located at offset angles of 1 degree from the nominal location of the Crab nebula. The obsutils.set_obs_list() will build an observation container from that pattern where each pointing will have a duration of 1800 seconds, cover the 0.1-100 TeV energy range and a field of view of 5°. The South_0.5h IRF from the Prod2 calibration database will be used. A model is then appended to the observation container using the obs.models() method. The obsutils.sim() function then simulates the event data, the obsutils.fit() function performs a maximum likelihood fit.

Note

The obsutils module is not yet fully developed and more convenience functions will be added in the future.

Access analysis results

Here are some examples that show how to access your analysis results in python. In the following it is assumed that you have a ctlike object called like which is setup and runs:

import gammalib
import ctools
like = ctools.ctlike(sim.obs())
like.run()
  • Best-fit parameters:

    The following command shows how to access the fit parameters and errors:

    obs = like.obs() # Get observations object
    obs.models()["Crab"]["Prefactor"].value() # This returns the actual fit value
    obs.models()["Crab"]["Prefactor"].error() # This returns the actual fit error
    
  • Open XML model file from python and print models on screen:

    models = gammalib.GModels("$CTOOLS/share/models/crab.xml")
    print(models)
    
  • Likelihood value

    like.opt().value() # Returns likelihood value
    
  • Curvature Matrix (aka Hessian)

    curvature = like.opt().curvature().invert() # Return GMatrix object
    print(curvature)
    

    To get the covariance matrix, the curvature matrix needs to be inverted:

    covariance = curvature.invert()
    print(covariance)
    

Modify and work on XML models

  • Fix a certain parameter
models["Crab"]["Prefactor"].fix()
  • Release a certain parameter
models["Crab"]["Prefactor"].free()
  • Check if parameter is free
if models["Crab"]["Prefactor"].is_free():
    ...
  • Set parameter value
models["Crab"]["Index"].value(-2.5)
  • Set parameter boundaries
models["Crab"]["Index"].min(-5.0)
models["Crab"]["Index"].max(-1.0)

or quicker

models["Crab"]["Index"].range(-5.0,-1.0)
  • Loop over models and parameters (e.g. fix all parameters)
for model in models: # Loop over models
    for par in model: # Loop over parameters
        par.fix() # fix parameter
  • Save modified XML model
models.save("my_model.xml")

Beyond the first steps

You now have learned the basics of using ctools and GammaLib within Python. To go beyond these initial steps you may check the Python scripts in the examples folder that provide useful analysis examples. Check the README.md file in that folder for an explanation of the scripts.