{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c6a20994",
   "metadata": {},
   "source": [
    "If you have not used gammapy before, you'll want to follow this link https://docs.gammapy.org/1.1/getting-started/install.html#installation to install the package first. Even if you will not be working with Gammapy catalogs with this analysis, it could be useful to go ahead and downoal the tutorial materials in this link https://docs.gammapy.org/dev/getting-started/index.html \n",
    "\n",
    "This notebook adapts the gammapy multi instrument joint 3D and 1D analysis to the context of a vegas fits file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb6fc206",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gammapy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "72e4b19e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy\n",
    "from astropy.io import fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e17342e",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nevis/milne/files/cmm2399/.conda/envs/gammapy-1.1/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "from astropy import units as u\n",
    "import matplotlib.pyplot as plt\n",
    "from gammapy.datasets import Datasets, FluxPointsDataset, SpectrumDatasetOnOff\n",
    "from gammapy.estimators import FluxPoints, FluxPointsEstimator\n",
    "from gammapy.maps import MapAxis\n",
    "from gammapy.modeling import Fit\n",
    "from gammapy.modeling.models import Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8074d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "#We'll need to edit some headers in order for gammapy to read of files, this will let us see the extension (1 SpectrumTable) we need to edit\n",
    "fits_image_filename = astropy.io.fits.open('/path/to/your.fits')\n",
    "fits.info('/path/to/your.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b08c245",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Before modifications:\")\n",
    "print()\n",
    "print(\"Extension 1:\")\n",
    "print(repr(fits.getheader('/path/to/your.fits', 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d841be38",
   "metadata": {},
   "source": [
    "You should look at the SED specifications for the headers and make sure to modify the right headers in accordance, but the changes I specify should work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4352e00",
   "metadata": {},
   "outputs": [],
   "source": [
    "with fits.open('/path/to/your.fits', 'update') as f:\n",
    "    for hdu in f:\n",
    "        hdu.header['TTYPE3'] = 'e_min'\n",
    "        hdu.header['TTYPE4'] = 'e_max'\n",
    "        hdu.header['TTYPE5'] = 'dnde'\n",
    "        hdu.header['TTYPE6'] = 'dnde_err'\n",
    "        hdu.header['TTYPE7'] = 'dnde_ul'\n",
    "        \n",
    "print(\"After modifications:\")\n",
    "print()\n",
    "print(\"Extension 1:\")\n",
    "print(repr(fits.getheader('/path/to/your.fits', 1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "041435df",
   "metadata": {},
   "source": [
    "Remember to repeat the previous steps with all the fits files you plan on using. You will probably need to do some more edits if your file looks does not look like the example hawc file in this tutorial https://docs.gammapy.org/1.1/tutorials/analysis-3d/analysis_mwl.html#hawc-1d-dataset-for-flux-point-fitting, which you should have (if you downloaded the tutorial dataset) at ./gammapy-data/1.1/hawc_crab/HAWC19_flux_points.fits\n",
    "\n",
    "you can open it using \n",
    "    fits_image_filename = astropy.io.fits.open()\n",
    "    fits.info()\n",
    "and pring the first extension (a table) with\n",
    "    print(\"Before modifications:\")\n",
    "    print()\n",
    "    print(\"Extension 1:\")\n",
    "    print(repr(fits.getheader('./gammapy-data/1.1/hawc_crab/HAWC19_flux_points.fits', 1)))\n",
    "\n",
    "Notice it has TS and is_UL columns, which the Vegas files do not. The following steps show how you can emulate this structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afcc0af7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "data = astropy.io.fits.open('/path/to/your.fits', mode='update', extname='SpectrumTable')\n",
    "table = Table.read('/path/to/your.fits', format='fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f093c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "table #Do this if you don't know how many rows of data your extension 1 table has"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf8e5b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This will calculate the TS \n",
    "import math\n",
    "sigval2 = []\n",
    "boolean = []\n",
    "for i in range(22):\n",
    "    x = table['Significance'][i] \n",
    "    y = x**2\n",
    "    sigvals.append(y)\n",
    "    if y < 4:\n",
    "        boolean.append(True)\n",
    "    else:\n",
    "        boolean.append(False)\n",
    "    #print(y)\n",
    "print(sigvals)\n",
    "print(boolean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81354039",
   "metadata": {},
   "outputs": [],
   "source": [
    "#make sure 12 and 13 are the right column positions for you. this will save the editted table to a new fits file, but if you wish to simply overwrite your previous table, use the same path as before and add overwrite=True.\n",
    "col12 = astropy.table.Column(name= 'TS', data=sigval2)\n",
    "col13 = astropy.table.Column(name='is_ul', data=boolean)\n",
    "table.add_column(col12)\n",
    "table.add_column(col13)\n",
    "table.write('/path/to/your/new.fits', format='fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31258009",
   "metadata": {},
   "outputs": [],
   "source": [
    "table #now you can see that the new columns were created"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdb8b3ce",
   "metadata": {},
   "source": [
    "Repeat these steps for your other fits file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13ee8fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import Table\n",
    "data = astropy.io.fits.open('/path/to/your/other.fits', mode='update', extname='SpectrumTable')\n",
    "table = Table.read('/path/to/your/other.fits', format='fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "06a9aa20",
   "metadata": {},
   "outputs": [],
   "source": [
    "table #check number of rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a118433",
   "metadata": {},
   "outputs": [],
   "source": [
    "#This will calculate the TS \n",
    "import math\n",
    "sigval2 = []\n",
    "boolean = []\n",
    "for i in range(22):\n",
    "    x = table['Significance'][i] \n",
    "    y = x**2\n",
    "    sigvals.append(y)\n",
    "    if y < 4:\n",
    "        boolean.append(True)\n",
    "    else:\n",
    "        boolean.append(False)\n",
    "    #print(y)\n",
    "print(sigvals)\n",
    "print(boolean)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12718bd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#make sure 12 and 13 are the right column positions for you. this will save the editted table to a new fits file, but if you wish to simply overwrite your previous table, use the same path as before and add overwrite=True.\n",
    "col12 = astropy.table.Column(name= 'TS', data=sigval2)\n",
    "col13 = astropy.table.Column(name='is_ul', data=boolean)\n",
    "table.add_column(col12)\n",
    "table.add_column(col13)\n",
    "table.write('/path/to/your/other/new.fits', format='fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a157e7ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "table #now you can see that the new columns were created"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "857e70e0",
   "metadata": {},
   "source": [
    "Now we can call the new files and actually create the plot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec36ca83",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = \"/path/to/your/new.fits\"\n",
    "flux_points_yourfile = FluxPoints.read(filename, hdu=1, sed_type=\"dnde\")\n",
    "\n",
    "dataset = FluxPointsDataset(data=flux_points_yourfile, name=\"yourfile\")\n",
    "\n",
    "print(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed1940bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "flux_points_yourfile.to_table() #check everything you need to be read is... there"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d32de3b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check things for your second file\n",
    "filename = \"/path/to/your/other/new.fits\"\n",
    "flux_points_yourotherfile = FluxPoints.read(filename, hdu=1, sed_type=\"dnde\")\n",
    "\n",
    "dataset = FluxPointsDataset(data=flux_points_yourotherfile, name=\"yourotherfile\")\n",
    "\n",
    "print(dataset)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e49c4bdc",
   "metadata": {},
   "source": [
    "Before moving on to creating the plot, change all of the dnde and dnde_err cells that are upper limits (True) to Null and all the dnde_ul cells that are not an upper limit (False) to Null. You can use the table.replace() command and manually do this, or create a boolean using is_UL. I did this manually with fv in my terminal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed913a3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now everything should be working well, and you can create your plot!\n",
    "fig, ax = plt.subplots(figsize=(8, 6))\n",
    "\n",
    "energy_bounds = [0.1, 100] * u.TeV\n",
    "sed_type = \"dnde\"\n",
    "\n",
    "flux_points_Martinez.plot(ax=ax, sed_type=sed_type, label=\"yourfile\")\n",
    "flux_points_Drake.plot(ax=ax, sed_type=sed_type, label=\"yourotherfile\")\n",
    "\n",
    "ax.set_xlim(energy_bounds)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gammapy-1.1",
   "language": "python",
   "name": "gammapy-1.1"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}