How to combine observations?

Generally, the CTA data you may want to analyse will not only be composed of a single observation (a.k.a. run) but of a list of observations that should be combined in a joint analysis. ctools has the capability to collect individual observations in a list and to perform for example a joint maximum likelihood fit of all observations in a single shot. Here is an example that illustrates how to do that.

We will start with the usual Python imports

import gammalib
import ctools
import cscripts

Simulating the dataset

Let’s start with the simulation of two 30 min long observations of the Crab nebula, each offset by 0.5° from the nebula in opposite directions:

# First pointing
obssim = ctools.ctobssim()
obssim['ra']        = 83.63
obssim['dec']       = 21.51
obssim['rad']       = 5.0
obssim['tmin']      = 0
obssim['tmax']      = 1800
obssim['emin']      = 0.1
obssim['emax']      = 100.0
obssim['caldb']     = 'prod2'
obssim['irf']       = 'South_0.5h'
obssim['inmodel']   = '$CTOOLS/share/models/crab.xml'
obssim['outevents'] = 'events1.fits'
obssim.execute()

# Second pointing
obssim = ctools.ctobssim()
obssim['ra']        = 83.63
obssim['dec']       = 22.51
obssim['rad']       = 5.0
obssim['tmin']      = 1800
obssim['tmax']      = 3600
obssim['emin']      = 0.1
obssim['emax']      = 100.0
obssim['caldb']     = 'prod2'
obssim['irf']       = 'South_0.5h'
obssim['inmodel']   = '$CTOOLS/share/models/crab.xml'
obssim['outevents'] = 'events2.fits'
obssim.execute()

This will produce the two event files events1.fits and events2.fits on disk.

As next step you have to create an observation definition XML file that collects both observations in a list. We’ll do it in Python, but for a simple observation definition file like this you could also create it using a text editor.

xml     = gammalib.GXml()
obslist = xml.append('observation_list title="observation library"')
for i in range(2):
    obs = obslist.append('observation name="Crab" id="%02d" instrument="CTA"' % (i+1))
    obs.append('parameter name="EventList" file="events%d.fits"' % (i+1))
xml.save('obs.xml')

Let’s peek at the observations definition file.

print(gammalib.GXml('obs.xml'))

=== GXml ===
GXmlDocument::version="1.0" encoding="UTF-8" standalone="no"
GXmlElement::observation_list title="observation library"
  GXmlElement::observation name="Crab" id="01" instrument="CTA"
    GXmlElement::parameter name="EventList" file="events1.fits"
  GXmlElement::observation name="Crab" id="02" instrument="CTA"
    GXmlElement::parameter name="EventList" file="events2.fits"

The file contains a single observation_list tag that contains two observation tags that each define an observation. Each observation has a name, an id and an instrument attribute. The name attribute can have any arbitrary value, and may be the same for all observations. However, the id attribute needs to be a unique character string for any given instrument. The instrument attribute is a case-sensitive string that identifies the instrument with which the observation was taken. Please make sure that the instrument string is set correctly so that ctools knows which instrument specific functions need to be called.

Note: the instrument string for a CTA observation is unsurprisingly CTA. In case that you want to analyse data from an existing Imaging Air Cherenkov Telescope you can also set the instrument string to HESS, MAGIC, or VERITAS. You may also combine observations from different telescopes for a joint analysis in an observation definition file. Please recall that instrument strings are case sensitive.

Combined likelihood analysis

Now you are ready to do a joint maximum likelihood analysis using ctlike.

like = ctools.ctlike()
like['inobs']    = 'obs.xml'
like['caldb']    = 'prod2'
like['irf']      = 'South_0.5h'
like['inmodel']  = '$CTOOLS/share/models/crab.xml'
like['outmodel'] = 'crab_results.xml'
like.run()

Instead of providing an event list or a counts cube, you now provide the filename of the observation definition XML file (here obs.xml) as input parameter. ctlike recognises this format and automatically performs a joint maximum likelihood analysis.

print(like.obs())

=== GObservations ===
 Number of observations ....: 2
 Number of models ..........: 2
 Number of observed events .: 48094
 Number of predicted events : 48093.9954936259

The fit converged quickly, the spectral parameters of the Crab nebula have now been constrained using the events from both observations. The computation time increases roughly linearly with the number of observations that are combined, although ctools implements parallel multi-core processing which will spread the likelihood computation for the different observations over all CPU cores that are available. Doing a joint unbinned analysis is thus an efficient solution if data from multiple observations should be combined.

print(like.opt())

=== GOptimizerLM ===
 Optimized function value ..: 320268.838
 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

Note: since you only invoked the run() method the output XML file was not saved on disk. You can save it now by invoking the save() method, or save directly along with running by invoking the execute() method.

Combining observations is not limited to unbinned data (i.e. event lists) but may also be applied to binned data (i.e. counts cubes). Using ctbin you can create counts cubes from both event lists which may then be combined in an observation definition XML file, that has CountsCube rather than EventList parameters describing each observation.

Feeding the observation definition XML file to ctlike will then lead to a joint binned analysis. In the joint binned analysis, the events of individual observations are not combined, but are kept separate in distinct counts cubes. This is not very efficient, as generally counts cubes for short duration observations are only sparsly populated and the likelihood computation has to loop over a hugh number of data space bins (though also here ctlike benefits from multi-core parallel processing). Though possible, a joint binned analysis is thus not the recommended method for combining observations. An alternative is to stack the events of all observations into a single counts cube. The following section describes how such a stacked analysis is done with ctools.

In principle, unbinned and binned observations may also be combined in a joint analysis, although this Use Case may be a bit academic.