{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# How to handle models from Python?" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this tutorial you will learn how to handle models from Python. This is done using the GammaLib classes `GModels`, `GModel`, `GModelSpatial`, `GModelSpectral`, and `GModelTemporal`. You can find all the information on the GammaLib classes in the [Doxygen documentation](http://cta.irap.omp.eu/gammalib-devel/doxygen/index.html).\n", "\n", "In GammaLib there are some conventions about units. If not otherwise specified:\n", "\n", "- energies are in MeV\n", "\n", "- photon fluxes are in photons/cm2/s\n", "\n", "- differential photon fluxes are in photons/cm2/s/MeV\n", "\n", "- energy fluxes are in erg/cm2/s\n", "\n", "To start we import the gammalib Python module." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import gammalib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reading and updating an existing XML model\n", "\n", "### Opening an XML model and parsing its content\n", "\n", "An existing XML model can be read into a model container (`GModels` class)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "container = gammalib.GModels('$CTOOLS/share/models/crab.xml')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The easiest way to inspect a model container is by using the function `print` that will display in the terminal its entire content." ] }, { "cell_type": "code", "execution_count": 3, "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.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n", " Index ....................: -2.48 +/- 0 [-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", "=== GCTAModelIrfBackground ===\n", " Name ......................: CTABackgroundModel\n", " Instruments ...............: CTA\n", " Observation identifiers ...: all\n", " Model type ................: \"PowerLaw\" * \"Constant\"\n", " Number of parameters ......: 4\n", " Number of spectral par's ..: 3\n", " Prefactor ................: 1 +/- 0 [0.001,1000] ph/cm2/s/MeV (free,scale=1,gradient)\n", " Index ....................: 0 +/- 0 [-5,5] (free,scale=1,gradient)\n", " PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n", " Number of temporal par's ..: 1\n", " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n" ] } ], "source": [ "print(container)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The model container contains two models, for a total of ten parameters. The first model is a sky model named `Crab`, the second a background model. Each model has a spatial, a spectral, and a temporal component. Each model component has parameters. For each parameter you have a value, an error (in our example it is always zero because no fit to the data was done yet), an allowed value range for fitting the parameter to data, and, whenever relevant, units. In the parenthesis you can se if the parameter is free or fixed, and some more technical information, i.e, if the model has an analytical gradient w.r.t. this parameter implemented (can be fit to data quickly under certain conditions), and if internally the parameter is scaled (parameters should be scaled so that value/scale is of order unity for a more efficient handling by the optimiser during a fit to the data). Furthemore each model can be restricted to a specific instrument or observation. In the example the only restriction is that the CTA background model is used only for CTA data (`Instruments ...............: CTA`).\n", "\n", "Model containers can be parsed using the model and parameter names if they are know to you. Below for example we directly access the properties of the Prefactor for the spectral model of the Crab source." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "value 5.7e-16\n", "error 0.0\n", "min value 1e-23\n", "max value 1e-13\n", "is free? True\n" ] } ], "source": [ "print('value', container['Crab']['Prefactor'].value())\n", "print('error', container['Crab']['Prefactor'].error())\n", "print('min value', container['Crab']['Prefactor'].min())\n", "print('max value', container['Crab']['Prefactor'].max())\n", "print('is free?', container['Crab']['Prefactor'].is_free())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can access specific model components through the `spatial`, `spectral`, and `temporal` methods, for example:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GModelSpectralPlaw ===\n", " Number of parameters ......: 3\n", " Prefactor ................: 5.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n", " Index ....................: -2.48 +/- 0 [-0,-5] (free,scale=-1,gradient)\n", " PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n" ] } ], "source": [ "print(container['Crab'].spectral())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can also blindly parse models and model parameters. Below we get names and values for the spectral parameters of all the models in the container." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "###### Name: Crab\n", "Prefactor = 5.7e-16\n", "Index = -2.48\n", "PivotEnergy = 300000.0\n", "###### Name: CTABackgroundModel\n", "Prefactor = 1.0\n", "Index = 0.0\n", "PivotEnergy = 1000000.0\n" ] } ], "source": [ "for s in range(container.size()):\n", " model = container[s]\n", " print('###### Name:', model.name())\n", " for i in range(model.spectral().size()):\n", " param = model.spectral()[i]\n", " print(param.name(),'=',param.value())\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally there are convenience methods to compute for a model component the flux and energy flux over an energy interval (spectral components), and the flux within a circular region (spatial models). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Crab photon flux: 1.8803968833520812e-11 cm-2 s-1\n", "Crab energy flux: 6.426073324618545e-11 erg cm-2 s-1\n" ] } ], "source": [ "# calculate Crab spectral model flux in the 1-10 TeV energy range\n", "\n", "# define energy bounds\n", "emin = gammalib.GEnergy(1.,'TeV')\n", "emax = gammalib.GEnergy(10.,'TeV')\n", "\n", "# extract model from container'\n", "crab = container['Crab']\n", "\n", "# photon flux (cm-2 s-1)\n", "flux = crab.spectral().flux(emin,emax)\n", "print('Crab photon flux:', flux, 'cm-2 s-1')\n", "\n", "# energy flux (cm-2 s-1)\n", "eflux = crab.spectral().eflux(emin,emax)\n", "print('Crab energy flux', eflux, 'erg cm-2 s-1')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Crab flux in region: 1.0 (relative value)\n" ] } ], "source": [ "# calculate Crab spatial model flux in 0.2 deg region centred on the source\n", "\n", "# region centre = source position\n", "centre = crab.spatial().dir()\n", "\n", "# circular region\n", "reg = gammalib.GSkyRegionCircle(centre,0.2)\n", "\n", "# flux in the circular region (relative value)\n", "flux = crab.spatial().flux(reg)\n", "print('Crab flux in region:',flux,'(relative value)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Spatial models can be normalised to unity. This is always the case for analytical models including point sources. Therefore we do expect to have a spatial model flux of 1 for a region encompassing the source." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Updating a model\n", "\n", "Once you have a model container you can update its parameters, change model components and attributes, delete existing models or append new ones.\n", "\n", "In the first example we change the position of the Crab source by 0.1 deg toward positive R.A. and we free the spatial parameters, but we restrict the allowed range around the known position." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== 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.7331 +/- 0 [79,89] deg (free,scale=1)\n", " DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)\n", " Number of spectral par's ..: 3\n", " Prefactor ................: 5.7e-16 +/- 0 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)\n", " Index ....................: -2.48 +/- 0 [-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" ] } ], "source": [ "# change R.A.\n", "crab['RA'].value(crab['RA'].value() + .1)\n", "# free spatial parameters\n", "crab['RA'].free()\n", "crab['DEC'].free()\n", "# restric range\n", "crab['RA'].min(79.)\n", "crab['RA'].max(89.)\n", "crab['DEC'].min(21.5)\n", "crab['DEC'].max(22.5)\n", "\n", "# inspect model to verify changes\n", "print(crab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the second example we change the spectral model from a power law to an exponentially-cutoff power law, keeping the power law parameters to the original values, and adding a cutoff at 50 TeV." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GModelSky ===\n", " Name ......................: Crab\n", " Instruments ...............: all\n", " Observation identifiers ...: all\n", " Model type ................: PointSource\n", " Model components ..........: \"PointSource\" * \"ExponentialCutoffPowerLaw\" * \"Constant\"\n", " Number of parameters ......: 7\n", " Number of spatial par's ...: 2\n", " RA .......................: 83.7331 +/- 0 [79,89] deg (free,scale=1)\n", " DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)\n", " Number of spectral par's ..: 4\n", " Prefactor ................: 5.7e-16 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)\n", " Index ....................: -2.48 +/- 0 [-10,10] (free,scale=1,gradient)\n", " CutoffEnergy .............: 50000000 +/- 0 [0.1,infty[ MeV (free,scale=1,gradient)\n", " PivotEnergy ..............: 300000 MeV (fixed,scale=1,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" ] } ], "source": [ "expplaw = gammalib.GModelSpectralExpPlaw()\n", "expplaw['Prefactor'].value(crab['Prefactor'].value())\n", "expplaw['Index'].value(crab['Index'].value())\n", "expplaw['PivotEnergy'].value(crab['PivotEnergy'].value())\n", "expplaw['CutoffEnergy'].value(50.e6) # value in MeV\n", "crab.spectral(expplaw)\n", "\n", "# inspect model to verify changes\n", "print(crab)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the last example we remove the IRF background model, and append a cube background model from another container." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "=== GModels ===\n", " Number of models ..........: 2\n", " Number of parameters ......: 11\n", "=== GModelSky ===\n", " Name ......................: Crab\n", " Instruments ...............: all\n", " Observation identifiers ...: all\n", " Model type ................: PointSource\n", " Model components ..........: \"PointSource\" * \"ExponentialCutoffPowerLaw\" * \"Constant\"\n", " Number of parameters ......: 7\n", " Number of spatial par's ...: 2\n", " RA .......................: 83.7331 +/- 0 [79,89] deg (free,scale=1)\n", " DEC ......................: 22.0145 +/- 0 [21.5,22.5] deg (free,scale=1)\n", " Number of spectral par's ..: 4\n", " Prefactor ................: 5.7e-16 +/- 0 [0,infty[ ph/cm2/s/MeV (free,scale=1,gradient)\n", " Index ....................: -2.48 +/- 0 [-10,10] (free,scale=1,gradient)\n", " CutoffEnergy .............: 50000000 +/- 0 [0.1,infty[ MeV (free,scale=1,gradient)\n", " PivotEnergy ..............: 300000 MeV (fixed,scale=1,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 ......................: Background\n", " Instruments ...............: CTA\n", " Observation identifiers ...: all\n", " Model type ................: \"PowerLaw\" * \"Constant\"\n", " Number of parameters ......: 4\n", " Number of spectral par's ..: 3\n", " Prefactor ................: 1 +/- 0 [0.01,100] ph/cm2/s/MeV (free,scale=1,gradient)\n", " Index ....................: 0 +/- 0 [-5,5] (free,scale=1,gradient)\n", " PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n", " Number of temporal par's ..: 1\n", " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n" ] } ], "source": [ "# remove existing background model\n", "container.remove('CTABackgroundModel')\n", "\n", "# open example background cube model\n", "bkgcube_container = gammalib.GModels('$CTOOLS/share/models/bkg_cube.xml')\n", "\n", "# append first model in new container to old container\n", "container.append(bkgcube_container[0])\n", "\n", "# inspect model to verify changes\n", "print(container)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Save model changes to an XML file\n", "\n", "You can use directly the model container in Python to pass it to ctools or cscripts, but you can also write it to an XML file." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "container.save('my_crab.xml')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating a new model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "New models can be created directly in Python. First, create a new empty model container." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "my_container = gammalib.GModels()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To create a new source we need to define at least spatial and spectral components (if not specified the temporal component is taken to be constants). Refer to the [Doxygen documentation](http://cta.irap.omp.eu/gammalib-devel/doxygen/index.html) for finding out all the models available in GammaLib and how to use them.\n", "\n", "We first create a Gaussian spatial component." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "# define source direction\n", "srcdir = gammalib.GSkyDir()\n", "# set R.A. and Dec\n", "srcdir.radec_deg(54.,-19.)\n", "\n", "# Gaussian spatial component\n", "spatial = gammalib.GModelSpatialRadialGauss(srcdir,0.5) # centre and radius in deg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a power-law spectral component." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "# define pivot energy\n", "pivot = gammalib.GEnergy(1,'TeV')\n", "\n", "# power law\n", "spectral = gammalib.GModelSpectralPlaw(1.e-18,-2.5,pivot) # differential photon flux, index, pivot energy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create the source and append it to the model container." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "