{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to combine observations?" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "
\"Download
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "We will start with the usual Python imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gammalib\n", "import ctools\n", "import cscripts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulating the dataset\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# First pointing\n", "obssim = ctools.ctobssim()\n", "obssim['ra'] = 83.63\n", "obssim['dec'] = 21.51\n", "obssim['rad'] = 5.0\n", "obssim['tmin'] = 0\n", "obssim['tmax'] = 1800\n", "obssim['emin'] = 0.1\n", "obssim['emax'] = 100.0\n", "obssim['caldb'] = 'prod2'\n", "obssim['irf'] = 'South_0.5h'\n", "obssim['inmodel'] = '$CTOOLS/share/models/crab.xml'\n", "obssim['outevents'] = 'events1.fits'\n", "obssim.execute()\n", "\n", "# Second pointing\n", "obssim = ctools.ctobssim()\n", "obssim['ra'] = 83.63\n", "obssim['dec'] = 22.51\n", "obssim['rad'] = 5.0\n", "obssim['tmin'] = 1800\n", "obssim['tmax'] = 3600\n", "obssim['emin'] = 0.1\n", "obssim['emax'] = 100.0\n", "obssim['caldb'] = 'prod2'\n", "obssim['irf'] = 'South_0.5h'\n", "obssim['inmodel'] = '$CTOOLS/share/models/crab.xml'\n", "obssim['outevents'] = 'events2.fits'\n", "obssim.execute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will produce the two event files `events1.fits` and `events2.fits` on disk.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "xml = gammalib.GXml()\n", "obslist = xml.append('observation_list title=\"observation library\"')\n", "for i in range(2):\n", " obs = obslist.append('observation name=\"Crab\" id=\"%02d\" instrument=\"CTA\"' % (i+1))\n", " obs.append('parameter name=\"EventList\" file=\"events%d.fits\"' % (i+1))\n", "xml.save('obs.xml')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's peek at the observations definition file." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GXml ===\n", "GXmlDocument::version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"\n", "GXmlElement::observation_list title=\"observation library\"\n", " GXmlElement::observation name=\"Crab\" id=\"01\" instrument=\"CTA\"\n", " GXmlElement::parameter name=\"EventList\" file=\"events1.fits\"\n", " GXmlElement::observation name=\"Crab\" id=\"02\" instrument=\"CTA\"\n", " GXmlElement::parameter name=\"EventList\" file=\"events2.fits\"\n" ] } ], "source": [ "print(gammalib.GXml('obs.xml'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "**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.\n", "\n", "## Combined likelihood analysis\n", "\n", "Now you are ready to do a joint maximum likelihood analysis using `ctlike`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "like = ctools.ctlike()\n", "like['inobs'] = 'obs.xml'\n", "like['caldb'] = 'prod2'\n", "like['irf'] = 'South_0.5h'\n", "like['inmodel'] = '$CTOOLS/share/models/crab.xml'\n", "like['outmodel'] = 'crab_results.xml'\n", "like.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GObservations ===\n", " Number of observations ....: 2\n", " Number of models ..........: 2\n", " Number of observed events .: 48094\n", " Number of predicted events : 48093.9954936259\n" ] } ], "source": [ "print(like.obs())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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**." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GOptimizerLM ===\n", " Optimized function value ..: 320268.838\n", " Absolute precision ........: 0.005\n", " Acceptable value decrease .: 2\n", " Optimization status .......: converged\n", " Number of parameters ......: 10\n", " Number of free parameters .: 4\n", " Number of iterations ......: 2\n", " Lambda ....................: 1e-05\n" ] } ], "source": [ "print(like.opt())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**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. \n", "\n", "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.\n", "\n", "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.\n", "\n", "In principle, unbinned and binned observations may also be combined in a joint analysis, although this Use Case may be a bit academic." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 2 }