{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6e01ed58",
   "metadata": {},
   "source": [
    "# Recept: Detector tijd offsets\n",
    "\n",
    "Deze notebooks werken alleen met Python 3.\n",
    "\n",
    "Tijdsverschillen tussen detectoren"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8807c460",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b986ad8e",
   "metadata": {},
   "source": [
    "## Events\n",
    "\n",
    "1 dag data van station 501"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a601ba6",
   "metadata": {},
   "outputs": [],
   "source": [
    "FILENAME = 'events.h5'\n",
    "station = 501"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05afc903",
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "start = datetime(2015, 5, 1)\n",
    "end = datetime(2015, 5, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e21d4f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tables\n",
    "data = tables.open_file(FILENAME, 'a')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "679892d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sapphire import download_data\n",
    "if '/events' not in data:\n",
    "    print('downloading events')\n",
    "    download_data(data, '/', station, start, end)\n",
    "else:\n",
    "    print('events already downloaded.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e24a528f",
   "metadata": {},
   "outputs": [],
   "source": [
    "events = data.root.events\n",
    "events"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56e9332c",
   "metadata": {},
   "source": [
    "## Bepalen van tijdsverschil"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2ce33c",
   "metadata": {},
   "outputs": [],
   "source": [
    "t1 = events.col('t1')\n",
    "t3 = events.col('t3')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d430f247",
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = t3 - t1\n",
    "_ = plt.hist(dt, bins = 50, histtype='step')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "602e417a",
   "metadata": {},
   "source": [
    "In het histogram valt op dat er *vreemde* pieken te zien zijn rond dt = 0 en dt\n",
    "= 1000 en dt = -1000.\n",
    "\n",
    "Dit komt omdat er veel events zijn waarbij detector 1 en 3 geen deeltjes hebben\n",
    "detecteerd, of niet allebei. Als er geen aankomsttijd van een deeltje\n",
    "gereconstrueerd kon worden is de aankomsttijd `t = -999`. Allereerst filteren we\n",
    "de tijdsverschillen daarom op `t > 0.`. Vervolgens filteren we de data op\n",
    "`pulshoogte > 200 (ADC)`, zodat er we alleen events overhouden waarbij de\n",
    "deeltjes dichtheid in de detector voldoende groot is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5aac619d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ph = events.col('pulseheights')\n",
    "ph1 = ph[:, 0]\n",
    "ph3 = ph[:, 2]\n",
    "dt_1_3 = dt.compress((t1 > 0.) & (t3 > 0.) & (ph1 > 200.) & (ph3 > 200.))\n",
    "_ = plt.hist(dt_1_3, bins=50, histtype='step')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b952d487",
   "metadata": {},
   "source": [
    "We zien `data` rond dt = 0 en ruis daarom heen. We plotten nu alleen\n",
    "tijdsverschillen kleiner dan ongeveer 100 ns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad8c4800",
   "metadata": {},
   "outputs": [],
   "source": [
    "bins=np.arange(-100, 100, 5)\n",
    "n, _, _ = plt.hist(dt_1_3, bins=bins, histtype='step')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bb394b6",
   "metadata": {},
   "source": [
    "Het tijdsverschillen zijn ongeveer normaal verdeeld, maar het gemiddelde is niet\n",
    "0. De verschuiving van het gemiddelde is de `tijd offset`, of hier 'detector\n",
    "tijd offset'\n",
    "\n",
    "## Detector tijd offset bepalen uit het histogram\n",
    "\n",
    "Het gemiddelde is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cb68547",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(dt_1_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd20e9c",
   "metadata": {},
   "source": [
    "We kunnen ook een normale verdeling fitten:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46648017",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss(x, N, mu, sigma):\n",
    "    \"\"\"Gaussian distribution\"\"\"\n",
    "    return N * np.exp(-0.5 * ((x - mu) / sigma) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "499d44c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "x = (bins[:-1] + bins[1:]) / 2\n",
    "popt, pcov = curve_fit(gauss, x, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24728bd4",
   "metadata": {},
   "source": [
    "De parameter `popt` bevat de parameters uit de fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "965e9027",
   "metadata": {},
   "outputs": [],
   "source": [
    "N, mu, sigma = popt\n",
    "mu, sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4889542f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(dt_1_3, bins=bins, histtype='step')\n",
    "plt.plot(bins, gauss(bins, N, mu, sigma), 'r--')\n",
    "plt.legend(['fit','dt'])\n",
    "plt.xlabel('t3 - t1 [ns]')\n",
    "plt.ylabel('counts')\n",
    "plt.xlim((-60, 70))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22775bdb",
   "metadata": {},
   "source": [
    "## Detector tijd offsets ophalen mbv de API\n",
    "\n",
    "De publieke database (https://data.hisparc.nl) bepaalt dagelijks de detector tijd\n",
    "offsets voor alle stations.\n",
    "\n",
    "De informatie is via de API beschikbaar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c13ae46",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sapphire import Station\n",
    "s501 = Station(501)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6de29ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sapphire.transformations.clock import datetime_to_gps\n",
    "ts = datetime_to_gps(datetime(2016,5,1))\n",
    "ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "970ad91b",
   "metadata": {},
   "outputs": [],
   "source": [
    "s501.detector_timing_offset(ts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2d35db6",
   "metadata": {},
   "outputs": [],
   "source": [
    "o1, o2, o3, o4 = s501.detector_timing_offset(ts)\n",
    "print(\"De detector tijd offset tussen detector 3 en 1 op is: %2.1f ns (t = %d)\" % ((o3 - o1), ts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "350949c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "data.close()"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}