In this notebook you will learn how to use the ctools and cscripts from Python instead of typing the commands in the console.
ctools provides two Python modules that allow using all tools and scripts as Python classes. To use ctools from Python you have to import the ctools
and cscripts
modules into Python. You should also import the gammalib
module, as ctools without GammaLib is generally not very useful.
Warning: Always import gammalib
before you import ctools
and cscripts
.
In [1]:
import gammalib
import ctools
import cscripts
In [2]:
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.5
sim['dec'] = 22.8
sim['rad'] = 5.0
sim['tmin'] = '2020-01-01T00:00:00'
sim['tmax'] = '2020-01-01T01:00:00'
sim['emin'] = 0.03
sim['emax'] = 150.0
sim.execute()
The first line generates an instance of the ctobssim tool as a Python class. User parameters are then set using the [ ]
operator. After setting all parameters the execute()
method is called to execute the ctobssim
tool. On output the events.fits
FITS file is created. Until now everything is analogous to running the tool from the command line, but in Python you can easily combine the different blocks into more complex workflows.
Remember that you can consult the manual of each tool to find out how it works and to discover all the input parameters that you can set.
In [3]:
!ctobssim --help
In a Jupyter notebook a code line starting with !
is executed in the shell, so you can do the operation above just from the command line.
One of the advantages of using ctools from Python is that you can run a tool using
In [4]:
sim.run()
The main difference to the execute()
method is that the run()
method will not write the output (i.e., the simulated event list) 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. You can always save to disk later using the save()
method.
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:
In [5]:
print(sim.obs())
There is one CTA observation in the container and to visualise that observation type:
In [6]:
print(sim.obs()[0])
The observation contains a CTA event list that is implement by the GammaLib class GCTAEventList
. You can access the event list using the events()
method. 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]
, the energies and the terrestrial times (TT) of all events. Let's peek at the first events of the list:
In [7]:
events = sim.obs()[0].events()
for event in events[:1]:
print(event)
We can use this feature to inspect some of the event properties, for example look at their energy spectrum. For this we will use the Python packages matplotlib. If you do not have matplotlib you can use another plotting package of your choice or skip this step.
In [8]:
import matplotlib.pyplot as plt
%matplotlib inline
#this will visualize plots inline
ax = plt.subplot()
ax.set_yscale('log')
ax.set_xlabel('Log10(Energy [TeV])')
energies = []
for event in events:
energies.append(event.energy().log10TeV())
n, bins, patches = plt.hist(energies)
In [9]:
like = ctools.ctlike(sim.obs())
like.run()
This is very compact. Where do we define the model fit to the data? Where are the user parameters? ctlike
doesn’t in fact need any parameters as all the relevant information is already contained in the observation container produced by the ctobssim class. Indeed, you constructed the ctlike instance by using the ctobssim observation container as constructor argument.
An observation container, implemented by the GObservations
class of GammaLib, is the fundamental brick of any ctools analysis. An observation container can hold more than events, e.g., in this case it also holds the model that was used to generate the events.
Many tools and scripts handle observation containers, and accept them upon construction and return them after running the tool via the obs()
method. 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). Also, if the script crashes the information is lost.
To check how the fit went you can inspect the optimiser used by ctlike
by typing:
In [10]:
print(like.opt())
You see that 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 can be access using the models()
method of the observation container:
In [11]:
print(like.obs().models())
For example, in this way we can fetch the minimum of the optimized function (the negative of the natural logarithm of the likelihood) to compare different model hypotheses.
In [12]:
like1 = like.opt().value()
print(like1)
Suppose you want to repeat the fit by optimising also the position of the point source. This is easy from Python, as you can modify the model and fit interactively. Type the following:
In [13]:
like.obs().models()['Crab']['RA'].free()
like.obs().models()['Crab']['DEC'].free()
like.run()
print(like.obs().models())
Now we can quantify the improvement of the model by comparing the new value of the optimized function with the previous one.
In [14]:
like2 = like.opt().value()
ts = -2.0 * (like2 - like1)
print(ts)
The test statistic (TS) is expected to be distributed as a $\chi^2_n$ with a number of degrees of freedom $n$ equal to the additional number of degrees of freedom of the (second) more complex model with respect to the (first) more parsimonious one, for our case two degrees of freedom. The chance probability that the likelihood improved that much because of pure statistical fluctuations is the integral from TS to infinity of the chi square with two degrees of freedom. In this case the improvement is negligible (i.e., the chance probability is very high), as expected since in the model the source is already at the true position.
By default, tools and scripts run from Python will not generate a log file. The reason for this is that Python scripts are often used to build ctools analysis pipelines and workflows, and one generally does not want that such a script pollutes the workspace with log files. You can however instruct a ctool or cscript to generate a log file as follows:
In [15]:
like.logFileOpen()
like.run()
This creates a log file named ctlike.log
in your local work space that you can visualise using any editor.