{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8cbe167",
   "metadata": {},
   "source": [
    "# Recept Richtingsreconstructie\n",
    "\n",
    "Deze notebooks werken alleen met Python 3.\n",
    "\n",
    "Richtingsrecontructie op basis van aankomsttijden van deeltjes is mogelijk als\n",
    "er 3 tijden zijn gemeten. In het geval van HiSPARC meetstations, is dat mogelijk\n",
    "bij stations met 4 detectoren.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf2c268b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tables\n",
    "from sapphire import download_data\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13573cb7",
   "metadata": {},
   "source": [
    "In het recept *download data* is reeds data gedownload in het bestand 'data.h5':\n",
    "\n",
    "(Als dit bestand nog niet bestaat wordt het nu gemaakt. We downloaden de\n",
    "benodigde data dan verderop alsnog)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aea4760a",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILENAME = 'data5.h5'\n",
    "data = tables.open_file(FILENAME, 'a')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4be44268",
   "metadata": {},
   "source": [
    "In het recept *download data* zijn reeds de events van station 501 van de eerste\n",
    "twee dagen van januari 2016 gedownload, in de tabel 's501'.\n",
    "\n",
    "Als de data nog niet is gedownload, dan downloaden we de data hier alsnog:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "998af45b",
   "metadata": {},
   "outputs": [],
   "source": [
    "if '/s501' not in data:\n",
    "    download_data(data, '/s501', 501, start=datetime(2016, 1, 1), end=datetime(2016,1,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a223d2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a43d6cbd",
   "metadata": {},
   "source": [
    "We gebruiken de SAPPHiRE class `ReconstructESDEvents`:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "633e7d9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sapphire import ReconstructESDEvents\n",
    "rec = ReconstructESDEvents(data, '/s501', 501)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aef561ba",
   "metadata": {},
   "source": [
    "De function `reconstruct_directions()` reconstrueert de richtingen (hoeken) van\n",
    "elk event uit de aankomsttijden per detector en slaat de uitkomst op in:\n",
    "- `rec.theta`: een array met de zenith hoek per event\n",
    "- `rec.phi`: een array met de azimuth hoek per event\n",
    "\n",
    "De hoeken kunnen ook meteen opgeslagen worden in het HDF5 bestand. Dit wordt in\n",
    "de volgende paragraaf toegelicht.\n",
    "\n",
    "Eerst reconstrueren we de aankomstrichting van de events:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf68b0ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec.reconstruct_directions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5d3b7f",
   "metadata": {},
   "source": [
    "Hieronder staan de eerste twintig zenithoeken. 'nan' (NaN) betekent:\n",
    "Not-a-number: De reconstructie was niet mogelijk. Ofwel er was onvoldoende\n",
    "informatie, bijvoorbeeld slechts twee aankomsttijden. Ofwel de oplossing was\n",
    "niet fysisch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc47adcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec.theta[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3728ecb",
   "metadata": {},
   "source": [
    "Met behulp van de functie `numpy.isnan()` kunnen we de NaNs verwijderen:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72008379",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import isnan\n",
    "zenith = [a for a in rec.theta if not isnan(a)]\n",
    "azimuth = [a for a in rec.phi if not isnan(a)]\n",
    "print(\"Er zijn %d events succesvol gereconstrueerd.\" % len(zenith))            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cca3a0",
   "metadata": {},
   "source": [
    "Nu bevat de array `zenith` slechts hoeken (in radialen):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a66e1795",
   "metadata": {},
   "outputs": [],
   "source": [
    "zenith[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec8d70d5",
   "metadata": {},
   "source": [
    "## Plotten\n",
    "\n",
    "We maken een polar-plot van de hoeken theta en phi (zenit, azimut) van de\n",
    "(succesvol) gereconstrueerde events, en een histrogram van de zenithoeken.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e869d608",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65f46a8d",
   "metadata": {},
   "source": [
    "We rekenen de hoeken om naar graden:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc4e112d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import degrees\n",
    "zenith = degrees(zenith)\n",
    "azimuth = degrees(azimuth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a55f5a5",
   "metadata": {},
   "source": [
    "Polar plot van zenit en azimuth:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7687240d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(polar=True)\n",
    "ax.scatter(zenith, azimuth)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0030d2e3",
   "metadata": {},
   "source": [
    "Histogram van de zenithoeken:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32f382ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import arange\n",
    "plt.hist(zenith, bins=arange(0,90., 5.), histtype='step')\n",
    "plt.xlabel('zenith angle (degrees)')\n",
    "plt.ylabel('counts')\n",
    "plt.title('Zenith histrogram station 501.')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6476c9a8",
   "metadata": {},
   "source": [
    "# Opslaan van reconstructies in het HDF5 bestand\n",
    "\n",
    "De SAPPHiRE reconstructie class `ReconstructESDEvents` kan de reconstructies ook\n",
    "direct opslaan in het HDF5 bestand waar de events in zijn opgeslagen:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e426701",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = ReconstructESDEvents(data, '/s501', 501)\n",
    "rec.reconstruct_and_store()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c80518b",
   "metadata": {},
   "source": [
    "Er waren nu twee progressbars te zien: De eerste voor richting reconstructie en\n",
    "de tweede voor core recontructie. We richten ons hier op de richtingen en laten\n",
    "de core reconstructie buiten beschouwing.\n",
    "\n",
    "naast `rec.theta` en `rec.phi` is er nu ook een nieuwe groep\n",
    "`/s501/reconstructions` in het hdf5 bestand:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9b1b484",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec_tabel = data.root.s501.reconstructions\n",
    "print(rec_tabel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d3208d0",
   "metadata": {},
   "source": [
    "De gereconstrueerde hoeken zijn nu ook opgeslagen op disk, in het HDF5 bestand.\n",
    "\n",
    "De kolommen 'zenith' en 'azimuth' bevatten de gereconstrueerde hoeken (of NaN):\n",
    "\n",
    "## Opgave:\n",
    "\n",
    "Maak een histogram van de zenithoeken en een polarplot van de azimuth en\n",
    "zenithhoeken, zoals in het voorbeeld hierboven. Gebruik de `zenith` en `azimuth`\n",
    "kolommen uit de groep `/s501/reconstructions` uit het HDF5 bestand.\n",
    "\n",
    "*Hints*:\n",
    "- gebruik: `.col('zenith')` om een kolom in te laden.\n",
    "- gebruik  `.compress(~isnan(...))` om te NaNs verwijderen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50fe9a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "zenith = rec_tabel.col('zenith')\n",
    "azimuth = rec_tabel.col('azimuth')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9874043d",
   "metadata": {},
   "outputs": [],
   "source": [
    "zenith = zenith.compress(~isnan(zenith))\n",
    "azimuth = azimuth.compress(~isnan(azimuth))\n",
    "print(\"Er zijn %d events succesvol gereconstrueerd.\" % len(zenith))           "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c312f64f",
   "metadata": {},
   "outputs": [],
   "source": [
    "zenith = degrees(zenith)\n",
    "azimuth = degrees(azimuth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a579c568",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plt.subplot(polar=True)\n",
    "ax.scatter(zenith, azimuth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1eac3769",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(zenith, bins=arange(0,90., 5.), histtype='step')\n",
    "plt.xlabel('zenith angle (degrees)')\n",
    "plt.ylabel('counts')\n",
    "plt.title('Zenith histrogram station 501.')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}