{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to perform a stacked analysis?" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "
\"Download
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A stacked analysis is a binned analysis where all data from multiple observations are stacked into a single counts cube.\n", "\n", "As usual we start with Python imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gammalib\n", "import ctools\n", "import cscripts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data and observation definition XML file `obs.xml` used here are those produced in the previous tutorial [How to combine observations?](howto_combine_observations.ipynb)\n", "\n", "## Stacking Events\n", "\n", "The event stacking is done using the `ctbin` tool. Instead of providing to `ctbin` an event list you should specify the observation definition XML file on input. `ctbin` will then loop over all observations and collect all events into a single counts cube." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "evbin = ctools.ctbin()\n", "evbin['inobs'] = 'obs.xml'\n", "evbin['xref'] = 83.63 # deg; we center the cube at the position of the Crab nebula\n", "evbin['yref'] = 22.01 # deg\n", "evbin['proj'] = 'CAR'\n", "evbin['coordsys'] = 'CEL'\n", "evbin['binsz'] = 0.02 # deg/bin\n", "evbin['nxpix'] = 200 \n", "evbin['nypix'] = 200\n", "evbin['ebinalg'] = 'LOG'\n", "evbin['emin'] = 0.1 # TeV\n", "evbin['emax'] = 100.0 # TeV\n", "evbin['enumbins'] = 20\n", "evbin['outobs'] = 'cntcube.fits'\n", "evbin.execute()" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Computing the stacked response and background\n", "\n", "You now have a stacked counts cube `cntcube.fits` on disk. Before you can use that counts cube in a maximum likelihood analysis, you need to compute the stacked instrument response function and the background model that is needed for the analysis.\n", "\n", "For the instrument response function, you have to compute the total exposure for the stacked cube (i.e. the sum of the effective areas for each observation multiplied by the corresponding livetimes) and an effective point spread function (i.e. the point spread function of the different observations weighted by the corresponding exposures). Optionally, you can also compute an effective energy dispersion (i.e. the energy dispersion of the different observations weighted by the corresponding exposures). To get these informations you use the `ctexpcube`, `ctpsfcube` and `ctedispcube` tools:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "expcube = ctools.ctexpcube()\n", "expcube['inobs'] = 'obs.xml'\n", "expcube['caldb'] = 'prod2'\n", "expcube['irf'] = 'South_0.5h'\n", "expcube['incube'] = 'cntcube.fits' # exposure cube definition is copied from counts cube\n", "expcube['outcube'] = 'expcube.fits'\n", "expcube.execute()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "psfcube = ctools.ctpsfcube()\n", "psfcube['inobs'] = 'obs.xml'\n", "psfcube['caldb'] = 'prod2'\n", "psfcube['irf'] = 'South_0.5h'\n", "psfcube['incube'] = 'NONE'\n", "psfcube['xref'] = 83.63\n", "psfcube['yref'] = 22.01\n", "psfcube['proj'] = 'CAR'\n", "psfcube['coordsys'] = 'CEL'\n", "psfcube['binsz'] = 1.0 # deg/bin; the PSF only varies slowly\n", "psfcube['nxpix'] = 10\n", "psfcube['nypix'] = 10\n", "psfcube['emin'] = 0.1\n", "psfcube['emax'] = 100.0\n", "psfcube['enumbins'] = 20\n", "psfcube['outcube'] = 'psfcube.fits'\n", "psfcube.execute()" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "edispcube = ctools.ctedispcube()\n", "edispcube['inobs'] = 'obs.xml'\n", "edispcube['caldb'] = 'prod2'\n", "edispcube['irf'] = 'South_0.5h'\n", "edispcube['incube'] = 'NONE'\n", "edispcube['xref'] = 83.63\n", "edispcube['yref'] = 22.01\n", "edispcube['proj'] = 'CAR'\n", "edispcube['coordsys'] = 'CEL'\n", "edispcube['binsz'] = 1.0 # deg/bin; the energy dispersion only varies slowly\n", "edispcube['nxpix'] = 10\n", "edispcube['nypix'] = 10\n", "edispcube['emin'] = 0.1\n", "edispcube['emax'] = 100.0\n", "edispcube['enumbins'] = 20\n", "edispcube['outcube'] = 'edispcube.fits'\n", "edispcube.execute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You have noticed that for `ctexpcube` you provided an input counts cube, while for the other tools you specified NONE. By providing an input counts cube you instructed `ctexpcube` to extract the geometry of the cube from the counts cube. This is a convenient trick to reduce the number of user parameters that you need to specify. You did however not apply this trick for `ctpsfcube` and `ctedispcube`. In fact, the point spread function and energy dispersion do not vary significantly on spatial scales of 0.02°, and using the counts cube definition for these cubes would lead to large response cube files with a spatial precision that is actually not needed (the point spread function and energy dispersion cubes are actually 4-dimensional data cubes, hence their size increases quickly for a large number of spatial pixels). Therefore, you have specified a larger image scale of 1° for both cubes and only a small number of 10x10 spatial pixels, leading to point spread function and energy dispersion cubes of modest size (a few MB).\n", "\n", "You provided the `obs.xml` file that defines all observations on input so that the tools know which observations were combined in the `ctbin` run. As final step of the analysis preparation, you need to generate a background cube using the `ctbkgcube` tool." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "bkgcube = ctools.ctbkgcube()\n", "bkgcube['inobs'] = 'obs.xml'\n", "bkgcube['caldb'] = 'prod2'\n", "bkgcube['irf'] = 'South_0.5h'\n", "bkgcube['incube'] = 'cntcube.fits'\n", "bkgcube['inmodel'] = '$CTOOLS/share/models/crab.xml'\n", "bkgcube['outcube'] = 'bkgcube.fits'\n", "bkgcube['outmodel'] = 'model.xml'\n", "bkgcube.execute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The usage of `ctbkgcube` is very similar to that of `ctexpcube`, yet it takes the model definition XML file as an additional input parameter. You used here the usual `$CTOOLS/share/models/crab.xml` model file that is shipped with the ctools. `ctbkgcube` provides on output the background cube file `bkgcube.fits` and the model definition XML file `model.xml` that can be used for further analysis. Having a look at the `model.xml` file illustrates how the background modelling works:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GXml ===\n", "GXmlDocument::version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"\n", "GXmlElement::source_library title=\"source library\"\n", " GXmlElement::source name=\"Crab\" type=\"PointSource\"\n", " GXmlElement::spectrum type=\"PowerLaw\"\n", " GXmlElement::parameter name=\"Prefactor\" value=\"5.7\" error=\"0\" scale=\"1e-16\" min=\"1e-07\" max=\"1000\" free=\"1\"\n", " GXmlElement::parameter name=\"Index\" value=\"2.48\" error=\"0\" scale=\"-1\" min=\"0\" max=\"5\" free=\"1\"\n", " GXmlElement::parameter name=\"PivotEnergy\" value=\"0.3\" scale=\"1000000\" min=\"0.01\" max=\"1000\" free=\"0\"\n", " GXmlElement::spatialModel type=\"PointSource\"\n", " GXmlElement::parameter name=\"RA\" value=\"83.6331\" scale=\"1\" min=\"-360\" max=\"360\" free=\"0\"\n", " GXmlElement::parameter name=\"DEC\" value=\"22.0145\" scale=\"1\" min=\"-90\" max=\"90\" free=\"0\"\n", " GXmlElement::source name=\"BackgroundModel\" type=\"CTACubeBackground\" instrument=\"CTA,HESS,MAGIC,VERITAS\"\n", " GXmlElement::spectrum type=\"PowerLaw\"\n", " GXmlElement::parameter name=\"Prefactor\" value=\"1\" error=\"0\" scale=\"1\" min=\"0.01\" max=\"100\" free=\"1\"\n", " GXmlElement::parameter name=\"Index\" value=\"0\" error=\"0\" scale=\"1\" min=\"-5\" max=\"5\" free=\"1\"\n", " GXmlElement::parameter name=\"PivotEnergy\" value=\"1\" scale=\"1000000\" free=\"0\"\n" ] } ], "source": [ "print(gammalib.GXml('model.xml'))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Crab source component is the same that is also present in `$CTOOLS/share/models/crab.xml` and is not modified. The background component, however, has been replaced by a model of type `CTACubeBackground`. This model is a 3-dimensional data cube that describes the expected background rate as function of spatial position and energy. The data cube is multiplied by a power law spectrum that allows to adjust the normalization and slope of the background spectrum in the fit. This power law could be replaced by any spectral model that is found as an appropriate multiplicator to the background cube.\n", "\n", "**There is no constraint on providing the same spatial binning or the same energy binning for an exposure cube, a PSF cube, an energy dispersion cube, a background cube and a counts cube**. ctools interpolates internally all response cubes hence any arbitrary appropriate binning may be used. Using the same binning for the exposure cube, the background cube and the counts cube is only a convenience.\n", "\n", "## Likelihood fitting\n", "\n", "Now you have all files at hand to perform a stacked maximum likelihood analysis using the `ctlike` tool:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "like = ctools.ctlike()\n", "like['inobs'] = 'cntcube.fits'\n", "like['expcube'] = 'expcube.fits'\n", "like['psfcube'] = 'psfcube.fits'\n", "like['bkgcube'] = 'bkgcube.fits'\n", "like['inmodel'] = 'model.xml'\n", "like['outmodel'] = 'crab_results.xml'\n", "like.execute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`ctlike` uses as input observations the counts cube and therefore needs also the exposure cube, the PSF cube, and the background cube file names.\n", "\n", "The results of the `ctlike` run are shown below." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GOptimizerLM ===\n", " Optimized function value ..: 87667.270\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": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GModels ===\n", " Number of models ..........: 2\n", " Number of parameters ......: 10\n", "=== GModelSky ===\n", " Name ......................: Crab\n", " Instruments ...............: all\n", " Observation identifiers ...: all\n", " Model type ................: PointSource\n", " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n", " Number of parameters ......: 6\n", " Number of spatial par's ...: 2\n", " RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)\n", " DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)\n", " Number of spectral par's ..: 3\n", " Prefactor ................: 5.7351470584048e-16 +/- 7.15792803493375e-18 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n", " Index ....................: -2.45746587500157 +/- 0.010743485688059 [-0,-5] (free,scale=-1,gradient)\n", " PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n", " Number of temporal par's ..: 1\n", " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n", " Number of scale par's .....: 0\n", "=== GCTAModelCubeBackground ===\n", " Name ......................: BackgroundModel\n", " Instruments ...............: CTA, HESS, MAGIC, VERITAS\n", " Observation identifiers ...: all\n", " Model type ................: \"PowerLaw\" * \"Constant\"\n", " Number of parameters ......: 4\n", " Number of spectral par's ..: 3\n", " Prefactor ................: 1.02157245578124 +/- 0.0113974414365358 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)\n", " Index ....................: 0.00968334444956244 +/- 0.00675463165472791 [-5,5] (free,scale=1,gradient)\n", " PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)\n", " Number of temporal par's ..: 1\n", " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n" ] } ], "source": [ "print(like.obs().models())" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "If you want to consider also the energy dispersion during the maximum likelihood fitting you set the `edisp` parameter to `True` and provide the energy dispersion cube to the `edispcube` parameter:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "like = ctools.ctlike()\n", "like['inobs'] = 'cntcube.fits'\n", "like['expcube'] = 'expcube.fits'\n", "like['psfcube'] = 'psfcube.fits'\n", "like['bkgcube'] = 'bkgcube.fits'\n", "like['edisp'] = True # Set to True (is False by default)\n", "like['edispcube'] = 'edispcube.fits' # Provide energy dispersion cube\n", "like['inmodel'] = 'model.xml'\n", "like['outmodel'] = 'crab_results.xml'\n", "# like.execute() # Uncomment if you have some time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The maximum likelihood computation including energy dispersion is more time consuming**, and in many situations the impact of the energy dispersion on the analysis results will be very small. So make sure that you really need energy dispersion before you are using it. Uncomment the `execute()` call in the previous cell to perform the likelihood analysis including energy dispersion." ] }, { "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 }