\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.10?urlpath=lab/tree/spectrum_pipe.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[spectrum_pipe.ipynb](../_static/notebooks/spectrum_pipe.ipynb) |\n",
"[spectrum_pipe.py](../_static/notebooks/spectrum_pipe.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Spectrum analysis with Gammapy (run pipeline)\n",
"\n",
"In this tutorial we will learn how to perform a 1d spectral analysis.\n",
"\n",
"We will use a \"pipeline\" or \"workflow\" class to run a standard analysis. If you're interested in implementation detail of the analysis in order to create a custom analysis class, you should read ([spectrum_analysis.ipynb](spectrum_analysis.ipynb)) that executes the analysis using lower-level classes and methods in Gammapy. \n",
"\n",
"In this tutorial we will use the folling Gammapy classes:\n",
"\n",
"- [gammapy.data.DataStore](..\/api/gammapy.data.DataStore.rst) to load the data to \n",
"- [gammapy.scripts.SpectrumAnalysisIACT](..\/api/gammapy.scripts.SpectrumAnalysisIACT.rst) to run the analysis\n",
"\n",
"We use 4 Crab observations from H.E.S.S. for testing."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"As usual, we'll start with some setup for the notebook, and import the functionality we need."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"from regions import CircleSkyRegion\n",
"\n",
"from gammapy.utils.energy import EnergyBounds\n",
"from gammapy.data import DataStore\n",
"from gammapy.scripts import SpectrumAnalysisIACT\n",
"from gammapy.catalog import SourceCatalogGammaCat\n",
"from gammapy.maps import Map\n",
"from gammapy.spectrum.models import LogParabola\n",
"from gammapy.spectrum import CrabSpectrum"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Select data\n",
"\n",
"First, we select and load some H.E.S.S. data (simulated events for now). In real life you would do something fancy here, or just use the list of observations someone send you (and hope they have done something fancy before). We'll just use the standard gammapy 4 crab runs."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[23523 23526 23559 23592]\n"
]
}
],
"source": [
"data_store = DataStore.from_dir(\"$GAMMAPY_DATA/hess-dl3-dr1/\")\n",
"mask = data_store.obs_table[\"TARGET_NAME\"] == \"Crab\"\n",
"obs_ids = data_store.obs_table[\"OBS_ID\"][mask].data\n",
"observations = data_store.get_observations(obs_ids)\n",
"print(obs_ids)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Configure the analysis\n",
"\n",
"Now we'll define the input for the spectrum analysis. It will be done the python way, i.e. by creating a config dict containing python objects. We plan to add also the convenience to configure the analysis using a plain text config file."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"crab_pos = SkyCoord.from_name(\"crab\")\n",
"on_region = CircleSkyRegion(crab_pos, 0.15 * u.deg)\n",
"\n",
"model = LogParabola(\n",
" alpha=2.3,\n",
" beta=0.01,\n",
" amplitude=1e-11 * u.Unit(\"cm-2 s-1 TeV-1\"),\n",
" reference=1 * u.TeV,\n",
")\n",
"\n",
"flux_point_binning = EnergyBounds.equal_log_spacing(0.7, 30, 5, u.TeV)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAEHCAYAAACwfMNTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAE7hJREFUeJzt3X+QXWV9x/H3h034kWyALAGkIXUvtAwEqkSDUmOJoDJU8FdVBhSEarG1iERjtUVHZaQdBUVtp3XsyM8SRTQRkVEhpikYMWkDhIQYlMpGi2QEgiYBY3DJt3+c55qbzWb37LN7du/Z/bxm7txzn3Puud8nyf3k/LjnOYoIzMyGap+xLsDM6snhMQIkaaxrGI661w/170Md63d4jIynx7qAYap7/VD/PtSufoeHmWVxeJhZFtX5bMuMGTOiu7t7rMugp6eHRqMx1mVkq3v9UP8+tHP9PT09sXnz5j02NCaNRTEjpbu7m9WrV491GWbj2ty5c/s9mOvdFjPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsyyVhoekWZKWS9ogab2kS1vmXSLpx6n9ypb2qyStljS/ytrMbHiqHs+jF1gYEfdJmgbcK2kpcDjweuAFEbFD0mEAko5N7zsFuB64q+L6zCxTpeEREZuATWl6m6QNwEzgIuCTEbEjzXs8vaUD2AkEULvRpM0mklE75iGpG5gDrAKOAf5M0ipJd0k6CSAi1gNTgBXAF/aynsmjUrCZ7aH1+zcqwxBK6gQWAwsiYqukScB04GTgJOAWSUdF4ZIB1iPgV817XHR1dY1C9WYTW09PD5KeSS9D0rSIiMq3PFJSLQYWRcSS1PwosCSFxX9T7KrMGGxdafnOiJgaEVPbdcBYs/Gk0WjQ/M6l719A9WdbBFwDbIiIq1tm3QqclpY5BtgXeLLKWsxsZFW92zIPOB9YJ2lNarsMuBa4VtKDwLPABc00M7N6qPpsywr2ftbkvCo/28yq5V+YmlkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4mFmWSsND0ixJyyVtkLRe0qV95n9AUkiakV7vI+lGSfdIOr7K2sxseKre8ugFFkbEccDJwMWSZkMRLMCrgZ+3LH86sAp4I7Cw4trMbBgqDY+I2BQR96XpbcAGYGaa/Vngg0C0vKUD2JkeqrI2MxueUTvmIakbmAOskvQ64BcR8UCfxe4A5gO3AVfvZT0dFZZpZgOQNLk5PWmUPrATWAwsoNiV+TDFLspuIqIXOGeA9QjYkp7p6uqqpF4z26WnpwdJz6SXIWlaRETlWx4pqRYDiyJiCXA00AAekLQROBK4T9LzBltXFDojYmpETG00GlWWbmZAo9Gg+Z1L37+Airc80hbCNcCGiLgaICLWAYe1LLMRmBsRT1ZZi5mNrKq3POYB5wOnSVqTHq+p+DPNbBRUuuURESsY5KxJRHRXWYOZVcO/MDWzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLKXG85A0D/g48Pz0HlGMCnhUdaWZWTsrOxjQNcD7gHuB56orx8zqomx4bImI71RaiZnVStnwWC7pKmAJsKPZ2Lyhk5lNPGXD46XpeW5LWwCnjWw5ZlYXpcIjIk6tuhAzq5dSp2olHSTpakmr0+Mzkg6qujgza19lf+dxLbANODs9tgLXVVWUmbW/ssc8jo6IN7W8vlzSmioKMrN6KLvlsV3Sy5sv0o/GtldTkpnVQdktj3cDN6TjHAKeAi6sqigza39lz7asAV4o6cD0emulVZlZ2xswPCSdFxE3SXp/n3YAIt353swmnsG2PKam52n9zIsRrsXMamTA8IiIL6bJ70XED1rnpYOmZjZBlT3b8i8l28xsghjsmMefAi8DDu1z3ONAoKPKwsysvQ12zGNfoDMt13rcYyvw5qqKMrP2N9gxj7uAuyRdHxE/G6WazKwGyh7z+I2kqyR9W9J/Nh+DvUnSLEnLJW2QtF7Span9KkkPSVor6RuSDm55z1Xp4rv5mX0ys1FQNjwWAQ8BDeByYCPwPyXe1wssjIjjgJOBiyXNBpYCJ0TEC4CfAP8AIOnY9L5TgItL1mZmY6BseBwSEdcAv4uIuyLiHRRhMKCI2NQcbSwitgEbgJkRcWdE9KbFVgJHpukOYCfFb0g0hH6Y2SgrGx6/S8+bJJ0paQ67vvClSOoG5gCr+sx6B/AdgIhYD0wBVgBf2Mt6fJanhiTt9rB6av3+lb0w7op0UdxCit93HEgxmnrZD+wEFgMLWq+LkfRhil2bRc22iLhkgPUI2JKe6erqKluCjaH+wkISEf6Rch309PQg6Zn0MiRNi4goe2Hc7WlyCzCkIQklTaYIjkURsaSl/QLgLOCVUfJfUVqus/l67ty5/tdnVrFGo8HmzZun9m0ve9OnQ4GLgO7W96RjHwO9TxT3fNnQehGdpDOADwHzI+I3ZWows/ZSdrflm8D3ge8xtJs+zQPOB9a1jDx2GfDPwH7A0rRJuzIi/mYI6zWzMVY2PKZExIeGuvKIWEH/Z02+PdR1WX1FxB7HPXy8o/7Knm25XdJrKq3ExrWI2O1h9Vc2PC6lCJDtkrZK2ibJo4mZTWBlz7b0NxiQmU1gg12Sf2xEPCTpRf3N971qzSauwbY8FlKcov1MP/N8r1qzCWywS/IvSs++V62Z7Waw3Za/GGh+6y9GzWxiGWy35bUDzAvA4WE2QQ222/KXo1WImdVLqd95SPqnPqN9TZd0RXVlmVm7K/sjsT+PiF83X0TErwD/4tRsAisbHh2S9mu+kHQAxYVtZlZC38GQxsOASGUvjLsJWCbpOooDpe8AbqisKrMJoO4DIpX9efqVktYCr6K4SvYTEXFHpZWZjRPjYSujP2W3PKAYvLg3Ir4naUoaimxbVYWZjRf9DUkwHpQ923IR8HWgeePrmcCtVRVlZu2v7AHTiylGBdsKEBEPA4dVVZTZeNPfsY06H++A8rstOyLi2eaml6RJFAdOzaykuodFX2W3PO6SdBlwgKRXA18DvlVdWWbW7sqGx98DTwDrgL+mGIP0I1UVZWbtr+yp2p2SbgVujYgnKq7JzGpgwC0PFT4u6UmKG13/WNITkj46OuWZWbsabLdlAcVZlpMi4pCI6AJeCsyTVPp2k2Y2/gwWHm8Hzo2InmZDRDwCnJfmmdkENVh4TI6IJ/s2puMek6spyczqYLDweDZznpmNc4OdbXnhXm7uJGD/Cuoxs5oYbBjCjtEqxMzqpeyPxMzMduPwMLMsDg8zy+LwMLMsDg8zy+LwMLMsDg8zy1JpeEiaJWm5pA2S1ku6NLV3SVoq6eH0PD217yPpRkn3SDq+ytrMbHiq3vLoBRZGxHHAycDFkmZTDC60LCL+GFiWXgOcDqwC3ggsrLg2MxuGSsMjIjZFxH1pehvF7RtmAq9n102jbgDekKY7gJ3pMf7GqjcbR0btmIekbmAOxZbF4RGxCYqAYddI7HcA84HbgKv3sh7/ZN5sjLR+/4Zy06fhfGAnsBhYEBFb93YDnIjoBc4ZYD0CtqRnurq6KqjWzFr19PQg6Zn0MtIN36LyLQ9JkymCY1FELEnNv5R0RJp/BPB4mXVFoTMipkbE1EajUU3RZvZ7jUaD5ncuff8Cqj/bIuAaYENEtO6G3AZckKYvAL5ZZR1mNvKq3m2ZB5wPrJO0JrVdBnwSuEXSO4GfA2+puA4zG2GVhkdErGDvZ01eWeVnm1m1/AtTM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLA4PM8vi8DCzLJWGh6RrJT0u6cGWthMlrZS0RtJqSS9J7ftIulHSPZKOr7IuMxu+qrc8rgfO6NN2JXB5RJwIfDS9BjgdWAW8EVhYcV1mNkyVhkdE3A081bcZODBNHwQ8lqY7gJ3poSrrMrPhmzQGn7kAuEPSpynC62Wp/Q7gJuDtwLvGoC4zG4KxCI93A++LiMWSzgauAV4VEb3AOYO9WVJHRDxXdZFmtqfW758iouoP6wZuj4gT0ustwMEREZIEbImIAwdYReu6BDzdfN3V1TVl8+bNI1+0mf3eIYccwlNPPfWb9DKAaRERY3Gq9jFgfpo+DXi47BujMLX5aDQalRRoZrs0Gg1avnedkbY4Kt1tkfQV4BXADEmPAh8DLgI+L2kS8Ft8fMOslioNj4g4dy+zXlzl55pZ9fwLUzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPL4vAwsywODzPLUvmtF6o0Y8aM6O7uHrXP6+npYTyN2O7+tL926FNPT09s3rx5jw2NWofHaJP0TERMHes6Ror70/7auU/ebTGzLA4PG098g/RR5N2WIZCkGEd/YO5P+2vnPjk8zCyLd1sSSRdKOmus6zCri0pvN9muJG0EtgHPAb0RMTfNOlvSGcAvI+ITAy2blvs80AF8KSI+mdovBE4FtgObgMnACcDZEfHsaPSvlaQOYDXwi4g4a2+1t1vdfUnaH7gb2I/i3+3XI+JjaV4d+zMLuBF4HrAT+PeI+HyaV4/+RMSEewAbgRl92i4E3pamvzrIsh3AT4GjgH2BB4DZLet5a5pelp4vA+aMUV/fD3wZuH2g2tut7n76IaAzTU8GVgEn17g/RwAvStPTgJ+kumvTH++27G5Leh7sQNBLgP+NiEeiSP2bgde3zN+anp9Iz89S/I85qiQdCZwJfKmleaDa26Lu/kTh6fRycnoE9e3Ppoi4L01vAzYAM6lRfyZqeARwp6R7Jb0rY9mZwP+1LPNoams3nwM+SLFZ3FSX2vcgqUPSGuBxYGlErKLG/WmS1A3Modiaqk1/JuQxD2BeRDwm6TBgqaSHIuL65syIOGegZen/9wSR3rvHeiLi0xX0YUDp4O/jEXGvpFe0zupn8WiXugcSEc8BJ0o6GPiGpBOocX8AJHUCi4EFEbFVUm36MyG3PCLisfT8OPANik3FoSz7KDCrZbEjgceqqjfTPOB16YDvzcBpkm6iHrUPKCJ+DfwXcAY17o+kyRTBsSgilqTm+vRnrA8cjfYDmApMa5m+BzhjKMtSbLE9AjTYdVDr+LHu2wB9fgW7DpjWqvaWPhwKHJymDwC+D5xV4/6I4mzL5/q016Y/E3G35XCKTV4o/qK+HBHfHeqykt4D3EFxdPzaiFhfdeEjISJ6a1r7EcAN6dTzPsAtEXE71PbvYh5wPrAuHccBuCwivl2X/vgXpmaWZUIe8zCz4XN4mFkWh4eZZXF4mFkWh4eZZXF4mFkWh4eZZXF4jEOSnpO0RtKDkr6VrgVB0h9I+nqJ9z+9l/Y3SJo9yHsfkPSVvMpHRtl+2vA4PMan7RFxYkScADwFXAzFdToR8eZhrPcNFGNL9EvScRT/pk6RNGa3CxiBfloJDo/x74ekS7oldUt6ME1PkXSLpLWSvipplaTmiGpI+se0FbFS0uGSXga8DrgqbdUc3c9nvRX4D+DOtGxzXe+V9KP0WTentk5J10lal9rflNpPl/RDSfdJ+lq66hRJGyVdntrXSTo2tc9P9ayRdL+kaX36uX/L59wv6dTUfqGkJZK+K+lhSVeO8J/7+DfWF9f4MfIP4On03AF8jXThH9ANPJimPwB8MU2fAPQCc9PrAF6bpq8EPpKmrwfePMDn/gR4PnA6cFtL+2PAfmm6eXHbp2i5KAyYDsygGGpwamr7EPDRNL0RuCRN/y3F8HwA36IYNgGgk+IapNZ+LgSuS9PHAj8H9qcYmesR4KD0+mfArLH+u6vTw1se49MB6WKrzUAXsLSfZV5Ocak+EfEgsLZl3rPA7Wn6Xoov44AknQQ8ERE/A5YBL5I0Pc1eCyySdB5FSAG8CvjX5vsj4lcUwwrOBn6Q6r+AIoyampett9b0A+BqSe+lCKZedvdyiq0hIuIhipA4Js1bFhFbIuK3wI/6fJYNwuExPm2PiBMpvgz7ko559DHQDZJ+F+m/aoqBn8tcfX0ucGwaP+SnwIHAm9K8MymC4sXAvZImpc/ve1WmKEYIOzE9ZkfEO1vm7+hbUxQDT/8VxWX6K5u7MyX7uaNlumw/LXF4jGMRsQV4L/CBNPBMqxXA2QDpDMqflFjlNorBencjaR/gLcALIqI7Iropxt08N82bFRHLKYZEPJhi9+JO4D0t65gOrATmSfqj1DZF0jEMQNLREbEuIj5FMUp83/C4G3hbWvYY4A+BH5foqw3C4THORcT9FAPKnNNn1r8Bh0paS3FsYS27BoDem5uBv0sHHlsPmJ5CcWuHX7S03U2xCzITuEnSOuB+4LNRjAR2BTA9nU5+ADg1Ip6gOBbxlVTXSvYMg74WtKxjO/CdfvrZkT7/q8CFEbGj70ps6DyexwSVBtWZHBG/TUGwDDgm2uCeJlYP3sebuKYAy9PujIB3OzhsKLzlYWZZfMzDzLI4PMwsi8PDzLI4PMwsi8PDzLL8P0hjCDuCzcLHAAAAAElFTkSuQmCC\n",
"text/plain": [
"