docs/tutorials/ctapipe_handson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Started with ctapipe\n",
    "\n",
    "This hands-on was presented at the Paris CTA Consoritum meeting (K. Kosack)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: load and loop over data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctapipe.io import EventSource\n",
    "from ctapipe import utils\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = utils.get_dataset_path(\"gamma_prod5.simtel.zst\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = EventSource(path, max_events=5)\n",
    "\n",
    "for event in source:\n",
    "    print(event.count, event.index.event_id, event.simulation.shower.energy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event.r1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for event in EventSource(path, max_events=5):\n",
    "    print(event.count, event.r1.tel.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event.r0.tel[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0tel = event.r0.tel[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0tel.waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r0tel.waveform.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "note that this is ($N_{channels}$, $N_{pixels}$, $N_{samples}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolormesh(r0tel.waveform[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "brightest_pixel = np.argmax(r0tel.waveform[0].sum(axis=1))\n",
    "print(f\"pixel {brightest_pixel} has sum {r0tel.waveform[0,1535].sum()}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(r0tel.waveform[0,brightest_pixel], label=\"channel 0 (high-gain)\")\n",
    "plt.plot(r0tel.waveform[1,brightest_pixel], label=\"channel 1 (low-gain)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "@interact\n",
    "def view_waveform(chan=0, pix_id=brightest_pixel):\n",
    "    plt.plot(r0tel.waveform[chan, pix_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "try making this compare 2 waveforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Explore the instrument description\n",
    "This is all well and good, but we don't really know what camera or telescope this is... how do we get instrumental description info?\n",
    "\n",
    "Currently this is returned *inside* the event (it will soon change to be separate in next version or so)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray = source.subarray "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.peek()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.tel[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.tel[2].camera"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.tel[2].optics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tel = subarray.tel[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tel.camera"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "tel.optics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tel.camera.geometry.pix_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tel.camera.geometry.to_table()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tel.optics.mirror_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctapipe.visualization import CameraDisplay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = CameraDisplay(tel.camera.geometry)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = CameraDisplay(tel.camera.geometry)\n",
    "disp.image = r0tel.waveform[0,:,10]  # display channel 0, sample 0 (try others like 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " ** aside: ** show demo using a CameraDisplay in interactive mode in ipython rather than notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Apply some calibration and trace integration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctapipe.calib import CameraCalibrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calib = CameraCalibrator(subarray=subarray)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for event in EventSource(path, max_events=5):\n",
    "    calib(event) # fills in r1, dl0, and dl1\n",
    "    print(event.dl1.tel.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event.dl1.tel[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl1tel = event.dl1.tel[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl1tel.image.shape # note this will be gain-selected in next version, so will be just 1D array of 1855"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dl1tel.peak_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CameraDisplay(tel.camera.geometry, image=dl1tel.image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for Hillas Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctapipe.image import hillas_parameters, tailcuts_clean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = dl1tel.image\n",
    "mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5)\n",
    "mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CameraDisplay(tel.camera.geometry, image=mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cleaned = image.copy()\n",
    "cleaned[~mask] = 0 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = CameraDisplay(tel.camera.geometry, image=cleaned)\n",
    "disp.cmap = plt.cm.coolwarm\n",
    "disp.add_colorbar()\n",
    "plt.xlim(0.5, 1.0)\n",
    "plt.ylim(-1.0, 0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params = hillas_parameters(tel.camera.geometry, cleaned)\n",
    "print(params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "disp = CameraDisplay(tel.camera.geometry, image=cleaned)\n",
    "disp.cmap = plt.cm.coolwarm\n",
    "disp.add_colorbar()\n",
    "plt.xlim(0.5, 1.0)\n",
    "plt.ylim(-1.0, 0.0)\n",
    "disp.overlay_moments(params, color='white', lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4:  Let's put it all together: \n",
    "- loop over events, selecting only telescopes of the same type (e.g. LST:LSTCam)\n",
    "- for each event, apply calibration/trace integration\n",
    "- calculate Hillas parameters \n",
    "- write out all hillas paremeters to a file that can be loaded with Pandas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "first let's select only those telescopes with LST:LSTCam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.telescope_types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subarray.get_tel_ids_for_type(\"LST_LST_LSTCam\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's write out program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = utils.get_dataset_path(\"gamma_prod5.simtel.zst\") \n",
    "source = EventSource(data) # remove the max_events limit to get more stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for event in source:\n",
    "    calib(event)\n",
    "    \n",
    "    for tel_id, tel_data in event.dl1.tel.items():\n",
    "        tel = source.subarray.tel[tel_id]\n",
    "        mask = tailcuts_clean(tel.camera.geometry, tel_data.image)\n",
    "        if np.count_nonzero(mask) > 0:\n",
    "            params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctapipe.io import HDF5TableWriter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with HDF5TableWriter(filename='hillas.h5', group_name='dl1', overwrite=True) as writer:\n",
    "    \n",
    "    source = EventSource(data, allowed_tels=[1,2,3,4],  max_events=10)\n",
    "    for event in source:\n",
    "        calib(event)\n",
    "    \n",
    "        for tel_id, tel_data in event.dl1.tel.items():\n",
    "            tel = source.subarray.tel[tel_id]\n",
    "            mask = tailcuts_clean(tel.camera.geometry, tel_data.image)\n",
    "            params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])\n",
    "            writer.write(\"hillas\", params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We can now load in the file we created and plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls *.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "\n",
    "hillas = pd.read_hdf(\"hillas.h5\", key='/dl1/hillas')\n",
    "hillas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = hillas.hist(figsize=(8,8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do this yourself, chose a larger file to loop over more events to get better statistics"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}