docs/tutorials/theta_square.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make a theta-square plot\n",
    "\n",
    "This is a basic example to analyze some events and make a $\\Theta^2$ plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:35.515499Z",
     "start_time": "2018-06-15T12:49:34.968051Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:37.807612Z",
     "start_time": "2018-06-15T12:49:35.520552Z"
    }
   },
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "from astropy.coordinates.angle_utilities import angular_separation\n",
    "from astropy.coordinates import SkyCoord, AltAz\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from ctapipe.containers import ImageParametersContainer\n",
    "from ctapipe.io import EventSource\n",
    "from ctapipe.visualization import CameraDisplay\n",
    "from ctapipe.instrument import CameraGeometry\n",
    "from ctapipe.calib import CameraCalibrator\n",
    "from ctapipe.reco import HillasReconstructor\n",
    "from ctapipe.image import hillas_parameters, tailcuts_clean\n",
    "from ctapipe.utils import datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:37.887391Z",
     "start_time": "2018-06-15T12:49:37.818824Z"
    }
   },
   "source": [
    "Get source events in MC dataset. Here we stop at 10 events, just to make this example run fast, but for real use, one would need more statistics. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:37.887391Z",
     "start_time": "2018-06-15T12:49:37.818824Z"
    }
   },
   "outputs": [],
   "source": [
    "filename = datasets.get_dataset_path(\"gamma_test_large.simtel.gz\")\n",
    "source = EventSource(filename, allowed_tels={1, 2, 3, 4}, max_events=10)\n",
    "subarray = source.subarray\n",
    "reco = HillasReconstructor(subarray=subarray)\n",
    "calib = CameraCalibrator(subarray=source.subarray)\n",
    "horizon_frame = AltAz()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:47.500199Z",
     "start_time": "2018-06-15T12:49:37.893169Z"
    }
   },
   "outputs": [],
   "source": [
    "off_angles = []\n",
    "for event in source:\n",
    "\n",
    "    # calibrating the event\n",
    "    calib(event)\n",
    "    \n",
    "    hillas_params = {}\n",
    "    \n",
    "    # pointing direction of the telescopes\n",
    "    telescope_pointings = {}\n",
    "    \n",
    "    subarray = source.subarray\n",
    "    \n",
    "    # get hillas params for each event in different telescopes\n",
    "    for tel_id in event.dl0.tel.keys():\n",
    "        \n",
    "        event.dl1.tel[tel_id].parameters = ImageParametersContainer()\n",
    "\n",
    "        # telescope pointing direction\n",
    "        telescope_pointings[tel_id] = SkyCoord(\n",
    "            alt=event.pointing.tel[tel_id].altitude,\n",
    "            az=event.pointing.tel[tel_id].azimuth,\n",
    "            frame=horizon_frame\n",
    "        )  \n",
    "\n",
    "        # Camera Geometry required for hillas parametrization\n",
    "        camgeom = subarray.tel[tel_id].camera.geometry\n",
    "\n",
    "        # note the [0] is for channel 0 which is high-gain channel\n",
    "        image = event.dl1.tel[tel_id].image\n",
    "\n",
    "        # Cleaning  of the image\n",
    "        cleaned_image = image\n",
    "        # create a clean mask of pixels above the threshold\n",
    "        cleanmask = tailcuts_clean(\n",
    "            camgeom, image, picture_thresh=10, boundary_thresh=5\n",
    "        )\n",
    "        # set all rejected pixels to zero\n",
    "        cleaned_image[~cleanmask] = 0\n",
    "\n",
    "        # Calulate hillas parameters\n",
    "        # It fails for empty pixels\n",
    "        try:\n",
    "            hillas = hillas_parameters(camgeom, cleaned_image)\n",
    "        except Exception as e:\n",
    "            print(e)\n",
    "            pass\n",
    "        \n",
    "        if hillas.width.value > 0:\n",
    "            hillas_params[tel_id] = hillas\n",
    "\n",
    "    if len(hillas_params) < 2:\n",
    "        continue\n",
    "\n",
    "    array_pointing = SkyCoord(\n",
    "        az=event.pointing.array_azimuth,\n",
    "        alt=event.pointing.array_altitude,\n",
    "        frame=horizon_frame\n",
    "    )\n",
    "    \n",
    "    reco = HillasReconstructor(subarray=subarray)\n",
    "    reco_result = reco._predict(\n",
    "        event,\n",
    "        hillas_params,\n",
    "        source.subarray,\n",
    "        array_pointing,\n",
    "        telescope_pointings,\n",
    "    )\n",
    "\n",
    "    # get angular offset between reconstructed shower direction and MC\n",
    "    # generated shower direction\n",
    "    mcshower = event.simulation.shower\n",
    "    off_angle = angular_separation(mcshower.az, mcshower.alt, reco_result.az, reco_result.alt)\n",
    "\n",
    "    # Appending all estimated off angles\n",
    "    off_angles.append(off_angle.to(u.deg).value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:47.507369Z",
     "start_time": "2018-06-15T12:49:47.502642Z"
    }
   },
   "source": [
    "calculate theta square for angles which are not nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:47.507369Z",
     "start_time": "2018-06-15T12:49:47.502642Z"
    }
   },
   "outputs": [],
   "source": [
    "off_angles = np.array(off_angles)\n",
    "thetasquare = off_angles[np.isfinite(off_angles)]**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:48.264122Z",
     "start_time": "2018-06-15T12:49:47.511172Z"
    }
   },
   "source": [
    "## Plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-06-15T12:49:48.264122Z",
     "start_time": "2018-06-15T12:49:47.511172Z"
    }
   },
   "outputs": [],
   "source": [
    "plt.hist(thetasquare, bins=10, range=[0,0.4])\n",
    "plt.xlabel(r'$\\theta^2$ (deg)')\n",
    "plt.ylabel(\"# of events\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "again, this plot is not beautiful since we have such low stats"
   ]
  }
 ],
 "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"
  },
  "toc": {
   "nav_menu": {
    "height": "13px",
    "width": "253px"
   },
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}