docs/notebooks/comparison_with_EventDisplay.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison with EventDisplay"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Purpose of this notebook:**\n",
    "\n",
    "Compare IRF and Sensitivity as computed by pyirf and EventDisplay on the same DL2 results\n",
    "\n",
    "**Notes:**\n",
    "\n",
    "The following results correspond to:\n",
    "\n",
    "- Paranal site\n",
    "- Zd 20 deg, Az 180 deg\n",
    "- 50 h observation time\n",
    "\n",
    "**Resources:**\n",
    "\n",
    "_EventDisplay_ DL2 data, https://forge.in2p3.fr/projects/cta_analysis-and-simulations/wiki/Eventdisplay_Prod3b_DL2_Lists\n",
    "\n",
    "\n",
    "Download and unpack the data using \n",
    "\n",
    "```bash\n",
    "$ curl -fL -o data.zip https://nextcloud.e5.physik.tu-dortmund.de/index.php/s/Cstsf8MWZjnz92L/download\n",
    "$ unzip data.zip\n",
    "$ mv eventdisplay_dl2 data\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents\n",
    "\n",
    "* [Optimized cuts](#Optimized-cuts)\n",
    "    - [Direction cut](#Direction-cut)\n",
    "* [Differential sensitivity from cuts optimization](#Differential-sensitivity-from-cuts-optimization)\n",
    "* [IRFs](#IRFs)\n",
    "    - [Effective area](#Effective-area)\n",
    "    - [Point Spread Function](#Point-Spread-Function)\n",
    "        + [Angular resolution](#Angular-resolution)\n",
    "    - [Energy dispersion](#Energy-dispersion)\n",
    "        + [Energy resolution](#Energy-resolution)\n",
    "    - [Background rate](#Background-rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "## Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import uproot\n",
    "from astropy.io import fits\n",
    "import astropy.units as u\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.table import QTable\n",
    "from matplotlib.ticker import ScalarFormatter\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (9, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "## Input data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "### _EventDisplay_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "The input data provided by _EventDisplay_ is stored in _ROOT_ format, so _uproot_ is used to transform it into _numpy_ objects. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "# Path of EventDisplay IRF data in the user's local setup\n",
    "# Please, empty the indir_EventDisplay variable before pushing to the repo\n",
    "indir = \"../../data/\"\n",
    "irf_file_event_display = \"DESY.d20180113.V3.ID0NIM2LST4MST4SST4SCMST4.prod3b-paranal20degs05b-NN.S.3HB9-FD.180000s.root\"\n",
    "\n",
    "irf_eventdisplay = uproot.open(os.path.join(indir, irf_file_event_display))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "## _pyirf_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following is the current IRF + sensititivy output FITS format provided by this software.\n",
    "\n",
    "Run `python examples/calculate_eventdisplay_irfs.py` after downloading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyirf_file = '../../pyirf_eventdisplay.fits.gz'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimized cuts\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Direction cut\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from astropy.table import QTable\n",
    "\n",
    "\n",
    "rad_max = QTable.read(pyirf_file, hdu='RAD_MAX')[0]\n",
    "\n",
    "\n",
    "theta_cut_ed = irf_eventdisplay['ThetaCut;1']\n",
    "plt.errorbar(\n",
    "    10**theta_cut_ed.edges[:-1],\n",
    "    theta_cut_ed.values,\n",
    "    xerr=np.diff(10**theta_cut_ed.edges),\n",
    "    ls='',\n",
    "    label='EventDisplay',\n",
    ")\n",
    "\n",
    "plt.errorbar(\n",
    "    0.5 * (rad_max['ENERG_LO'] + rad_max['ENERG_HI'])[1:-1].to_value(u.TeV),\n",
    "    rad_max['RAD_MAX'].T[1:-1, 0].to_value(u.deg),\n",
    "    xerr=0.5 * (rad_max['ENERG_HI'] - rad_max['ENERG_LO'])[1:-1].to_value(u.TeV),\n",
    "    ls='',\n",
    "    label='pyirf',\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.ylabel('θ-cut / deg²')\n",
    "plt.xlabel(r'$E_\\mathrm{reco} / \\mathrm{TeV}$')\n",
    "plt.xscale('log')\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.table import QTable\n",
    "\n",
    "\n",
    "gh_cut = QTable.read(pyirf_file, hdu='GH_CUTS')[1:-1]\n",
    "\n",
    "\n",
    "plt.errorbar(\n",
    "    0.5 * (gh_cut['low'] + gh_cut['high']).to_value(u.TeV),\n",
    "    gh_cut['cut'],\n",
    "    xerr=0.5 * (gh_cut['high'] - gh_cut['low']).to_value(u.TeV),\n",
    "    ls='',\n",
    "    label='pyirf',\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.ylabel('G/H-cut')\n",
    "plt.xlabel(r'$E_\\mathrm{reco} / \\mathrm{TeV}$')\n",
    "plt.xscale('log')\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential sensitivity from cuts optimization\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# [1:-1] removes under/overflow bins\n",
    "sensitivity = QTable.read(pyirf_file, hdu='SENSITIVITY')[1:-1]\n",
    "\n",
    "# make it print nice\n",
    "sensitivity['reco_energy_low'].info.format = '.3g'\n",
    "sensitivity['reco_energy_high'].info.format = '.3g'\n",
    "sensitivity['reco_energy_center'].info.format = '.3g'\n",
    "sensitivity['n_signal'].info.format = '.1f'\n",
    "sensitivity['n_signal_weighted'].info.format = '.1f'\n",
    "sensitivity['n_background_weighted'].info.format = '.1f'\n",
    "sensitivity['n_background'].info.format = '.1f'\n",
    "sensitivity['relative_sensitivity'].info.format = '.2g'\n",
    "sensitivity['flux_sensitivity'].info.format = '.3g'\n",
    "\n",
    "\n",
    "sensitivity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,8))\n",
    "\n",
    "# Get data from event display file\n",
    "h = irf_eventdisplay[\"DiffSens\"]\n",
    "bins = 10**h.edges\n",
    "x = 0.5 * (bins[:-1] + bins[1:])\n",
    "width = np.diff(bins)\n",
    "y = h.values\n",
    "yerr = np.sqrt(h.variances)\n",
    "\n",
    "fig, (ax_sens, ax_ratio) = plt.subplots(\n",
    "    2, 1,\n",
    "    gridspec_kw={'height_ratios': [4, 1]},\n",
    "    sharex=True,\n",
    ")\n",
    "\n",
    "ax_sens.errorbar(\n",
    "    x,\n",
    "    y, \n",
    "    xerr=width/2,\n",
    "    yerr=yerr,\n",
    "    label=\"EventDisplay\",\n",
    "    ls=''\n",
    ")\n",
    "\n",
    "unit = u.Unit('erg cm-2 s-1')\n",
    "\n",
    "\n",
    "e = sensitivity['reco_energy_center']\n",
    "s = (e**2 * sensitivity['flux_sensitivity'])\n",
    "w = (sensitivity['reco_energy_high'] - sensitivity['reco_energy_low'])\n",
    "\n",
    "ax_sens.errorbar(\n",
    "    e.to_value(u.TeV),\n",
    "    s.to_value(unit),\n",
    "    xerr=w.to_value(u.TeV) / 2,\n",
    "    ls='',\n",
    "    label='pyirf'\n",
    ")\n",
    "\n",
    "ax_ratio.errorbar(\n",
    "    e.to_value(u.TeV), s.to_value(unit) / y,\n",
    "    xerr=w.to_value(u.TeV)/2,\n",
    "    ls=''\n",
    ")\n",
    "ax_ratio.set_yscale('log')\n",
    "ax_ratio.set_xlabel(\"Reconstructed energy / TeV\")\n",
    "ax_ratio.set_ylabel('pyirf / eventdisplay')\n",
    "ax_ratio.grid()\n",
    "ax_ratio.yaxis.set_major_formatter(ScalarFormatter())\n",
    "\n",
    "ax_ratio.set_ylim(0.2, 5.0)\n",
    "ax_ratio.set_yticks([0.2, 0.5, 1, 2, 5])\n",
    "ax_ratio.set_yticks([], minor=True)\n",
    "\n",
    "# Style settings\n",
    "ax_sens.set_title('Minimal Flux Needed for 5σ Detection in 50 hours')\n",
    "ax_sens.set_xscale(\"log\")\n",
    "ax_sens.set_yscale(\"log\")\n",
    "ax_sens.set_ylabel(rf\"$(E^2 \\cdot \\mathrm{{Flux Sensitivity}}) /$ ({unit.to_string('latex')})\")\n",
    "\n",
    "ax_sens.grid(which=\"both\")\n",
    "ax_sens.legend()\n",
    "fig.tight_layout(h_pad=0)\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IRFs\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Effective area\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Data from EventDisplay\n",
    "h = irf_eventdisplay[\"EffectiveAreaEtrue\"]\n",
    "\n",
    "x = 0.5 * (10**h.edges[:-1] + 10**h.edges[1:])\n",
    "xerr = 0.5 * np.diff(10**h.edges)\n",
    "y = h.values\n",
    "yerr = np.sqrt(h.variances)\n",
    "\n",
    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls='', label=\"EventDisplay\")\n",
    "\n",
    "for name in ('', '_NO_CUTS', '_ONLY_GH', '_ONLY_THETA'):\n",
    "\n",
    "    area = QTable.read(pyirf_file, hdu='EFFECTIVE_AREA' + name)[0]\n",
    "\n",
    "    \n",
    "    plt.errorbar(\n",
    "        0.5 * (area['ENERG_LO'] + area['ENERG_HI']).to_value(u.TeV)[1:-1],\n",
    "        area['EFFAREA'].to_value(u.m**2).T[1:-1, 0],\n",
    "        xerr=0.5 * (area['ENERG_LO'] - area['ENERG_HI']).to_value(u.TeV)[1:-1],\n",
    "        ls='',\n",
    "        label='pyirf ' + name,\n",
    "    )\n",
    "\n",
    "# Style settings\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"True energy / TeV\")\n",
    "plt.ylabel(\"Effective collection area / m²\")\n",
    "plt.grid(which=\"both\")\n",
    "plt.legend()\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Point Spread Function\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psf_table = QTable.read(pyirf_file, hdu='PSF')[0]\n",
    "# select the only fov offset bin\n",
    "psf = psf_table['RPSF'].T[:, 0, :].to_value(1 / u.sr)\n",
    "\n",
    "offset_bins = np.append(psf_table['RAD_LO'], psf_table['RAD_HI'][-1])\n",
    "phi_bins = np.linspace(0, 2 * np.pi, 100)\n",
    "\n",
    "\n",
    "\n",
    "# Let's make a nice 2d representation of the radially symmetric PSF\n",
    "r, phi = np.meshgrid(offset_bins.to_value(u.deg), phi_bins)\n",
    "\n",
    "# look at a single energy bin\n",
    "# repeat values for each phi bin\n",
    "center = 0.5 * (psf_table['ENERG_LO'] + psf_table['ENERG_HI'])\n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=(15, 5))\n",
    "axs = [fig.add_subplot(1, 3, i, projection='polar') for i in range(1, 4)]\n",
    "\n",
    "\n",
    "for bin_id, ax in zip([10, 20, 30], axs):\n",
    "    image = np.tile(psf[bin_id], (len(phi_bins) - 1, 1))\n",
    "    \n",
    "    ax.set_title(f'PSF @ {center[bin_id]:.2f} TeV')\n",
    "    ax.pcolormesh(phi, r, image)\n",
    "    ax.set_ylim(0, 0.25)\n",
    "    ax.set_aspect(1)\n",
    "    \n",
    "fig.tight_layout()\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Profile\n",
    "center = 0.5 * (offset_bins[1:] + offset_bins[:-1])\n",
    "xerr = 0.5 * (offset_bins[1:] - offset_bins[:-1])\n",
    "\n",
    "for bin_id in [10, 20, 30]:\n",
    "    plt.errorbar(\n",
    "        center.to_value(u.deg),\n",
    "        psf[bin_id],\n",
    "        xerr=xerr.to_value(u.deg),\n",
    "        ls='',\n",
    "        label=f'Energy Bin {bin_id}'\n",
    "    )\n",
    "    \n",
    "#plt.yscale('log')\n",
    "plt.legend()\n",
    "plt.xlim(0, 0.25)\n",
    "plt.ylabel('PSF PDF / sr⁻¹')\n",
    "plt.xlabel('Distance from True Source / deg')\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Angular resolution\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data from EventDisplay\n",
    "h = irf_eventdisplay[\"AngRes\"]\n",
    "x = 0.5 * (10**h.edges[:-1] + 10**h.edges[1:])\n",
    "xerr = 0.5 * np.diff(10**h.edges)\n",
    "y = h.values\n",
    "yerr = np.sqrt(h.variances)\n",
    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls='', label=\"EventDisplay\")\n",
    "\n",
    "# pyirf\n",
    "\n",
    "ang_res = QTable.read(pyirf_file, hdu='ANGULAR_RESOLUTION')[1:-1]\n",
    "\n",
    "plt.errorbar(\n",
    "    0.5 * (ang_res['true_energy_low'] + ang_res['true_energy_high']).to_value(u.TeV),\n",
    "    ang_res['angular_resolution'].to_value(u.deg),\n",
    "    xerr=0.5 * (ang_res['true_energy_high'] - ang_res['true_energy_low']).to_value(u.TeV),\n",
    "    ls='',\n",
    "    label='pyirf'\n",
    ")\n",
    "\n",
    "\n",
    "# Style settings\n",
    "plt.xlim(1.e-2, 2.e2)\n",
    "plt.ylim(2.e-2, 1)\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"True energy / TeV\")\n",
    "plt.ylabel(\"Angular Resolution / deg\")\n",
    "plt.grid(which=\"both\")\n",
    "plt.legend(loc=\"best\")\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Energy dispersion\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "edisp = QTable.read(pyirf_file, hdu='ENERGY_DISPERSION')[0]\n",
    "\n",
    "e_bins = edisp['ENERG_LO'][1:]\n",
    "migra_bins = edisp['MIGRA_LO'][1:]\n",
    "\n",
    "plt.title('pyirf')\n",
    "plt.pcolormesh(e_bins.to_value(u.TeV), migra_bins, edisp['MATRIX'].T[1:-1, 1:-1, 0].T, cmap='inferno')\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.colorbar(label='PDF Value')\n",
    "\n",
    "plt.xlabel(r'$E_\\mathrm{True} / \\mathrm{TeV}$')\n",
    "plt.ylabel(r'$E_\\mathrm{Reco} / E_\\mathrm{True}$')\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Energy resolution\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data from EventDisplay\n",
    "h = irf_eventdisplay[\"ERes\"]\n",
    "x = 0.5 * (10**h.edges[:-1] + 10**h.edges[1:])\n",
    "xerr = np.diff(10**h.edges) / 2\n",
    "y = h.values\n",
    "yerr = np.sqrt(h.variances)\n",
    "\n",
    "# Data from pyirf\n",
    "bias_resolution = QTable.read(pyirf_file, hdu='ENERGY_BIAS_RESOLUTION')[1:-1]\n",
    "\n",
    "# Plot function\n",
    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls='', label=\"EventDisplay\")\n",
    "plt.errorbar(\n",
    "    0.5 * (bias_resolution['true_energy_low'] + bias_resolution['true_energy_high']).to_value(u.TeV),\n",
    "    bias_resolution['resolution'],\n",
    "    xerr=0.5 * (bias_resolution['true_energy_high'] - bias_resolution['true_energy_low']).to_value(u.TeV),\n",
    "    ls='',\n",
    "    label='pyirf'\n",
    ")\n",
    "plt.xscale('log')\n",
    "\n",
    "# Style settings\n",
    "plt.xlabel(r\"$E_\\mathrm{True} / \\mathrm{TeV}$\")\n",
    "plt.ylabel(\"Energy resolution\")\n",
    "plt.grid(which=\"both\")\n",
    "plt.legend(loc=\"best\")\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Background rate\n",
    "[back to top](#Table-of-contents)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyirf.utils import cone_solid_angle\n",
    "\n",
    "# Data from EventDisplay\n",
    "h = irf_eventdisplay[\"BGRate\"]\n",
    "x = 0.5 * (10**h.edges[:-1] + 10**h.edges[1:])\n",
    "width = np.diff(10**h.edges)\n",
    "xerr = width / 2\n",
    "y = h.values\n",
    "yerr = np.sqrt(h.variances)\n",
    "\n",
    "# pyirf data\n",
    "\n",
    "bg_rate = QTable.read(pyirf_file, hdu='BACKGROUND')[0]\n",
    "\n",
    "reco_bins = np.append(bg_rate['ENERG_LO'], bg_rate['ENERG_HI'][-1])\n",
    "\n",
    "# first fov bin, [0, 1] deg\n",
    "fov_bin = 0\n",
    "rate_bin = bg_rate['BKG'].T[:, fov_bin]\n",
    "\n",
    "\n",
    "# interpolate theta cut for given e reco bin\n",
    "e_center_bg = 0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI'])\n",
    "e_center_theta = 0.5 * (rad_max['ENERG_LO'] + rad_max['ENERG_HI'])\n",
    "theta_cut = np.interp(e_center_bg, e_center_theta, rad_max['RAD_MAX'].T[:, 0])\n",
    "\n",
    "# undo normalization\n",
    "rate_bin *= cone_solid_angle(theta_cut)\n",
    "rate_bin *= np.diff(reco_bins)\n",
    "\n",
    "\n",
    "# Plot function\n",
    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls='', label=\"EventDisplay\")\n",
    "\n",
    "plt.errorbar(\n",
    "    0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI']).to_value(u.TeV)[1:-1],\n",
    "    rate_bin.to_value(1 / u.s)[1:-1],\n",
    "    xerr=np.diff(reco_bins).to_value(u.TeV)[1:-1] / 2,\n",
    "    ls='',\n",
    "    label='pyirf',\n",
    ")\n",
    "\n",
    "# Style settings\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(r\"$E_\\mathrm{Reco} / \\mathrm{TeV}$\")\n",
    "plt.ylabel(\"Background rate / (s⁻¹ TeV⁻¹) \")\n",
    "plt.grid(which=\"both\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.yscale('log')\n",
    "\n",
    "None # to remove clutter by mpl objects"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}