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 predicted events : 0
=== GCTAObservation ===
Name ......................:
Identifier ................:
Instrument ................: CTA
Event file ................: events.fits
Event type ................: EventList
Statistics ................: Poisson
Ontime ....................: 1800 s
Livetime ..................: 1710 s
Deadtime correction .......: 0.95
User energy range .........: undefined
=== GCTAPointing ===
Pointing direction ........: (RA,Dec)=(83.63,22.01)
=== GCTAResponseIrf ===
Caldb mission .............: cta
Caldb instrument ..........: prod2
Response name .............: South_0.5h
Energy dispersion .........: Not used
Save energy range .........: undefined
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 ..........: 45370
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.5743, DEC=21.9746 [-0.000617357,-0.000902194] Energy=328.153 GeV Time=-3.15576e+08 s (TT)
Dir=RA=83.5756, DEC=22.0127 [4.78931e-05,-0.000879886] Energy=353.264 GeV Time=-3.15576e+08 s (TT)
Dir=RA=83.6811, DEC=22.0034 [-0.000114612,0.000826825] Energy=550.35 GeV Time=-3.15576e+08 s (TT)
...
Note
If you inspect the event list in detail you may recognise that all events appear twice in the event list. This is because we ran ctobssim twice (a first time when calling the execute() method and a second time when calling the run() method). ctobssim will append simulated events to any pre-existing events in an observation container. To get a single simulation you should clear the observation container and re-run ctobssim:
>>> sim.obs().clear()
>>> sim.run()
>>> print(sim.obs()[0].events())
=== GCTAEventList ===
Number of events ..........: 22685
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
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 ..: 152299.205
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 ..........: "SkyDirFunction" * "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.81591e-16 +/- 1.01725e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
Index ....................: -2.50618 +/- 0.0153538 [-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 ................: 0.998788 +/- 0.0118567 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)
Index ....................: 0.00990034 +/- 0.00735768 [-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.6334 +/- 0.000917669 [-360,360] deg (free,scale=1)
DEC ......................: 22.0144 +/- 0.000836554 [-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 predicted events : 22685
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 predicted events : 0
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 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 degrees. The standard dummy CTA calibration information 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.
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)
models["Crab"]["Prefactor"].fix()
models["Crab"]["Prefactor"].free()
if models["Crab"]["Prefactor"].is_free(): ...
models["Crab"]["Index"].value(-2.5)
models["Crab"]["Index"].min(-5.0) models["Crab"]["Index"].max(-1.0)or quicker
models["Crab"]["Index"].range(-5.0,-1.0)
for model in models: # Loop over models for par in model: # Loop over parameters par.fix() # fix parameter
models.save("my_model.xml")
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 file in that folder for an explanation of the scripts.