notebook_starting_kit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> <a href=http://www.datascience-paris-saclay.fr/>Paris Saclay Center for Data Science</a> </h1>\n",
    "\n",
    "<h2> RAMP on charged particle tracking in 2D with a possible future LHC Silicon detector </h2>\n",
    "\n",
    "<i> Thomas Boser (CHALEARN), Isabelle Guyon (LRI UPSud/INRIA), Mikhail Hushchyn (YSDA/Russia), Balázs Kégl (LAL/Orsay), David Rousseau (LAL/Orsay), Yetkin Yılmaz (LAL/Orsay) </i>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2> Introduction </h2>\n",
    "\n",
    "Tracking is one of the most important tasks in a high-energy physics experiment, as it provides high-precision position and momentum information belonging to charged particles. Such information is crucial for a diverse variety of physics studies - from Standard Model tests to new particle searches - which requires robust low-level optimization without information loss which can be further refined for a narrower and more specific physics context.\n",
    "\n",
    "Through the history of high-energy physics, there have been many different types of tracking detectors with very different design principles: from bubble chambers to time-projection chambers, from proportional counters to spark chambers... Although each of these have provided a different data topology, they all relied on simple basics: the small energy deposit of particles in well-defined locations, with particles bending in an externally applied magnetic field.\n",
    "\n",
    "In this challenge, we focus on the topology of silicon detectors in which there are few locations along the polar axis, however a very high precision along the azimuth. Such topology helps to reduce the tracking problem to layer-by-layer azimuth determination, however, we hope to open room for further innovation as well.\n",
    "\n",
    "\n",
    "<h2> Quick Glossary </h2>\n",
    "\n",
    "We use the following vocabulary in the context of this challenge. It may be slightly different from the general high-energy physics context. Some terms are not used in the workbook, but kept here in case they are used in the discussions.\n",
    "\n",
    "**event**: a recorded (or simulated) collision in which many particles are present, the basic unit of particle physics data\n",
    "\n",
    "**pixel**: the smallest detector unit\n",
    "\n",
    "**hit**: a pixel through which a particle has passed and left signal in a given event\n",
    "\n",
    "**cluster**: a set of hits, belonging (or predicted to be belonging) to the trajectory of a single particle\n",
    "\n",
    "**reconstruct**: same thing as predict, but may also refer to further derived quantities\n",
    "\n",
    "**track**: a reconstructed particle, that is, a cluster but also including additional derived information, such as the overall curvature and angle\n",
    "\n",
    "**vertex**: the point close to the center of the detector, from which the particles have originated\n",
    "\n",
    "**impact parameter**: the shortest distance between a track and the origin\n",
    "\n",
    "**(momentum/angle/etc) resolution** : width of the normal distribution of the difference between a predicted quantity and the true value of that quantity\n",
    "\n",
    "<h2> Objectives </h2>\n",
    "\n",
    "The main objective of the challenge is the optimal matching of the hits in an event. The positions of the hits in the detector are provided as input data, and the user is expected to implement a clustering scheme (trained on data if needed), so that every hit is assigned to a cluster_id.\n",
    "\n",
    "The value of the cluster_id itself is not relevant for the task, what is relevant is which hits are clustered together, and whether this clustering corresponds well to the input particles. The score function that describes this is included in the notebook, and details will be mentioned there.\n",
    "\n",
    "<h2> Application </h2>\n",
    "\n",
    "The user is expected to implement the class *clusterer.py*, which contains the __init__, fit, and predict_single_event functions.\n",
    "\n",
    "**\\__init\\__** is where parameters should be set.\n",
    "\n",
    "**fit** is the training function (not to be confused with track-fitting), where the algorithm has access to the ground-truth. This function is to be run once on an input array that contains a set of training events. The user is able to implement any event-level or particle-level segmentation of the input array in order to set up the training in any desired way.\n",
    "\n",
    "**predict_single_event** is the function to reconstruct the hit clusters (tracks), returning an array of predicted (reconstructed) ids associated to each hit in the input array. This function takes only the hits from a single event as input, with the event_id dropped, and the RAMP machinery takes care of running this function on each event.\n",
    "\n",
    "The **task** in the challenge is to implement this class in a way that the predict_single_event function returns a numpy array of assigned cluster_ids. At any level of this task, machine-learning techniques can be employed for sub-tasks defined by the user.\n",
    "\n",
    "<h2> Detector </h2>\n",
    "\n",
    "Image from the Atlas Experiment:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://atlasexperiment.org/photos/atlas_photos/selected-photos/events/atlas2009-collision-atlantis-141749-405315-hits-web.png\" alt=\"Drawing\" style=\"width: 700px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data provided to the user is a list of hit positions from a simple toy detector model that mimics the Atlas detector design (which is generic enough for recent silicon-based tracking detectors). The detector has an onion-like geometry with 9 layers surrounding the origin with polar distances of R = [39,85,155,213,271,405,562,762,1000] cm.\n",
    "These layers have a very small thickness compared to the distances, therefore the thickness can be neglected.\n",
    "\n",
    "Each layer is segmented in azimuth with high granularity. There are ($2\\pi$R/pitch)+1 pixels in every layer, where pitch is 0.025 cm for layers 0-4 and 0.05 cm for layers 5-9.\n",
    "\n",
    "Every \"pixel\" corresponds to the smallest detector unit defined by **layer** and **iphi** (azimuthal index).\n",
    "\n",
    "<h2> Simulation </h2>\n",
    "\n",
    "The challenge uses a toy model for particle generation and simulation, in which a Poisson distribution is sampled to determine the number of particles in each event, with an average of 10 particles per event.\n",
    "\n",
    "The particles are sampled uniformly in azimuth and momentum, with bounds on the momentum. Each particle originates from a vertex that is also randomly sampled from a narrow normal distribution around the origin. The proper dimensions of the momentum and position and determination of these values for the tracks are beyond the scope of the challenge.\n",
    "\n",
    "The particles generated this way are simulated in a uniform magnetic field. The detector material is considered to cause multiple-scattering, and this is implemented as a random rotation of the particle momentum at every detector position, sampled from a narrow normal distribution that roughly corresponds to the material of the Atlas tracker.\n",
    "\n",
    "In addition, some hit inefficiency is simulated by a random drop of some hits (with 3% probability), and a particle stopping probability of 1% is applied at each layer to simulate effects of hadronic interactions. Keeping these in mind, the algorithms might be desired to be able to handle cases when the particle doesn't have a hit on every layer.\n",
    "\n",
    "Since the detector has a very high granularity in azimuth, the cases where two particles pass through a single pixel are neglected (less than 0.2% probability).\n",
    "\n",
    "<h2> Submission </h2>\n",
    "\n",
    "Submissions are to be made in the RAMP sandbox:\n",
    "http://ramp.studio/events/HEP_tracking/sandbox\n",
    "\n",
    "You can choose to either edit the code here directly, or upload your local version of the clusterer.py. You cannot upload any additional files. All code, (additional classes, functions, or constants etc) must be written inside the clusterer.py.\n",
    "\n",
    "Don't forget to save and submit. You will make several submissions, make sure they are named in a later recognizable way.\n",
    "\n",
    "Submission are run on a private set of events, different from the one provided, hence slightly different score is expected. You’re supposed to develop your software on your own laptop/platform and infrequently submit to RAMP platform (once you submitted on RAMP you have to wait 15 min before the next submission). Your submission will run in the backend within a few hours, and you will be able to see the score on the larger sample.\n",
    "\n",
    "It is very <b><span style=\"color:red\">important that you test your submission files before submitting them</span></b>, for this, please see the **unit test** instructions at the very end of this notebook.\n",
    "\n",
    "<h2> Beyond the challenge </h2>\n",
    "\n",
    "The following aspects of tracking fall out of scope of the challenge:\n",
    "\n",
    "* Track fitting\n",
    "* Particle efficiency\n",
    "* Fake tracks\n",
    "* Momentum resolution\n",
    "* Vertex finding and impact parameter resolution\n",
    "* Tracking in 3D\n",
    "\n",
    "Fake combinatorial tracks do affect the score in an indirect way, by reducing the efficiency due to one-to-one assignment between true and predicted cluster_ids. If there are fake tracks - that is, clusters that contain hits belonging to many different particles - that means the good cluster attached to each of these particles will be missing some of the hits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "\n",
    "* numpy>=1.10.0  \n",
    "* matplotlib>=1.5.0 \n",
    "* pandas>=0.17.0  \n",
    "* scikit-learn>=0.17 (different syntaxes for v0.17 and v0.18)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploratory data analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = \"public_train_hacked.csv\"\n",
    "df = pd.read_csv(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data consists of the following columns:\n",
    "event_id, cluster_id, layer, iphi, x, y\n",
    "\n",
    "**event_id** (integer): The event number; particles that have the same event_id belong to the same event. Notice that the event ids are not consecutive since this is a randomly sampled subset of the full sample.\n",
    "**cluster_id** (integer): The true cluster_id corresponds to the particle that created the hits. The predicted cluster_id will correspond to the label of the cluster that is reconstructed. The value itself of the reconstructed cluster is irrelevant, only how the hits are grouped with those ids matters. The only requirement on the predicted cluster_ids is that they are non-negative for properly reconstructed clusters. A value of -1 means that hit is not associated to any cluster.\n",
    "**layer** (integer): The detector layer that the hit sits in, representing how far the hit is from the origin. The distances of the various layers were mentioned above in *Detector*.\n",
    "**iphi** (integer): The azimuthal id of the pixel along a given layer. The phi position of the hit can be determined form this value based on how many segmentations that layer has, which was described in *Detector*.\n",
    "**x** (float): The x position of the hit, in cm.\n",
    "**y** (float): The y position of the hit, in cm.\n",
    "\n",
    "Notice that the information provided in data is redundant: the **layer & iphi** and **x & y** values are translatable. It is provided this way for the sake of convenience, so be careful when you need to feed algorithms with independent variables.\n",
    "\n",
    "The array indices of the above variables may change depending on whether the event_id has been dropped or not, therefore make sure you are careful in using the index for the variable you want to pick up."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>event_id</th>\n",
       "      <th>cluster_id</th>\n",
       "      <th>layer</th>\n",
       "      <th>iphi</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>66</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1443</td>\n",
       "      <td>23.472228</td>\n",
       "      <td>31.145698</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>3122</td>\n",
       "      <td>51.614467</td>\n",
       "      <td>67.534782</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>5666</td>\n",
       "      <td>94.656900</td>\n",
       "      <td>122.739852</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>70</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>7764</td>\n",
       "      <td>130.515322</td>\n",
       "      <td>168.329292</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>53</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>9854</td>\n",
       "      <td>166.531067</td>\n",
       "      <td>213.795238</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>60</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>5</td>\n",
       "      <td>7324</td>\n",
       "      <td>250.419293</td>\n",
       "      <td>318.300452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>10102</td>\n",
       "      <td>349.898441</td>\n",
       "      <td>439.789815</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>7</td>\n",
       "      <td>13594</td>\n",
       "      <td>478.431276</td>\n",
       "      <td>593.083059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>8</td>\n",
       "      <td>17680</td>\n",
       "      <td>634.064698</td>\n",
       "      <td>773.280000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>9295</td>\n",
       "      <td>36.958474</td>\n",
       "      <td>-12.452760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>20281</td>\n",
       "      <td>80.732137</td>\n",
       "      <td>-26.595526</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>37004</td>\n",
       "      <td>147.381280</td>\n",
       "      <td>-47.997483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>69</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "      <td>50865</td>\n",
       "      <td>202.641758</td>\n",
       "      <td>-65.614921</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "      <td>64730</td>\n",
       "      <td>257.932600</td>\n",
       "      <td>-83.137078</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>5</td>\n",
       "      <td>48390</td>\n",
       "      <td>385.801871</td>\n",
       "      <td>-123.214919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>6</td>\n",
       "      <td>67181</td>\n",
       "      <td>535.840039</td>\n",
       "      <td>-169.468147</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>7</td>\n",
       "      <td>91140</td>\n",
       "      <td>727.313411</td>\n",
       "      <td>-227.286608</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>8</td>\n",
       "      <td>119688</td>\n",
       "      <td>955.690631</td>\n",
       "      <td>-294.372923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0</td>\n",
       "      <td>2332</td>\n",
       "      <td>2.959582</td>\n",
       "      <td>38.887541</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>5093</td>\n",
       "      <td>6.188221</td>\n",
       "      <td>84.774441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>9286</td>\n",
       "      <td>11.314854</td>\n",
       "      <td>154.586462</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "      <td>12754</td>\n",
       "      <td>15.716876</td>\n",
       "      <td>212.419349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>16216</td>\n",
       "      <td>20.268475</td>\n",
       "      <td>270.240983</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "      <td>12094</td>\n",
       "      <td>31.443203</td>\n",
       "      <td>403.777569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>64</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>6</td>\n",
       "      <td>16744</td>\n",
       "      <td>45.549340</td>\n",
       "      <td>560.151102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>7</td>\n",
       "      <td>22631</td>\n",
       "      <td>65.319563</td>\n",
       "      <td>759.195202</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>8</td>\n",
       "      <td>29586</td>\n",
       "      <td>91.372164</td>\n",
       "      <td>995.816814</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>9322</td>\n",
       "      <td>37.168451</td>\n",
       "      <td>-11.811276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>54</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>20301</td>\n",
       "      <td>80.887183</td>\n",
       "      <td>-26.120178</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>2</td>\n",
       "      <td>37015</td>\n",
       "      <td>147.466204</td>\n",
       "      <td>-47.735927</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559362</th>\n",
       "      <td>90099</td>\n",
       "      <td>3</td>\n",
       "      <td>5</td>\n",
       "      <td>37988</td>\n",
       "      <td>-9.124192</td>\n",
       "      <td>-404.897208</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559361</th>\n",
       "      <td>90099</td>\n",
       "      <td>3</td>\n",
       "      <td>6</td>\n",
       "      <td>52776</td>\n",
       "      <td>-9.599398</td>\n",
       "      <td>-561.918011</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559377</th>\n",
       "      <td>90099</td>\n",
       "      <td>3</td>\n",
       "      <td>7</td>\n",
       "      <td>71656</td>\n",
       "      <td>-8.049829</td>\n",
       "      <td>-761.957479</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559390</th>\n",
       "      <td>90099</td>\n",
       "      <td>3</td>\n",
       "      <td>8</td>\n",
       "      <td>94191</td>\n",
       "      <td>-2.849989</td>\n",
       "      <td>-999.995939</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559398</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>6720</td>\n",
       "      <td>-15.359496</td>\n",
       "      <td>-35.848094</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559381</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>14632</td>\n",
       "      <td>-33.795538</td>\n",
       "      <td>-77.992702</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559371</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>26620</td>\n",
       "      <td>-63.042657</td>\n",
       "      <td>-141.600224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559394</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>3</td>\n",
       "      <td>36497</td>\n",
       "      <td>-88.546447</td>\n",
       "      <td>-193.722809</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559388</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>4</td>\n",
       "      <td>46333</td>\n",
       "      <td>-114.974331</td>\n",
       "      <td>-245.401514</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559350</th>\n",
       "      <td>90099</td>\n",
       "      <td>4</td>\n",
       "      <td>5</td>\n",
       "      <td>34433</td>\n",
       "      <td>-180.313378</td>\n",
       "      <td>-362.645951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559374</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>0</td>\n",
       "      <td>7809</td>\n",
       "      <td>11.273995</td>\n",
       "      <td>-37.334930</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559375</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>17017</td>\n",
       "      <td>24.515285</td>\n",
       "      <td>-81.387965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559370</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>2</td>\n",
       "      <td>31063</td>\n",
       "      <td>45.470863</td>\n",
       "      <td>-148.180298</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559364</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>3</td>\n",
       "      <td>42729</td>\n",
       "      <td>63.500558</td>\n",
       "      <td>-203.314237</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559402</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>54416</td>\n",
       "      <td>82.029876</td>\n",
       "      <td>-258.286855</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559387</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "      <td>40753</td>\n",
       "      <td>126.948003</td>\n",
       "      <td>-384.589657</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559365</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>6</td>\n",
       "      <td>56697</td>\n",
       "      <td>183.045992</td>\n",
       "      <td>-531.355027</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559348</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>7</td>\n",
       "      <td>77125</td>\n",
       "      <td>260.065898</td>\n",
       "      <td>-716.246975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559396</th>\n",
       "      <td>90099</td>\n",
       "      <td>5</td>\n",
       "      <td>8</td>\n",
       "      <td>101613</td>\n",
       "      <td>359.982503</td>\n",
       "      <td>-932.959055</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559363</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>3229</td>\n",
       "      <td>-18.664275</td>\n",
       "      <td>34.243902</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559366</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>1</td>\n",
       "      <td>6999</td>\n",
       "      <td>-39.831853</td>\n",
       "      <td>75.089436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559355</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>2</td>\n",
       "      <td>12723</td>\n",
       "      <td>-71.752685</td>\n",
       "      <td>137.391965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559359</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>3</td>\n",
       "      <td>17454</td>\n",
       "      <td>-97.940287</td>\n",
       "      <td>189.147298</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559400</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>22178</td>\n",
       "      <td>-123.971610</td>\n",
       "      <td>240.981410</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559389</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>5</td>\n",
       "      <td>16527</td>\n",
       "      <td>-183.262241</td>\n",
       "      <td>361.164715</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559382</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>6</td>\n",
       "      <td>22867</td>\n",
       "      <td>-251.312901</td>\n",
       "      <td>502.678651</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559379</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>7</td>\n",
       "      <td>30896</td>\n",
       "      <td>-335.893091</td>\n",
       "      <td>683.973560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559369</th>\n",
       "      <td>90099</td>\n",
       "      <td>6</td>\n",
       "      <td>8</td>\n",
       "      <td>40379</td>\n",
       "      <td>-433.298019</td>\n",
       "      <td>901.250701</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559378</th>\n",
       "      <td>90099</td>\n",
       "      <td>7</td>\n",
       "      <td>0</td>\n",
       "      <td>2857</td>\n",
       "      <td>-10.047652</td>\n",
       "      <td>37.683480</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7559376</th>\n",
       "      <td>90099</td>\n",
       "      <td>7</td>\n",
       "      <td>1</td>\n",
       "      <td>6192</td>\n",
       "      <td>-21.059450</td>\n",
       "      <td>82.349861</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>7559409 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         event_id  cluster_id  layer    iphi           x           y\n",
       "66              0           0      0    1443   23.472228   31.145698\n",
       "67              0           0      1    3122   51.614467   67.534782\n",
       "16              0           0      2    5666   94.656900  122.739852\n",
       "70              0           0      3    7764  130.515322  168.329292\n",
       "53              0           0      4    9854  166.531067  213.795238\n",
       "60              0           0      5    7324  250.419293  318.300452\n",
       "59              0           0      6   10102  349.898441  439.789815\n",
       "46              0           0      7   13594  478.431276  593.083059\n",
       "44              0           0      8   17680  634.064698  773.280000\n",
       "49              0           1      0    9295   36.958474  -12.452760\n",
       "17              0           1      1   20281   80.732137  -26.595526\n",
       "30              0           1      2   37004  147.381280  -47.997483\n",
       "69              0           1      3   50865  202.641758  -65.614921\n",
       "10              0           1      4   64730  257.932600  -83.137078\n",
       "22              0           1      5   48390  385.801871 -123.214919\n",
       "0               0           1      6   67181  535.840039 -169.468147\n",
       "3               0           1      7   91140  727.313411 -227.286608\n",
       "11              0           1      8  119688  955.690631 -294.372923\n",
       "72              0           2      0    2332    2.959582   38.887541\n",
       "35              0           2      1    5093    6.188221   84.774441\n",
       "43              0           2      2    9286   11.314854  154.586462\n",
       "40              0           2      3   12754   15.716876  212.419349\n",
       "4               0           2      4   16216   20.268475  270.240983\n",
       "15              0           2      5   12094   31.443203  403.777569\n",
       "64              0           2      6   16744   45.549340  560.151102\n",
       "1               0           2      7   22631   65.319563  759.195202\n",
       "52              0           2      8   29586   91.372164  995.816814\n",
       "39              0           3      0    9322   37.168451  -11.811276\n",
       "54              0           3      1   20301   80.887183  -26.120178\n",
       "5               0           3      2   37015  147.466204  -47.735927\n",
       "...           ...         ...    ...     ...         ...         ...\n",
       "7559362     90099           3      5   37988   -9.124192 -404.897208\n",
       "7559361     90099           3      6   52776   -9.599398 -561.918011\n",
       "7559377     90099           3      7   71656   -8.049829 -761.957479\n",
       "7559390     90099           3      8   94191   -2.849989 -999.995939\n",
       "7559398     90099           4      0    6720  -15.359496  -35.848094\n",
       "7559381     90099           4      1   14632  -33.795538  -77.992702\n",
       "7559371     90099           4      2   26620  -63.042657 -141.600224\n",
       "7559394     90099           4      3   36497  -88.546447 -193.722809\n",
       "7559388     90099           4      4   46333 -114.974331 -245.401514\n",
       "7559350     90099           4      5   34433 -180.313378 -362.645951\n",
       "7559374     90099           5      0    7809   11.273995  -37.334930\n",
       "7559375     90099           5      1   17017   24.515285  -81.387965\n",
       "7559370     90099           5      2   31063   45.470863 -148.180298\n",
       "7559364     90099           5      3   42729   63.500558 -203.314237\n",
       "7559402     90099           5      4   54416   82.029876 -258.286855\n",
       "7559387     90099           5      5   40753  126.948003 -384.589657\n",
       "7559365     90099           5      6   56697  183.045992 -531.355027\n",
       "7559348     90099           5      7   77125  260.065898 -716.246975\n",
       "7559396     90099           5      8  101613  359.982503 -932.959055\n",
       "7559363     90099           6      0    3229  -18.664275   34.243902\n",
       "7559366     90099           6      1    6999  -39.831853   75.089436\n",
       "7559355     90099           6      2   12723  -71.752685  137.391965\n",
       "7559359     90099           6      3   17454  -97.940287  189.147298\n",
       "7559400     90099           6      4   22178 -123.971610  240.981410\n",
       "7559389     90099           6      5   16527 -183.262241  361.164715\n",
       "7559382     90099           6      6   22867 -251.312901  502.678651\n",
       "7559379     90099           6      7   30896 -335.893091  683.973560\n",
       "7559369     90099           6      8   40379 -433.298019  901.250701\n",
       "7559378     90099           7      0    2857  -10.047652   37.683480\n",
       "7559376     90099           7      1    6192  -21.059450   82.349861\n",
       "\n",
       "[7559409 rows x 6 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.sort_values(by=['event_id','cluster_id','layer','iphi'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y = df.drop(['layer','iphi','x','y'], axis=1).values\n",
    "X = df.drop(['cluster_id'], axis=1).values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Looping over events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_events = y[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "events = np.unique(y_events)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize the hits in a single event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ievent = 1\n",
    "event_indices = X[:,0]==ievent\n",
    "X_event = X[event_indices,1:] # event_id column dropped\n",
    "pixelx = X_event[:,2]\n",
    "pixely = X_event[:,3]\n",
    "particles = y[event_indices][:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Notice that** the information provided has some redundancy. You may either use **layer and iphi** (integers) or **x and y** (floats)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# event_id column dropped already\n",
    "X_detector = X_event[:,0:2]\n",
    "X_xy = X[:,2:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hits from all particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGAdJREFUeJzt3W2QXNV54PH/g6RhxghYadMIW4YZJTIgb4X1aEtaZ+2K\nBwMy5gMQJwtoPxgvgwsiO2tnd8sI54NUSSrB3mVtjEu24ogYbwVJOCGGbLESqGBclaqsNcGyFZeE\nrKx3BiSDNBsrslgLIcHZD31H6hn1vPfLvd3/X1UXrdO3u8+dHvqZc57znBspJSRJ7e2CZndAktR8\nBgNJksFAkmQwkCRhMJAkYTCQJFGjYBARWyLiSETsrWhbFBHPRsSBiNgZEZdWPPZARByMiP0Rsaai\nfWVE7I2IH0fEl2vRN0nS1Go1Mvgz4CPj2tYDu1JKVwPPAw8ARMR7gduBFcBHgU0REdlzvgb0p5Su\nAq6KiPGvKUmqg5oEg5TS3wDHxjXfCjyW3X8MuC27fwuwLaV0JqU0BBwEVkfE5cDFKaXB7LhvVTxH\nklRH9cwZXJZSOgKQUnoNuCxrXwq8UnHc4axtKXCoov1Q1iZJqrNGJpDd90KScmp+HV/7SEQsSSkd\nyaaAjmbth4ErKo57d9Y2Uft5IsLAIkmzkFKKau21HBlEdhv1NPCJ7P5dwFMV7XdGREdELAOWA7uz\nqaTjEbE6Syh/vOI550kpNey2YcOGhr5fo2+eX3FvrXxunl/tb5OpycggIh4H+oB/HhEvAxuAB4Fv\nR8TdwDDlFUSklPZFxBPAPuA0sC6d6+WngG8CncAzKaUdteifJGlyNQkGKaV/N8FDN0xw/B8Df1yl\n/UXgV2vRJ0nS9FmBPA19fX3N7kJdeX7F1crnBp5fI8VU80h5FBGpiP2WpGaKCFIDEsiSpIIyGEiS\nDAaSJIOBJAmDgSQJg4EkCYOBJAmDgSQJg4FazMjICIODg4yMjDS7K1KhGAzUMrZu3U539zXceON9\ndHdfw9at25vdJakw3I5CLWFkZITu7ms4efIF4FpgL11d1zE8/BKlUqnZ3ZNywe0o1PKGhobo6Oih\nHAgArmXBgm6Ghoaa1ympQAwGagk9PT28+eYQsDdr2cvp08P09PQ0r1NSgRgM1BJKpRJbtmyiq+s6\nLrlkJV1d17FlyyaniKRpMmegljIyMsLQ0BA9PT0GAmmcyXIGBgMVgl/y0tyZQFahuWRUqj9HBso1\nl4xKtePIQIXlklGpMQwGyjWXjEqNYTBQrrlkVGoMcwYqBFcTSXPn0lJJkglkSdLkDAaSJIOBJMlg\nIEnCYCBJwmAgScJgIEnCYCCdNTIywuDgICMjI83uitRwBgMJt8mWrEBW23ObbLULK5ClSbhNtmQw\nkNwmW8JgoDYyUYLYbbIlcwZqE1u3bqe/fx0dHeVRwJYtm1i79o4xx7hNtlqdW1irrZkglspMIKut\nmSCWpmYwUMszQSxNzWCglmeCWJqaOQO1DRPEancmkCVJJpAlSZMzGEiS6h8MImIoIn4YEXsiYnfW\ntigino2IAxGxMyIurTj+gYg4GBH7I2JNvfsnSWrMyOBtoC+l1JtSWp21rQd2pZSuBp4HHgCIiPcC\ntwMrgI8CmyKi6vyWJKl2GhEMosr73Ao8lt1/DLgtu38LsC2ldCalNAQcBFYjSaqrRgSDBDwXEYMR\ncU/WtiSldAQgpfQacFnWvhR4peK5h7M2KXe8MppayfwGvMcHUkqvRkQJeDYiDlAOEJVmvE5048aN\nZ+/39fXR19c3lz5KMzKdje+kZhsYGGBgYGBaxza0ziAiNgCvA/dQziMciYjLgRdSSisiYj2QUkpf\nyI7fAWxIKX1v3OtYZ6CmceM7FVXT6gwi4h0RsTC7fxGwBvh74GngE9lhdwFPZfefBu6MiI6IWAYs\nB3bXs4/STLnxnVpRvaeJlgB/FREpe68/Tyk9GxF/BzwREXcDw5RXEJFS2hcRTwD7gNPAOocAarbx\n21iM3fiuPDJw4zsVndtRSJOYKDcw2r5gQTenTw+bM1AhuDeRNAtT5Qbc+E5FM1kwaMRqIqmQRnMD\nJ0+enxsolUpnb1IrcG8iaQJeFEftxGAgTcCL4qidmDOQpmBuQK3CBLIkyYvbSJImZzCQJBkMJEkG\nA6mu3OZaRWEwkOpk69btdHdfw4033kd39zVs3bq92V2SJuRqIqkO3OZaeeR2FFKDjNYkHDt2bNKt\nLKS8MRhINVK5w+mpUz/h7bcTbnOtonCaSKqBatNCCxZ8gPnzO9zmWrnhNJFUZ9V2OO3qeg/f/vaD\nLFq0yK0slHsGA6kGJrr6WW9vr0FAheDSUqkG3OFURWfOQKohdzhVnrlrqSTJXUslSZMzGEiSDAZS\nM7mRnfLCYCA1iRvZKU9MIEsNNjIywp49e7jttrVuZKeGMoEs5cToaOBjH/sPnDy5mHIggMqN7KRm\ncGQgNcjY/YveCVwNDODIQI3i3kRSDpy/f9HXgF/joouW8/bbh6xYVlM5TSTAVS2NMHb/IoAVdHZ2\n8OST/4Xh4Ze44YYP+xmoaQwGclVLg1Tbv+jRR79Ob28vmzd/gyuvvMrPQE1jzqDNeXnGxqvcv2jX\nrufp71+XJZP/kfLU0Qo/A9WFq4k0odF5bFe1NE6pVGLVqlUAWSB4AThIOZm8Dninn4EazmDQ5s6f\nx/byjI1SLRBDN/Ccn4EazmDQ5tyHv3mqBWI4QGfnp9iyZROACWU1jDkDAY3bh9/9/sfaunU7/f3r\nzl4n+fOf/0/85m/+Bn/5l3/FH/3RQ3R0lAOG109WLXg9A+XC6BefX3BjjU8o3333fbzxxpvA32JS\nX7VkMFDTuWppaud+Ro8ADwEvnn3skktWsmvX5rOJZ2k2XE2kpnPV0tTO/YxuBIYwqa9GMhioIeq9\naqkVKqjP/YxeBTYBfcByk/pqCIOBGqKeq5ZapYJ67M/oC3R2Jv7gD+5mePglcyuqO3MGaqharyZq\nxVzERD8jV2JprswZKDdGq29r9WVWy1xEXqaaqv2MNm/+BldcsZzrr+8v9OhH+WUwUKHVKheR56mm\nzZu/wX33fYZTp5Zx4sRhTp68n/7+dU0PWmotThOp8MYXbs20fqEWU031msIZGRnhiiuu4tSp757t\nG1zHwoXv4vnnH3WpqWbEaSK1tLVr72B4+CV27do8q2TrXKea6jmqKPdtGWP3L1rKqVP/h4ULF9bs\nfSRHBmp7cxkZTPXcuY4Yqr0+vJ+OjiVccME/8eijX3elkabNkYE0ibkse51sVDHZiGG6yerKvi1c\n+D7g/cAFvPnmYt54I7jrrnvMHag2Ukq5ugE3AS8BPwbun+CYJNXa0aNH0+7du9PRo0dn9JyursUJ\nfpggJfhh6upanPbt21e1/ejRo+nxx7elrq7F6dJLV6aursXp8ce3TdmXo0ePpq985SsJusa8Jrwj\n7dy5s2Y/A7W27Luz6ndvrkYGEXEB8FXgI8C/ANZGxDXN7ZXaxWyWvU40qnj99derjhj27Nlz9oI2\nx4+/yMmTL5y3MqjaiKJUKrFkyRJgKWPzB++syblL85vdgXFWAwdTSsMAEbENuJXySEHKpbVr7+CG\nGz48JjcwMjJSseS1PNd/+vQwAB0dPZw8ef600ujzRoNF+Zi99Pdfx89//nM++9nPAW+Oec2Ojv9L\nb29vw89ZrSdvwWAp8ErFvw9RDhBSrpVKpTEjitERQ3//dWOWvPb29lYNEqN1EaM5iMpgMW/eu/jM\nZz6XLS/dT3nPosV0dv6MRx/dbDWyaiJvwWDaNm7cePZ+X18ffX19TeuLVE21EQNQNUiMPja2iG40\nWLxMR8evcOrUtVnbh7noog/y5JPbWLNmDVBOSO/ZsweA3t5eA4QAGBgYYGBgYHoHT5RMaMaN8lKJ\nHRX/Xk+VJDImkFVwkyWrRxPMl1zSm7q6Fqevf/1PxiWjX0gXXnhJ2rdv39njOzouTbA8wTvSggUL\nqyalJSZJIOeqziAi5gEHgOsp7+O7G1ibUto/7riUp35LtTa+PmG0yjqlS3njjdfo6loOHOZLX3qQ\n3/3d9ePqEPro7Ey8/PKPHSFojEJd6SwibgIeplwDsSWl9GCVYwwGajv79++nt/ffjNma4sILP8S8\neUv5xS9+VHHkSi666P/x5JOPsGjRInc51VmTBYPc5QxSSjuAq5vdDylvXn/9dTo7fznLHQBcm11P\n+idU5hhgiFOn3uSWW+6go2MZZ87MfL8mtZ9c1RlIqm5kZIRjx46dt0PrmTMv8/DD/5WOjl8HlgO/\nxgUXvMGZM29z6tR3OXHi+1VrGaTxDAZSzo0Wod1++wOcOfMmHR2/PqbA7d57P8mhQwfZuXMT27f/\nGfPnzwOuorI4bf78K73etCaVu5zBdJgzULuotlFdZ+eHeOqp7VWXkA4ODnL99f2cOHEYOPecCy/8\nEK+8ci6h7FXT2pMb1UkFVW0jvI6OZSxatKjql3hPTw9nzhwG7geuA/4l8H4efviLZ4/P84V81DyO\nDKQcm8322qPLUOfNexenT5dzCvfe+8lZv55aR6FWE0k6Z6JtLSb74p6o8nlkZIRnnnmG+fPHbnZX\nuTeS2pcjAylnqs3nz3WOf3S0MH9+NydOvARsBD5H5cgAMI/Q4gpVdDYdBgO1qtEv7XL9wFBN6gMm\nulrawoW/wltv/ZQtWzYB1Px9lT8GA6kA6jWfPzg4yI033sfx4y+ebbv44l4eeeSz3HzzzQDmEdqE\nq4mkApjsEppzMXYnVBgtVhsNBOU8QnfN31fFYjCQcqLal3bltQ5ma6Krse3a9Tzd3dfwO7/zUJZH\nOP99p3utZhWf00RSk1Umh3ftep7+/nVjVg7Vau6+8n1g/NTQF4GNXHzx1Zw587J5hBZlzkDKqWoJ\n42rLQmutWh5h4cJf5atf/c/mEVqYOQMphyqvd3z8+ItnN5QDWLVqVV2/dKtNSb311k+5+eabKZVK\ndctfKL8MBlKTNPMLd6I8QvXLb0Kt8hfKL6eJpCbJw9YQkxWzjU5h1SN/oeYwZyDlSKMSxrXg7qat\nxWAg5USzEsYSGAykXMjDtJDam6uJpBxohxU6FqkVl8FAapBWX6HjRXOKzWkiqYFadYWOU2DF4MVt\npJyY6MIzRTc6BXbypBfNKSqDgdRgpVKp5b4gx06BlUcGrTQF1g7MGUias6kqmpV/5gwkzUo9Ls+p\n+rLOQFJN1ePynKo/g4GkmnHlUHFZdCY1QLsUXLVD8Vw7MhhINdBOBVetXjzXrpwmkuaoHadNWrV4\nrtVZdCbVUTsWXNWyeM4VSPngNJE0R+06bVIqleZ8ec52ml7LO6eJpBpw2mTm2nF6rdmcJpLqrFX3\nHKqndpxeyzODgVQjrbjnUD25n1G+mDOQ1BTuZ5Qv5gwkNZWriRrH7SgkNZVf+PngdhSSmsblo8Xg\nyEBS3bh8NF8cGUiz1C6bz9WLm9oVh8FAmoDTG3PXrtXZReQ0kVSF0xu1Y3V2fliBLM2Q1bG1Y3V2\nMRgMpCqsjq0tq7Pzz5yBVIXVscVhkr826pYziIgNwCeBo1nT51NKO7LHHgDuBs4An0kpPZu1rwS+\nCXQCz6SUPjvBa5szUENYLJVvo/mIjo7ySM58xOSaUoGcBYMTKaX/Nq59BfA4sAp4N7ALeE9KKUXE\n94BPp5QGI+IZ4OGU0s4qr20wkNqcSf6Za2adQbU3vRXYllI6k1IaAg4CqyPicuDilNJgdty3gNvq\n3D9JBWUNQ23VOxh8OiJ+EBF/GhGXZm1LgVcqjjmctS0FDlW0H8raJOk81jDU1pxWE0XEc8CSyiYg\nAb8HbAJ+P5v++UPgIeCeubxfpY0bN56939fXR19fX61eWlIBjCb5+/uvG1PD4BTROQMDAwwMDEzr\n2IYUnUVEN/DXKaVrI2I9kFJKX8ge2wFsAIaBF1JKK7L2O4EPpZR+u8rrmTOQBJjkn4mm5AyyHMCo\njwE/yu4/DdwZER0RsQxYDuxOKb0GHI+I1RERwMeBp+rVP0mtoVQqsWrVKgPBHNWz6OyLEfE+4G1g\nCLgXIKW0LyKeAPYBp4F1FX/mf4qxS0t31LF/kgrAv/wbw72JJOWWdQS15ZXOJBWOdQS15/UMJBWO\ndQSNZTBQW3D/muKxjqCxDAZqeV6kppjcLLCxzBmopTnvXHyuJqodL26jtuVFaorPayE0htNEamnO\nO0vTYzBQS3PeWZoecwZqC847SxadSZKw6EySptTutSgGA0ltz1oUp4kktbl2qkVxmkiSJuAeSGUG\nA0ltzVqUMoOBpLZmLUqZOQNJoj1qUawzkCSZQJYkTc5gIKkttHtR2VQMBpJankVlUzNnIKmltVNR\n2VTMGaiwHNprriwqmx6DgXLLob1qwaKy6XGaSLnk0F61tHXrdvr717FgQTenTw+zZcsm1q69o9nd\najivgazC8drFqqW1a+/ghhs+3PJFZXNhMFAujR3al0cGDu01F6VSySAwCXMGyiX3i5Eay5yBcq0d\n9ouRGsW9iSRJ1hlIkiZnMJAkGQwkSQYDSRIGA0kSBgNJEgYDSRIGA0kSBgNJEgYDSRIGA0kSBgNJ\nEgYDSRIGA0kSBgNJEgYDSRJzDAYR8VsR8aOIeCsiVo577IGIOBgR+yNiTUX7yojYGxE/jogvV7R3\nRMS27Dl/GxFXzqVvkqTpm+vI4O+B3wC+W9kYESuA24EVwEeBTRExenWdrwH9KaWrgKsi4iNZez/w\ns5TSe4AvA1+cY99qZmBgoNldqCvPr7ha+dzA82ukOQWDlNKBlNJBYPxl1G4FtqWUzqSUhoCDwOqI\nuBy4OKU0mB33LeC2iuc8lt3/C+D6ufStlvL0gdWD51dcrXxu4Pk1Ur1yBkuBVyr+fThrWwocqmg/\nlLWNeU5K6S3gnyJicZ36J0mqMH+qAyLiOWBJZROQgN9LKf11vTrG+aMNSVK9pJTmfANeAFZW/Hs9\ncH/Fv3cA/xq4HNhf0X4n8LXKY7L784Cjk7xf8ubNmzdvM79N9L065chgBir/kn8a+POI+BLl6Z/l\nwO6UUoqI4xGxGhgEPg58peI5dwHfA/4t8PxEb5RSctQgSTU0p2AQEbcBjwC/BPyPiPhBSumjKaV9\nEfEEsA84DaxL2Z/0wKeAbwKdwDMppR1Z+xbgv0fEQeAfKY8aJEkNEOe+oyVJ7artK5AnKpyLiO6I\n+EVEfD+7bap4rDCFc+1UGBgRGyLiUMVndlPFYzM61yKIiJsi4qWs7/c3uz+zERFDEfHDiNgTEbuz\ntkUR8WxEHIiInRFxacXxVT/HvIiILRFxJCL2VrTN+Hya8ntZiwRykW/A1cB7KOcoKpPg3cDeCZ7z\nPWBVdv8Z4CPZ/d8GNmX376Bca5HX81sB7KE8VdgD/APnRoqFOb9x57oB+I9V2md8rnm/Uf5D7h+y\n39MFwA+Aa5rdr1mcx0+ARePavgB8Lrt/P/Bgdv+9E32OebkBHwTeV/ndMZvzacbvZduPDNLEhXNU\nayta4dwk59dShYEVqn2OsznXvFsNHEwpDaeUTgPbKJ9n0QTnz1BU/p49xrnP5BaqfI6N6OR0pZT+\nBjg2rnlG59Os38u2DwZT6MmmG16IiA9mba1SONeqhYGfjogfRMSfVgzHZ3OueTf+nIrU90oJeC4i\nBiPinqxtSUrpCEBK6TXgsqx9os8x7y6b4fk05feylktLc2uWhXM/Ba5MKR3L5tq/ExHvnelbz7y3\nM9dOhYGTnSuwCfj9lFKKiD8EHgLuOf9VlCMfSCm9GhEl4NmIOED586zUaqtccnk+bREMUko3zuI5\np8mGeyml70fE/wauohy9r6g49N1ZGxWP/TQi5gGXpJR+Npe+T7OvMz4/Jj6P3J1fpRmc6zeA0UA4\nm3PNu8NAZQK/SH0/K6X0avbfkYj4DuVpnyMRsSSldCSbMjmaHV7Uz2um59OU83SaaKyzf+lGxC9F\nxAXZ/V+mXDj3k2yYdzwiVmc7sX4ceCp72mjhHExRONck4wsD78xWCC3jXGFgYc8v+x9t1MeAH2X3\nZ3OueTcILM9WvXVQrst5usl9mpGIeEdELMzuXwSsobwT8tPAJ7LD7mLs7995n2NDOz09wfn/r30i\nuz/l+TTt97LZ2fdm3ygnZl4BTgKvAv8zax/9Mvk+8HfAzRXP+VeUf2kPAg9XtF8IPJG1/y+gJ6/n\nlz32AOUVDPuBNUU8v3Hn+i1gL+WVNd+hPPc8q3Mtwg24CTiQ9X19s/szi/4vyz6rPdlnsD5rXwzs\nys7tWeCfTfU55uUGPE55ivkU8DLw74FFMz2fZvxeWnQmSXKaSJJkMJAkYTCQJGEwkCRhMJAkYTCQ\nJGEwkCRhMJAkAf8fLTn9yEHWrDkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1134df7d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(pixelx,pixely)\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is what the algorithm will see."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hits from a single particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD7CAYAAAB68m/qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEv9JREFUeJzt3W2MXGd5xvH/Hexth4JDoBtHSmAmiJdNVLlgSEBCLWPA\naQApyQcUgloEBFGFqIBoRWNDpbhfCkEqL5XqD5Rt5YRo4wBtCRKFEMXbShWBhSQsJS+4glkS05hD\nKVRUlrwJdz/McTLe2N5Ze3fPzLP/n7TKmWfO7Ll2Mr529jkvE5mJJGn8ndV0AEnS6rDQJakQFrok\nFcJCl6RCWOiSVAgLXZIKsampDUeEx0tK0mnIzDjReKPv0DNzrL5uvPHGxjOUnnnc8prZvOud+VSc\ncpGkQljoklQIC30Fut1u0xFWbNwyj1teMPN6GLe80EzmWG5OZs02HJFNbVuSxlVEkKO4U1SStHos\ndEkqhIWuE6qqirm5OaqqajqKpCFZ6HqamZn9tNtT7Nx5He32FDMz+5uOJGkI7hTVcaqqot2e4siR\nA8A2YJ5WawcLCw8xOTnZdDxpw3OnqIbW6/WYmOjQL3OAbWze3KbX6zUXStJQLHQdp9PpcPRoD5iv\nR+ZZXFyg0+k0F0rSUCx0HWdycpLp6b20WjvYsmU7rdYOpqf3Ot0ijQHn0HVCVVXR6/XodDqWuTRC\nTjWHbqFL0hhxp6gkbQAWuiQVYqhCj4jdEfH9iJiPiFsjYiIizomIOyPi4Yj4WkScvWT9gxHxYERc\ntnbxJUnHLFvoEdEG3gO8PDO30f/YurcBu4C7MvOlwN3A7nr9i4GrgYuANwJ7I+KE8z2SpNUzzDv0\n/wWOAr8VEZuAFnAIuBLYV6+zD7iqXr4CuC0zH8/MHnAQuHQ1Q0uSnm7ZQs/M/wH+Gvgx/SL/ZWbe\nBWzNzMP1Oo8B59YPOR94ZOBbHKrHJElraNNyK0TEC4EPAm3gl8DnI+IPgaXHHK74GMQ9e/Y8udzt\ndsfyU0kkaS3Nzs4yOzs71LrLHoceEVcDOzPzPfXttwOvBl4HdDPzcEScBxzIzIsiYheQmXlTvf5X\ngRsz85tLvq/HoUvSCp3pcegPA6+OiN+sd26+HngAuAN4Z73OO4Av1ct3ANfUR8JcCLwI+NYZ5Jck\nDWHZKZfM/G5E3Ax8B3gCuA/4DPBs4PaIuBZYoH9kC5n5QETcTr/0F4HrfSsuSWvPU/8laYx46r8k\nbQAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEhS5JhbDQJakQFrokFcJCl6RCWOiSVAgLXZIK\nYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAW\nukZKVVXMzc1RVVXTUaSxY6FrZMzM7KfdnmLnzutot6eYmdnfdCRprERmNrPhiGxq2xo9VVXRbk9x\n5MgBYBswT6u1g4WFh5icnGw6njQyIoLMjBPd5zt0jYRer8fERId+mQNsY/PmNr1er7lQ0pix0DUS\nOp0OR4/2gPl6ZJ7FxQU6nU5zoaQxY6FrJExOTjI9vZdWawdbtmyn1drB9PRep1ukFXAOXSOlqip6\nvR6dTscyl07gVHPoFrokjRF3ikrSBmChS1Ihhir0iDg7Ij4fEQ9GxPcj4lURcU5E3BkRD0fE1yLi\n7IH1d0fEwXr9y9YuviTpmGHfoX8a+EpmXgT8LvAQsAu4KzNfCtwN7AaIiIuBq4GLgDcCeyPihPM9\nkqTVs2yhR8QW4Pcy8x8AMvPxzPwlcCWwr15tH3BVvXwFcFu9Xg84CFy62sElSccb5h36hcDPIuIf\nIuLeiPhMRDwT2JqZhwEy8zHg3Hr984FHBh5/qB6TJK2hYQp9E7Ad+NvM3A78H/3plqXHHHoMoiQ1\naNMQ6zwKPJKZ365vf5F+oR+OiK2ZeTgizgN+Wt9/CHj+wOMvqMeeZs+ePU8ud7tdut3uisJLUulm\nZ2eZnZ0dat2hTiyKiH8F3pOZP4iIG4Fn1nf9PDNviogbgHMyc1e9U/RW4FX0p1q+Drx46VlEnlgk\nSSt3qhOLhnmHDvB+4NaI2Az8EHgX8Azg9oi4Fligf2QLmflARNwOPAAsAtfb3JK09jz1X5LGiKf+\nS9IGYKFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAl\nqRAWuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0SSqEha4Npaoq5ubmqKqq6SjSqrPQtWHMzOyn3Z5i\n587raLenmJnZ33QkaVVFZjaz4YhsatvaeKqqot2e4siRA8A2YJ5WawcLCw8xOTnZdDxpaBFBZsaJ\n7vMdujaEXq/HxESHfpkDbGPz5ja9Xq+5UNIqs9C1IXQ6HY4e7QHz9cg8i4sLdDqd5kJJq8xC14Yw\nOTnJ9PReWq0dbNmynVZrB9PTe51uUVGcQ9eGUlUVvV6PTqdjmWssnWoO3UKXpDHiTlFJ2gAsdEkq\nhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1Ihhi70iDgrIu6NiDvq2+dExJ0R8XBEfC0i\nzh5Yd3dEHIyIByPisrUILkk63kreoX8AeGDg9i7grsx8KXA3sBsgIi4GrgYuAt4I7I2IE56mKkla\nPUMVekRcALwJ+OzA8JXAvnp5H3BVvXwFcFtmPp6ZPeAgcOmqpJUkndSw79A/CXwIGLya1tbMPAyQ\nmY8B59bj5wOPDKx3qB6TJK2hTcutEBFvBg5n5v0R0T3Fqiu+dOKePXueXO52u3S7p/r2krTxzM7O\nMjs7O9S6y14+NyL+Cvgj4HGgBTwb+CfglUA3Mw9HxHnAgcy8KCJ2AZmZN9WP/ypwY2Z+c8n39fK5\nkrRCZ3T53Mz8cGa+IDNfCFwD3J2Zbwe+DLyzXu0dwJfq5TuAayJiIiIuBF4EfOsMfwZJ0jKWnXI5\nhY8Bt0fEtcAC/SNbyMwHIuJ2+kfELALX+1Zcktaen1gkSWPETyySpA3AQpekQljoklQIC12SCmGh\nS6epqirm5uaoqqrpKBJgoUunZWZmP+32FDt3Xke7PcXMzP6mI0ketiitVFVVtNtTHDlyANgGzNNq\n7WBh4SEmJyebjqfCediitIp6vR4TEx36ZQ6wjc2b2/R6veZCSVjo0op1Oh2OHu0B8/XIPIuLC3Q6\nneZCSVjo0opNTk4yPb2XVmsHW7Zsp9XawfT0Xqdb1Djn0KXTVFUVvV6PTqdjmWvdnGoO3UKXpDHi\nTlFJ2gAsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAW\nuiQVwkKXpEJY6JJUCAtdkgphoUtSISx0qQFVVTE3N0dVVU1HUUEsdGmdzczsp92eYufO62i3p5iZ\n2d90JBXCzxSV1lFVVbTbUxw5cgDYBszTau1gYeEhP2haQ/EzRaUR0ev1mJjo0C9zgG1s3tym1+s1\nF0rFsNClddTpdDh6tAfM1yPzLC4u0Ol0mgulYljo0jqanJxkenovrdYOtmzZTqu1g+npvU63aFUs\nO4ceERcANwNbgV8Df5eZfxMR5wD7gTbQA67OzF/Wj9kNXAs8DnwgM+88wfd1Dl0bVlVV9Ho9Op2O\nZa4VOdUc+jCFfh5wXmbeHxHPAr4DXAm8C/jvzPx4RNwAnJOZuyLiYuBW4BLgAuAu4MVL29tCl6SV\nO6Odopn5WGbeXy//CniQflFfCeyrV9sHXFUvXwHclpmPZ2YPOAhcekY/gSRpWSuaQ4+IDvAy4B5g\na2Yehn7pA+fWq50PPDLwsEP1mCRpDQ1d6PV0yxfoz4n/Clg6X+L8iSQ1aNMwK0XEJvplfktmfqke\nPhwRWzPzcD3P/tN6/BDw/IGHX1CPPc2ePXueXO52u3S73RWFl0rkDlMNmp2dZXZ2dqh1hzpTNCJu\nBn6WmX86MHYT8PPMvOkkO0VfRX+q5eu4U1QayszMft797uuZmOgfrz49vZe3ve2tTcfSCDnTo1xe\nA/wb8D360yoJfBj4FnA7/XfjC/QPW/xF/ZjdwLuBRTxsURqKlwXQME5V6MtOuWTmvwPPOMndbzjJ\nYz4KfHTohJKevCzAkSNPvyyAha5heKaoNCK8LIDOlIUujQgvC6Az5eVzpRHjUS46lTPaKbpWLHRJ\nWjmvhy5JG4CFLkmFsNAlqRAWujQmqqpibm6OqqqajqIRZaFLY2BmZj/t9hQ7d15Huz3FzMz+piNp\nBHmUizTivCSABnmUizTGjl0SoF/mMHhJAGmQhS6NOC8JoGFZ6NKI85IAGpZz6NKY8JIAAk/9l6Ri\nuFNUkjYAC12SCmGhS1IhLHRJKoSFLkmFsNAlqRAWulQIr8YoC10qgFdjFHhikTT2vBrjxuKJRVLB\nvBqjjrHQpTHn1Rh1jIUujTmvxqhjnEOXCuHVGDcGr7YoSYVwp6gkbQAWulQwTzbaWCx0qVCebLTx\nOIcuFciTjcrlHLq0wXiy0cZkoUsF8mSjjclClwrkyUYb05rNoUfE5cCn6P/SmM7Mm5bc7xy6tMY8\n2ag8635iUUScBfwAeD3wE2AOuCYzHxpYx0KXpBVqYqfopcDBzFzIzEXgNuDKNdqWJIm1K/TzgUcG\nbj9aj0mS1simJje+Z8+eJ5e73S7dbrexLJI0imZnZ5mdnR1q3bWaQ381sCczL69v7wJycMeoc+iS\ntHJNzKHPAS+KiHZETADXAHes0bYkSazRlEtmPhERfwLcyVOHLT64FtuSJPV5LRdJGiNey0WSNgAL\nXZIKYaFLUiEsdEkqhIUuSYWw0CWpEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmFsNAl\nqRAW+goM+zFQo2TcMo9bXjDzehi3vNBMZgt9BXxRrb1xywtmXg/jlhcsdEnSGbDQJakQjX4EXSMb\nlqQxd7KPoGus0CVJq8spF0kqhIUuSYVYk0KPiLdExH9ExBMRsX3Jfbsj4mBEPBgRlw2Mb4+I+Yj4\nQUR8amB8IiJuqx/zjYh4wVpkXpLxkoj4VkTcV//3laebfz1FxPvqXN+LiI+NQ+Y6x59FxK8j4rkD\nYyOXOSI+Xue5PyK+GBFbRjnviUTE5RHxUJ3nhqbzAETEBRFxd0R8v37tvr8ePyci7oyIhyPiaxFx\n9sBjTvh8N5D9rIi4NyLuGInMmbnqX8BLgRcDdwPbB8YvAu4DNgEd4D95ah7/m8Al9fJXgD+ol98L\n7K2X3wrcthaZl+Q/AFxWL78ROFAvX7zS/Ov1BXSBO4FN9e3fPt3nfJ1zXwB8FfgR8NxRzgy8ATir\nXv4Y8NFRf10syX9Wna0NbAbuB6aayjOQ6zzgZfXys4CHgSngJuDP6/EbgI8t93w3kP2DwOeAO+rb\njWZek3fomflwZh4Elu6JvZJ+IT+emT3gIHBpRJwHPDsz5+r1bgauGnjMvnr5C8Dr1yLzEv8FHPvN\n+hzgUL18BSvPv17eS//F8zhAZv6sHj+d53w9fRL40JKxkcycmXdl5q/rm/fQ/2UEo/26GHQpcDAz\nFzJzEbiN/nPdqMx8LDPvr5d/BTxI/7kd/Le/j6eeuxM+3+samv5fFsCbgM8ODDeaeb3n0M8HHhm4\nfageOx94dGD80XrsuMdk5hPALwb/NF8ju4BPRMSPgY8Du5dmqQ2Tf728BPj9iLgnIg5ExCvq8ZHN\nHBFXAI9k5veW3DWymQdcS/8dN4xHXnh6zqbzPE1EdICX0f+FuTUzD0O/9IFz69VO9nyvt2NvRgYP\nFWw086bTfWBEfB3YOjhE/wf7SGZ++UyDnWrTq/JNTp7/L4D3Ae/LzH+OiLcAfw/sXI3tnollMm8C\nzsnMV0fEJcDngReuf8rjLZP5w4zA8zpomNd1RHwEWMzMmQYiFisinkX/r/APZOavTnCuysgcYx0R\nbwYOZ+b9EdE9xarrmvm0Cz0zT+cf4iHg+QO3L6jHTjY++JifRMQzgC2Z+fPT2PZxTpU/Ij537P7M\n/EJEHPuT6nTyr5plMl8H/GO93ly9Q/p5dY7BHckjkTkifof+XOJ3IyLq7d8bEZc2mXm513VEvJP+\nn9mvGxhu9HWxAid7XhsXEZvol/ktmfmlevhwRGzNzMP19NVP6/FReF5fA1wREW8CWsCzI+IW4LFG\nM6/xDoMDwCsGbh/bMTABXMjxO4/uoT+nFPT/lL28Hr+ep3aKXsP67BT9DvDaevn1wNzp5l+vL+CP\ngb+sl18CLIx65iX5f0T/L4yRzQxcDnwfeN6S8ZHMe4L8z+CpnaIT9HeKXtRUniXZbgY+sWTsJuCG\nevlEOxif9nw3lP21PLVT9ONNZl6rH/Aq+vNFR+jvYPyXgft21z/Mg9RHktTjrwC+R39nwacHxn8D\nuL0evwforMP/oFfSPzrhPuAbwMtPN/86vqg2A7fUGb5N/QtplDMvyf9D6qNcRjVzvc0F4N76a+8o\n5z3Jz3A5/aNIDgK7ms5TZ3oN8ET9C+a++rm9HHgucFed907gOcs93w3lHyz0RjN76r8kFcIzRSWp\nEBa6JBXCQpekQljoklQIC12SCmGhS1IhLHRJKoSFLkmF+H+hQewv7SQnXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e2f3f10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "iparticle = 2\n",
    "hits_particle = (particles == iparticle)\n",
    "plt.scatter(pixelx[hits_particle],pixely[hits_particle])\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, ideally, our algorithm wants to assign specific ids to all hits in a way that corresponds to the particle. The perfect algorithm will split the event into various tracks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGetJREFUeJzt3X+M3Gd94PH3Z9c73sGJiSmDI4Wy6xPgmMKpAdkXjupu\nJewQ+COGqs1l+wfJsZxQAwjCSZBc/0jUq1rgxMXhUNDpOi0hVzbkqh5xT7kkWLAnISVkCQlJmuCa\nH7MkLjHTNnESs+sd28/9Md+1x+v94fXOj+/MvF/SyN999js7z5PZzGef5/P8iJQSkqT+NtDpCkiS\nOs9gIEkyGEiSDAaSJAwGkiQMBpIkmhQMIqIcEUci4smGsi0R8VBEHIyIByPitQ3fuyUiDkXEsxFx\nVUP5OyPiyYj4+4jY14y6SZJW16yewV8C71tUdjNwIKW0HfgOcAtARLwNuBbYAbwfuDMiInvOV4GJ\nlNJbgbdGxOKfKUlqgaYEg5TS94AXFxXvBe7Kru8CPphdXwPck1I6kVKqAIeAXRFxKXBxSmk6u+/r\nDc+RJLVQK3MGb0gpHQFIKb0AvCErvwx4ruG+w1nZZcDzDeXPZ2WSpBZrZwLZfS8kKac2tPBnH4mI\nrSmlI9kQ0K+y8sPAbzbc98asbLnyc0SEgUWSLkBKKZYqb2bPILLHgv3ADdn19cB9DeXXRUQhIrYB\nbwYezYaSjkbEriyh/OGG55wjpdS2x6233trW12v3w/Z176OX22b7mv9YSVN6BhHxDWAM+I2I+AVw\nK/B54H9FxEeAGeoziEgpPRMR9wLPADXgxnSmlh8HvgYMA/enlB5oRv0kSStrSjBIKf3BMt/avcz9\nfwb82RLljwHvaEadJEnnzxXI52FsbKzTVWgp29e9erltYPvaKVYbR8qjiEjdWG9J6qSIILUhgSxJ\n6lIGA0mSwUCSZDCQJGEwkCRhMJAkYTCQJGEwkCRhMFCPmZ+v8vLL08zPVztdFamrtHILa6mtjhyZ\n5ODBCSIKpDTP9u1ltm4d73S1pK7gdhTqCfPzVR55ZIRTp2ZPlw0MFLnyyhkKhVIHayblh9tRqOfN\nzVWIKJxVFjHE3FylMxWSuozBQD1heHiUlObPKkupxvDwaGcqJHUZg4F6QqFQYvv2MgMDRQYHNzMw\nUGT79rJDRNJ5MmegnjI/X2VursLw8KiBQFpkpZyBwUBdoVqFSgVGR6HkZ7x0QUwgq6tNTsLICOzZ\nU/93crLTNZJ6jz0D5Vq1Wg8As2dmjFIswsyMPQRprewZqGtVKlA4e8YoQ0P1cknNYzBQro2OwvzZ\nM0ap1erlkprHYKBcK5WgXK4PDW3eXP+3XHaISGo2cwbqCs4mktbPqaWSJBPIkqSVGQwkSQYDSZLB\nQJKEwUCShMFAkoTBQJKEwUA6bX6+yssvTzM/X+10VaS229DpCkh5cOTIJAcPThBRIKV5tm8vs3Xr\neKerJbWNK5DV9+bnqzzyyAinTp3ZJ3tgoMiVV854Wpp6iiuQpRXMzVWIOHuf7Igh5uYqnamQ1AEG\nA/W94eFRUjp7n+yUagwPj3amQlIHGAzUN6pVmJ6u/9uoUCixfXuZgYEig4ObGRgosn172SEi9RVz\nBuoLk5MwMVE/NW1+vn4mwvii/PD8fJW5uQrDw6MGAvUkt7BWX/McZanOBLL6mucoS6szGKjneY6y\ntDqDgXqe5yhLqzNnoL7hOcrqdyaQJUkmkCVJKzMYSJJaHwwiohIRP4qIxyPi0axsS0Q8FBEHI+LB\niHhtw/23RMShiHg2Iq5qdf0kSe3pGZwCxlJKV6SUdmVlNwMHUkrbge8AtwBExNuAa4EdwPuBOyNi\nyfEtSVLztCMYxBKvsxe4K7u+C/hgdn0NcE9K6URKqQIcAnYhSWqpdgSDBHw7IqYj4qNZ2daU0hGA\nlNILwBuy8suA5xqeezgrk3LHk9HUS9px0tl7Ukq/jIgS8FBEHKQeIBqteZ7obbfddvp6bGyMsbGx\n9dRRWhNPRlM3mJqaYmpq6rzubes6g4i4FXgV+Cj1PMKRiLgU+G5KaUdE3AyklNIXsvsfAG5NKX1/\n0c9xnYE6xpPR1K06ts4gIl4TERdl15uAq4CngP3ADdlt1wP3Zdf7gesiohAR24A3A4+2so7SWnky\nmnpRq4eJtgL/OyJS9lp/lVJ6KCJ+ANwbER8BZqjPICKl9ExE3As8A9SAG+0CqNMWb2PhyWjqRW5H\nIa1guUNxzuQMhkipZs5AXcG9iaQLsNqhOJ6Mpm6zUjBox2wiqSstHIrTGAwWDsUplepnJxsE1Cvc\nm0hahofiqJ8YDKRleCiO+ok5A2kVHoqjXmECWZLk4TaSpJUZDCRJBgNJksFAaim3uVa3cNGZ1CJu\nc61u4mwiqQXc5lp55HYUUpssrEnYunVhm+szwWBhm2uDgfLInIHUJJOT9Y3t9uyBXbtGqdXc5lrd\nw2AgNUG1Wt/qenYWjh6FI0dKfPGLZSKKDA5uZmCgyPbtZXsFyi2HiaQmWGqH04cfHqdY3M3ll7vN\ntfLPYCA1wXI7nG7bVmLzZoOA8s9hIqkJ3OFU3c6ppVITucOp8sxdSyVJ7loqSVqZwUCSZDCQOsmN\n7JQXTi2VOsSN7JQnJpClNqtW4ec/rzI7O0JKbmSn9jGBLOXEwv5FN95Y4dVXC2d9b2EjO6kTDAZS\nmzTuX/STn4wyOOhGdsoPg4HUJgv7FwEcPVrfyO748SLgRnbqPHMGqnPpbMtVq/UhosbN7LZurfLo\noxUuvXQUgLk5N7VT65gz0MoaN+IfGal/raZbav+i228vUSzu5Ac/OMDDD4/wox/t4ZFHRjhyxPdA\n7WXPoN8t9edqsQgzM/YQWqSxE3bgANx0U5W77hph40ZnFqm17BloeY0D2QuGhurlaolSCXburF9P\nTMAll1So1ZxZpM4yGPS75TbiHx3tRG36ykIcfuGFUTZscGaROstg0O/ciL9jFuLwwsyiubkix45t\nJqI+swhwqwq1jTkD1bVrNpGzls4yOVkfKhoagmKxyle+UuE97xnl5z8/QK02wcCAW1WoeTzPQPmw\n8MlXKNT/JC6XYdwPOBPKahcTyOq8xuW3R4/W/52YqJf3ORPKygODgdrDWUurMqGsTjIYqD1aPWup\nWoXp6a7uaayWUHaISK1kMFB7tHLWUo+soG78T/TYY+PccMMMs7MHePe7Z0weq+VMIKu9mj2bqAdX\nUC/3n2h+vureRVqXlRLInnSm9iqVmvshvTDQ3hgMFnIRa32dnEx7Xeo/0aFDkxw+XJ9qCk41VfM5\nTKTu1qxcRI6Hmr75zSo//ekEMMupU0c5dWqWgwcnXIympjIYqLs1IxfRjGmvLUpgV6vwp39a4cSJ\nRTOxcKqpmstgoO43Pl7PERw4UP93rQvZ1jvttYW9ikoFXnzx3Kmmp04dZ3Dwoqa9jmQCWVpPEnq1\n564zD7Hw46+8cpLPfnYCSGzcOAcUGRiAyy83d6Dz5wpkaSXrGWpaqVexUo/hPIeVFqr2yCPjfOYz\njzEwkIiAiFlSmuXv/s7cgZojdz2DiLga2Ec9UJVTSl9Y4h57Bmq+C/krfrmewWOPwbvetXSP4cCB\n1fdoWlSXahWeeGKa+fk9bNp09PRtx45t5h3vOMCb3rTzAhutftI1PYOIGAC+ArwP+C1gPCIu72yt\n1DcWNglay3DOcr2KV19dusfw+OOrJ6uX6FGUSvAbpy7idT+dY+ilM7du2FDjhRdG19VsCfK3zmAX\ncCilNAMQEfcAe4Efd7RW0krGx2H37rN7FdXq0lNeYeV1EY0zmxbumZiAl1/miptu4lhtgOIgPPHJ\nYarvDfbtK3P33S5A0/rlqmcAXAY81/D181mZlG+LexXL9RiuuGLldRFL5SAGB+FTnyJmZ7noxCyD\nx+G3vpT4jx9+jOuvH+/WhdbKmbz1DM7bbbfddvp6bGyMsbGxjtVFWtJSPQaoB4WFE21qtbOT1cst\noisU4Pjx00UbXrOR7/3PV9lyVf3r6rEqj//0cXgJrthxBSUjhICpqSmmpqbO695cJZAj4krgtpTS\n1dnXNwNpcRLZBLK63krJ6sbjz2o1uP12uOmms4aWTm7cyIuPP87rd+xg8qlJrv/P11P7Vg0GoBAF\nvvYXX2Pcg4O0SNecdBYRg8BB4L3AL4FHgfGU0rOL7jMYqLctDhZZgJgHTszO8olikXuA2//77Xz6\n2U8z91/m4MSZpxeLRWZmZuwh6CxdEwzg9NTSOzgztfTzS9xjMFDf+cdnn+VDV1zBj48f5x+zso3b\nNjL4/kF+Xf41nBlFYtPFm/ibv/4btmzZwujoqEFBQJftWppSegDY3ul6SHnz81df5anhYY425A6G\nfj3E/MXzcPLse4/PHeeavddQKBQ4UTtBuVx22EgryttsIklLqFarvPjii8wvSi6ffPkkX77mywx9\naAiGgI0wODTIiVMnOD53nFdefoXZ2VkmJiaodvEpcGo9g4GUc5OTk4yMjHDttddy4sQJCoUCmzdv\nplgsUi6X+di//hiH//IwD/7gQb65/5sM/sHgOX3+DUMbqHjetFaQu5zB+TBnoH5RrVYZGRlhtmEm\n0fDwMPfddx9XXHHuFNLpw9O896vv5ZUvvHJWQnnj8Eae+8Vzp++vVqtUKhXzCX2ma7ajkHS2SqVC\nYdEitEKhwJYtW5b8EB+9ZJQTxRP1dfsbgI3AENzx1TtO37/Q09izZw8jIyNM5uggH3WOPQMpx5bq\nGaw2bXTy6Ukm7ptgcHaQ2j/VuGP8Dj72bz52wT9PvaOrZhNJOqNUKlEul5mYmGBoaIharUa5XF7x\ng3v87ePs3rabyksVRi8ZpbQpGxo6VuX+R+9nw9AGOGtrpCEqlYrBoM/ZM5ByZqnx/PWO8U8+NcnE\n/gk2zG7glS++ArUz31voGQDmEXpcVy06Ox8GA/WqyclJJiYmKBQKzM/PN2V9QPVYlZF9I8yeyLoD\nTwH74eLixafXIABNf13lj8FA6gKtGs+fPjzNnrv3cPT4mUNxLqpdxFfe/RU+sOsDAOYR+oSziaQu\nsNTMoYXx/PUYvWSU+ZOLFqsVT/KBsQ/AaziTR2jy66q7GAyknBgdHT1nhXGtVmN04ayDC1TaVKK8\nt0xxQ5HNGzdT3FCkvLfMgZ8dYGTfCJ98+JO8MvvKkq9brVaZnp529XIfcJhI6rDG5PCBAwfOmTnU\nrLH76rHq6RlGgHmEPmTOQMqppRLGu3fvbvmsHvMI/cmcgZRD1WqViYkJZmdnOXr06OkN5QB27tzZ\n0g/dlfIIpVKpZfkL5ZfBQOqQTn7gLpdHWFig1qr8hfLLYSKpQ/KwNURjHmEhECxYGMJqRf5CnWHO\nQMqRdiWMm8HdTXuLexNJObFUwnhmZia3H7ilUil3dVJr2DOQ2iQPw0Lqb84mknKgH2bouEitexkM\npDbp9Rk6HprT3RwmktqoV2foOATWHUwgSzkxPj7elhXG7bYwBNYYDDw0p7sYDKQ268UZOr0+BNYP\nzBlIWreF4zmLxSKbN2+mWCyuejyn8sWcgaQLstTqZRep5ZsrkCU11cKZyoXBAvMn5ynvLTP+9u5P\nhPc6g4GkpjnnTGWguKHIzKdnztnfSPniojOpDfplwVXlpQqFwUWL5waHqLxU6UyF1BQGA6kJ+mnB\n1VJnIdRO1k6foKbu5DCRtE79uOBq8ulJJu6bYGhwiNrJmjmDLuGiM6mF+nHB1fjbx9m9bfeyZyGs\nhTOQ8sFhImmd+nXBVWlTiZ2X7VxXIOin4bW8c5hIaoJe3XOolfpxeK3THCaSWqxX9xxqpX4cXssz\ng4HUJL2451Ar9evwWl6ZM5DUEe5nlC/mDCR1lLOJ2sftKCR11FKb2qn93I5CUsdMPjXJyL4R9ty9\nh5F9I0w+7fTRPLJnIKll3NQuX+wZSBeoXzafaxU3teseBgNpGa6OXT83teseDhNJS3B1bPO4qV1+\nuAJZWiNXxzZPMze1U+sYDKQluDq2uUqbSgaBnDNnIC3B1bHdwyR/c7QsZxARtwL/AfhVVvSfUkoP\nZN+7BfgIcAL4VErpoaz8ncDXgGHg/pTSp5f52eYM1Baujs23hd1iC4UC8/Pz7ha7io6sQM6CwSsp\npf+6qHwH8A1gJ/BG4ADwlpRSiojvA59IKU1HxP3AHSmlB5f42QYDqc+Z5F+7Tq4zWOpF9wL3pJRO\npJQqwCFgV0RcClycUprO7vs68MEW109Sl1pI8jdaSPJr7VodDD4REU9ExJ9HxGuzssuA5xruOZyV\nXQY831D+fFYmSecwyd9c65pNFBHfBrY2FgEJ+CPgTuCPs+GfPwG+BHx0Pa/X6Lbbbjt9PTY2xtjY\nWLN+tKQusJDkX3zCnENEZ0xNTTE1NXVe97Zl0VlEjAB/m1L6lxFxM5BSSl/IvvcAcCswA3w3pbQj\nK78O+LcppT9c4ueZM5AEmORfi47kDLIcwILfBZ7OrvcD10VEISK2AW8GHk0pvQAcjYhdERHAh4H7\nWlU/Sb2hVCqxc+dOA8E6tXLR2Rcj4reBU0AF+BhASumZiLgXeAaoATc2/Jn/cc6eWvpAC+snqQt4\nFkJ7uDeRpNyafGqSif0TFAYLzJ+cd1+jdfKkM0ldx7MQms/zDCR1Hc9CaC+DgfqC+9d0H89CaC+D\ngXqeh9R0p9KmEuW9ZYobimzeuJnihiLlvWWHiFrEnIF6mvvXdD9nEzWPh9uob3lITffzLIT2cJhI\nPc39a6TzYzBQT/OQGun8mDNQX3D/GslFZ5IkXHQmSavq97UoBgNJfc+1KA4TSepz/bQWxWEiSVqG\nZynXGQwk9TXXotQZDCT1Ndei1JkzkCT6Yy2K6wwkSSaQJUkrMxhI6gvVY1WmD09TPdafi8pWYzCQ\n1PMmn5pkZN8Ie+7ew8i+ESaf7r9FZasxZyCpp1WPVRnZN8LsiYZFZRuKzHx6pu/OSTBnoK7V7/vF\naP0qL1UoDC5aVDY4ROWlSmcqlFMGA+WW+8WoGUYvGWX+5KJFZSdrjF4y2pkK5ZTDRMqlftovRq03\n+fQkE/dNMDQ4RO1kjfLeMuNvH+90tdrOM5DVdTy7WM00/vZxdm/bTeWlCqOXjPZdruB8GAyUS+4X\no2YrbSoZBFZgzkC55H4xUnuZM1Cu9cN+MVK7uDeRJMl1BpKklRkMJEkGA0mSwUCShMFAkoTBQJKE\nwUCShMFAkoTBQJKEwUCShMFAkoTBQJKEwUCShMFAkoTBQJKEwUCSxDqDQUT8XkQ8HREnI+Kdi753\nS0QciohnI+KqhvJ3RsSTEfH3EbGvobwQEfdkz3k4It60nrpJks7fensGTwEfAv5fY2FE7ACuBXYA\n7wfujIiF03W+CkyklN4KvDUi3peVTwD/nFJ6C7AP+OI669Y0U1NTna5CS9m+7tXLbQPb107rCgYp\npYMppUPA4mPU9gL3pJROpJQqwCFgV0RcClycUprO7vs68MGG59yVXf818N711K2Z8vSGtYLt6169\n3Dawfe3UqpzBZcBzDV8fzsouA55vKH8+KzvrOSmlk8BLEfG6FtVPktRgw2o3RMS3ga2NRUAC/iil\n9Letqhjn9jYkSa2SUlr3A/gu8M6Gr28GPtfw9QPAvwIuBZ5tKL8O+GrjPdn1IPCrFV4v+fDhw4eP\ntT+W+1xdtWewBo1/ye8H/ioibqc+/PNm4NGUUoqIoxGxC5gGPgx8ueE51wPfB34f+M5yL5RSstcg\nSU20rmAQER8E/hvweuD/RMQTKaX3p5SeiYh7gWeAGnBjyv6kBz4OfA0YBu5PKT2QlZeBuyPiEPBP\n1HsNkqQ2iDOf0ZKkftX3K5CXWzgXESMR8euI+GH2uLPhe12zcK6fFgZGxK0R8XzDe3Z1w/fW1NZu\nEBFXR8SPs7p/rtP1uRARUYmIH0XE4xHxaFa2JSIeioiDEfFgRLy24f4l38e8iIhyRByJiCcbytbc\nno78XjYjgdzND2A78BbqOYrGJPgI8OQyz/k+sDO7vh94X3b9h8Cd2fW/o77WIq/t2wE8Tn2ocBT4\nCWd6il3TvkVtvRX4zBLla25r3h/U/5D7SfZ7OgQ8AVze6XpdQDt+BmxZVPYF4LPZ9eeAz2fXb1vu\nfczLA/gd4LcbPzsupD2d+L3s+55BWn7hHEuVddvCuRXa11MLAxss9T5eSFvzbhdwKKU0k1KqAfdQ\nb2e3Cc4doWj8PbuLM+/JNSzxPrajkucrpfQ94MVFxWtqT6d+L/s+GKxiNBtu+G5E/E5W1isL53p1\nYeAnIuKJiPjzhu74hbQ17xa3qZvq3igB346I6Yj4aFa2NaV0BCCl9ALwhqx8ufcx796wxvZ05Pey\nmVNLc+sCF879A/CmlNKL2Vj7tyLibWt96bXXdu36aWHgSm0F7gT+OKWUIuJPgC8BHz33pyhH3pNS\n+mVElICHIuIg9fezUa/Ncslle/oiGKSU9lzAc2pk3b2U0g8j4qfAW6lH799suPWNWRkN3/uHiBgE\nNqeU/nk9dT/Puq65fSzfjty1r9Ea2vo/gIVAeCFtzbvDQGMCv5vqflpK6ZfZv9WI+Bb1YZ8jEbE1\npXQkGzL5VXZ7t75fa21PR9rpMNHZTv+lGxGvj4iB7PpfUF8497Osm3c0InZlO7F+GLgve9rCwjlY\nZeFchyxeGHhdNkNoG2cWBnZt+7L/0Rb8LvB0dn0hbc27aeDN2ay3AvV1Ofs7XKc1iYjXRMRF2fUm\n4CrqOyHvB27Ibrues3//znkf21rp8xOc+//aDdn1qu3p2O9lp7PvnX5QT8w8B8wCvwT+b1a+8GHy\nQ+AHwAcanvMu6r+0h4A7Gso3Avdm5Y8Ao3ltX/a9W6jPYHgWuKob27eorV8HnqQ+s+Zb1MeeL6it\n3fAArgYOZnW/udP1uYD6b8veq8ez9+DmrPx1wIGsbQ8Bl6z2PublAXyD+hDzceAXwL8Htqy1PZ34\nvXTRmSTJYSJJksFAkoTBQJKEwUCShMFAkoTBQJKEwUCShMFAkgT8f39S8+VYg9OqAAAAAElFTkSu\nQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1134cced0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cmap = ['k','b','y','g','r','k','b']\n",
    "for iparticle in np.unique(particles):\n",
    "    hits_particle = (particles == iparticle)\n",
    "    plt.scatter(pixelx[hits_particle],pixely[hits_particle],color=cmap[int(iparticle) % 7])\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matching the hits together"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Choose the tracking class that you want to use\n",
    "(in the submission, it has to be derived from BaseEstimator and named Clusterer in clusterer.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from hough1 import Clusterer\n",
    "tracker = Clusterer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Training the algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tracker.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Running the algorithm to reconstruct the tracks in data\n",
    "This will take some time to go through all events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_true = y\n",
    "y_predicted = tracker.predict_single_event(X_event)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare how well the tracks corresponds to the particles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>particle</th>\n",
       "      <th>track</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>47</th>\n",
       "      <td>0</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>49</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>45</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>1</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>37</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>34</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>2</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>46</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>39</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>2</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>48</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>33</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>3</td>\n",
       "      <td>7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>26</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>30</th>\n",
       "      <td>3</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>3</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>4</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4</td>\n",
       "      <td>9.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>51</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>31</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>36</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>38</th>\n",
       "      <td>4</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>35</th>\n",
       "      <td>5</td>\n",
       "      <td>6.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>28</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>29</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>5</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>5</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    particle  track\n",
       "0          0    3.0\n",
       "50         0    3.0\n",
       "47         0    9.0\n",
       "43         0    3.0\n",
       "7          0    3.0\n",
       "10         0   10.0\n",
       "11         0    3.0\n",
       "32         0    3.0\n",
       "23         0    3.0\n",
       "49         1    1.0\n",
       "45         1    1.0\n",
       "44         1    8.0\n",
       "40         1    1.0\n",
       "37         1    1.0\n",
       "34         1    1.0\n",
       "17         1    1.0\n",
       "52         1    1.0\n",
       "3          1    1.0\n",
       "15         2    0.0\n",
       "18         2    7.0\n",
       "46         2    0.0\n",
       "22         2    0.0\n",
       "14         2    0.0\n",
       "24         2    0.0\n",
       "8          2    0.0\n",
       "39         2    0.0\n",
       "16         2    0.0\n",
       "48         3    5.0\n",
       "4          3    5.0\n",
       "33         3    5.0\n",
       "12         3    7.0\n",
       "26         3    5.0\n",
       "21         3    5.0\n",
       "30         3    5.0\n",
       "27         3    6.0\n",
       "19         4    2.0\n",
       "20         4    8.0\n",
       "6          4    9.0\n",
       "42         4    2.0\n",
       "51         4    2.0\n",
       "31         4    2.0\n",
       "9          4    2.0\n",
       "36         4    2.0\n",
       "38         4    2.0\n",
       "35         5    6.0\n",
       "5          5    4.0\n",
       "25         5    4.0\n",
       "28         5    4.0\n",
       "29         5    4.0\n",
       "13         5   10.0\n",
       "2          5    4.0\n",
       "1          5    4.0\n",
       "41         5    4.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_result = pd.DataFrame({'particle':y_true[event_indices,1],'track':y_predicted}).sort_values(by='particle')\n",
    "df_result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tracks = y_predicted"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0 50 47 43  7 10 11 32 23 49 45 44 40 37 34 17 52  3 15 18 46 22 14 24  8\n",
      " 39 16 48  4 33 12 26 21 30 27 19 20  6 42 51 31  9 36 38 35  5 25 28 29 13\n",
      "  2  1 41]\n",
      "[0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 4\n",
      " 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5]\n",
      "[  3.   3.   9.   3.   3.  10.   3.   3.   3.   1.   1.   8.   1.   1.   1.\n",
      "   1.   1.   1.   0.   7.   0.   0.   0.   0.   0.   0.   0.   5.   5.   5.\n",
      "   7.   5.   5.   5.   6.   2.   8.   9.   2.   2.   2.   2.   2.   2.   6.\n",
      "   4.   4.   4.   4.  10.   4.   4.   4.]\n",
      "0\n",
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGtVJREFUeJzt3X1sXNd55/Hvw5cRJ5JoKc6Y3soJx4ljxt5kETuxKmyL\nLAGLjh2gkTe76zUL1PaaLtK8dJMu0MTe/CFtUWyboN0o3cBB0Mw2jtGM4bTYWC28iiLEXCBoaTON\n/FbbEyXBsJbj0LOpLDk2qTvkPPvHvSQvqeH7vNyZ+X2AAUeHd8hzMRQfnvOc8xxzd0REpLN1NbsD\nIiLSfAoGIiKiYCAiIgoGIiKCgoGIiKBgICIi1CgYmFnOzKbN7OlY214zO2FmBTP7jpldEvvcfWZ2\n2syeN7ObYu3Xm9nTZvYjMztai76JiMj6ajUy+Avggyva7gVOuvsQ8D3gPgAzuxa4DbgGuAW438ws\nes1XgDF3vxq42sxWfk0REamDmgQDd/8+cHZF8yHggej5A8Ct0fMPAw+5+5y7F4HTwH4zuxzY7e6T\n0XXfiL1GRETqqJ45g8vcfRrA3X8OXBa17wNejF33UtS2DzgTaz8TtYmISJ01MoGsuhciIgnVU8ev\nPW1mA+4+HU0BvRK1vwS8NXbdFVHbau0XMTMFFhGRLXB3q9Zey5GBRY8Fx4C7oud3Ao/E2m83s5SZ\nXQlcBTwRTSWdM7P9UUL5jthrLuLuDXscPny4od+v0Q/dX+s+2vnedH+1f6ylJiMDM/smMAxcamb/\nBBwG/hj4lpndDUwRriDC3Z8zs4eB54Ay8HFf6uUngK8DfcCj7n68Fv0TEZG11SQYuPtvrvKpg6tc\n/0fAH1Vp/wfgPbXok4iIbJx2IG/A8PBws7tQV7q/1tXO9wa6v0ay9eaRksjMvBX7LSLSTGaGNyCB\nLCIiLUrBQEREFAxERETBQEREUDAQEREUDEREBAUDERFBwUBERFAwkDZTKsHkZPhRRDZOwUDaRj4P\ng4MwMhJ+zOeb3SOR1qFyFNIWSqUwAMzMLLWl0zA1BZlM8/olkiQqRyFtr1iEVGp5W29v2C4i61Mw\nkLaQzUIQLG8rl8N2EVmfgoG0hUwGcrlwaqi/P/yYy2mKSGSjlDOQtlIqhVND2awCgchKa+UMFAyk\nJQSlgNniLH3ZPlKZ1PovEJGLrBUManLspUg9TeenKYwVsJThgTOUG2JgdKDZ3RJpKxoZSKIFpYCJ\nwQkqM5XFtq50FwemDmiEILJJWloqLWu2OIullv/sWq8xW5xtUo9E2pOCgSRaX7YPD5aPAr3s9GX7\nmtQjkfakYCCJlsqkGMoN0ZXuoru/m650F0O5IU0RidSYcgbSErSaSGT7tLRURESUQBYRkbUpGIiI\niIKBiIgoGIiICAoGIiKCgoGIiKBgICIiKBiILCqVYHIy/CjSaRQMRIB8HgYHYWQk/JjPN7tHIo2l\nHcjS8UqlMADMzCy1pdMwNaXT0qS9aAeyyBqKRUitKHfU2xu2i3QKBQPpeNksBMHytnI5bBfpFAoG\n0jGCUsD5yfMEpeW/+TMZyOXCqaH+/vBjLqcpIuksyhlIR9jIOcqlUjg1lM0qEEh7Uglr6Wg6R1kk\npASydDSdoyyyPgUDaXs6R1lkfQoG0vZ0jrLI+pQzkI6hc5Sl0ymBLCIiSiCLiMjaFAxERKT+wcDM\nimb2lJmdMrMnora9ZnbCzApm9h0zuyR2/X1mdtrMnjezm+rdPxERaczIoAIMu/t17r4/arsXOOnu\nQ8D3gPsAzOxa4DbgGuAW4H4zqzq/JSIitdOIYGBVvs8h4IHo+QPArdHzDwMPufucuxeB08B+RESk\nrhoRDBz4rplNmtk9UduAu08DuPvPgcui9n3Ai7HXvhS1iSSOTkaTdtLTgO/xa+7+spllgBNmViAM\nEHGbXid65MiRxefDw8MMDw9vp48im5LPw9hYeA5CEIRVTkdHm90rkeXGx8cZHx/f0LUN3WdgZoeB\nXwL3EOYRps3scuAxd7/GzO4F3N0/H11/HDjs7o+v+DraZyBNo5PRpFU1bZ+Bmb3JzHZFz3cCNwHP\nAMeAu6LL7gQeiZ4fA243s5SZXQlcBTxRzz6KbJZORpN2VO9pogHgf5uZR9/rL939hJn9AHjYzO4G\npghXEOHuz5nZw8BzQBn4uIYA0mwry1joZDRpRypHIbKG1Q7FWcgZ9PaGgUA5A2kFqk0ksgXrHYqj\nk9Gk1awVDBqxmkikJS0eihNLFC8cipPKpMhkFASkfag2kcgqdCiOdBIFA5FV6FAc6STKGYisQ4fi\nSLtQAllERHS4jYiIrE3BQEREFAxERETBQKSuVOZaWoWCgUid5PNhddORkfBjPt/sHomsTquJROpA\nZa4libSaSKRBglLA+cnznD5VVplraSmqTSRSI/EKp2cv9HChcoDwCPCQylxLkmlkIFIDQSmgMFag\nMlNh/tw8/bMX+H1/gXTa6e8Pp4hyOU0RSXJpZCBSA9UqnI6k/x93feuX/GLvbpW5lsRTMBCpgdUq\nnF5x3Q7eriAgLUDTRCI1oAqn0uq0tFSkhlThVJJMVUtFRET7DEREZG0KBiIiomAg0kwqZCdJoWAg\n0iQqZCdJogSySIMFpYAzpy7w7lt3MTOzlMtTITupt7USyNp0JtJAC/WLXujaTdfMe4j/F1woZKdg\nIM2gYCDSIPH6RZfxBnMs/wNNheykmZQzEGmQxfpFwB7K/D4vsIN5+ne6CtlJ0ylnIACUgoDi7CzZ\nvj4yKwvxS00EpYCJwQkqM5XFtvN9O3jLI+/nndf1AuE0kYraSb1o05msKT89zeDEBCNPPcXgxAT5\n6elmd6ktVatf9Kv/6+28/zrn2FffYHDQtbJImkYjgw5XCgIGJyaYqSz9tZru6mLqwAGNEOokXr/o\n7MmzPH73T7ltdj8X6F68RiuLpB40MpBVFWdnSdnyn41eM4qzs03qUftLZVL039APQGGswM9mU/Sw\n/I8bHZEpjaZg0OGyfX0EK0ZZZXeyfX1N6lHnWEgoX86sVhZJ0ykYdLhMKkVuaIh0Vxf93d2ku7rI\nDQ1piqgBFg7Eia8s2skc6bSTy4XXqFSFNIpyBgI0bjWRVi0tt7AJzXqNs0EPqc9dzXv+3W6O//Uc\nv/vf06RSRhCEy05HR5vdW2l1Os9AEiE/Pc1YoUDKjMCd3NAQowMDze5W0ymhLI2iBLI0XSkIGCsU\nmKlUODc/z0ylwlihQCkImt21plNCWZJAwUAaQquW1qeEsjSTgoE0RL1XLZWCgMnz51t6pLFeQllT\nRFJPyhlIwyzkDHrNKNcwZ9BOuYhqCeX3ffRSBQKpCSWQJTFqvZqoHXdQxxPKqczSPZRKql0k26ME\nsiRGJpXihv7+mv2irmUuIilTTQsJ5XggePCrcwy+1Tl4o6t2kdSFgoG0tFrlIpJcrO8fv/oKv/07\nxswF4/xrxswMjI1pM5rUloKBtLRa7KCuxbLXeo0qglLA9z915uKlpj2upaZSUzrpTFre6MAAB/fu\n3XIuYmGqaSbWtjDVtJGvVc8E9mxxll9JBcxdWD4VduEC7NpVk28hAmhkIG1iO7mI7Uw1rTeq2O6I\noS/bxyVzweJS0x3MAw7uvO99rtyB1IyCgXS87Uw1rZXAXisPsdEgsXAgzkj6F3xt5ynCNVPGhXIX\nMzPG3Xe5cgdSE4mbJjKzm4GjhIEq5+6fb3KXpANsdapptVHFru7uxRHDwvTTWKHAwb17OXn27LrT\nSvEluAOjA+w9uJcLD71B6j875dh13cE8p085mZt6t3H3IgkbGZhZF/Bl4IPAvwRGzexdze2VdIqt\nTDWtNqr45fx81RHDqddeWzdZXW1Ekcqk2DNUIRj6JVyydO0cYfkKke1K2shgP3Da3acAzOwh4BDw\nQlN7JbKGaqOKUhBUHTEAayar4zmI+Iji/Nwcv7fjJ3T9iUOPkfrCVdhjA3ym90dccd07GnOj0tYS\nNTIA9gEvxv59JmoTSbSVo4rVRgzX7d69ZrK6Wg6ix4xP/fjHzHiFC7sc+ir4Z37Enw88ziceePOy\nzWkiW5W0kcGGHTlyZPH58PAww8PDTeuLSDWr5SFyQ0MX1Wha+Fy1HERQqZAy40KsPd3XxTueGGLg\nbZeGDaUSpcIpinsge+V1ZHaqXoXA+Pg44+PjG7o2UbWJzOwAcMTdb47+fS/gK5PIqk0krW6tGk0r\nC/p98R3v4Pd+8pPl9ZcMTr9/iH07/wXk8+T/5E7GbimTmocgnSL3ka8z+m4djSbLtUyhOjPrBgrA\njcDLwBPAqLs/v+I6BQNpayuDxUKA6GGOoBLwWfsyN9r3eNelX8Q+8GkGPzbLTGxBUbonzdSnpzRC\nkGXWCgaJmiZy93kz+yRwgqWlpc+v8zKRtpNJpZaNGEYHBvjArgp/+4ObGWCKPX6OisOZ73+KuUu7\nSc2zLBj00sWpl0+xN72X7J6sgoKsK1HBAMDdjwNDze6HSNLsnj/DtV1TzM+fW2wLfqWXt/8iIOhe\nfu1sJeDQN3+DHb07CCpz5A7lNG0ka0raaiIRqSIISpTLZ6lULixvv2SeNx/9M3KP9pIuQ/8spL0b\nD8rMesC54DVm5mYYe2SM0uvaqiyrS9zIQESWm57OUyiMYZYCKkAv3d1p3MsMDeXo+cAoozd/hIOF\nUxTtVc7e81vcdus852L/u3u7eii+WtR0kaxKwUAkwYKgRKEwRqUyA9E2tK6uNNde+y12776OVCr6\n5Z7JkMncRGZyktKrOwi6l9c8Ks8FZPdkF/9der1E8dWi8gmySNNEIgk2O1uMRgRLzHrp7d27FAji\nslky5+bIPcLStFEZcjd+afGXfv6ZPINHBxl5cITBo4Pkn1XpU0nY0tKN0tJS6RRBUGJiYjAaGYS6\nutIcODBVPRhAeCbm2Bil/m6KbyqT/W9fIvNbHwXCEcHg0UFm5pa+npahdo6WWVoqIsulUhmGhnJR\nzqB3MU+waiAAGB2FgwfJFItkslnIRNeWShSffJRUV8/y2kjdvconiEYGIkkTBCVmZ4v09WUXf+lX\na9uUhdHCJT0Mjr1WdYMaoDxCm2uZHcgbpWAg7Sq+csg9YGgox8DANvcHlEowOAgz4Xgg/24YOwS9\nO3dTjvYg4DB2bIxUd4pgPtC+hDalYCDSAraUH9iIyUkYGYFzS5vVSgO7KD74ZbL/+kMAyiN0iLWC\ngVYTiSTEaiuHZmeL2/vC2SysOF4zc36eG977ITJvQPHvwjxC3EIeQTqHgoFIQvT1ZXFf/kvbvUxf\nX3Z7XziTgVwO0mno7w8/5nJw8iQMDpK943cJXn9t2UvK82Wye7KUXi8x+dKkdi93AE0TiTRZPDl8\n9uzJi1YObTtnsKBUgmIxHCmA8ggdSDkDkYSqljDeu/fg9lYObYTyCB1JOQORBIqXmpifP0elMkOh\nMAZAf/8N9QsEsHYeYWeG4qtFUt3L8xfKI7Q3BQORJqlbwngjVssjRBvUsnuyBPMr6htFeQRpTwoG\nIk1St4TxRo2OwtRUmEiemgr/HcnszJA7lCPdk6Z/Rz/pnjS5QzlNEbUx5QxEGqxhCeMaUHXT9qIE\nskhCNC1hLIKCgUgi1G2HscgGaTWRSAI0NWHcINqk1roUDEQapOkJ4zrToTmtTdNEIg20lDNIZsJ4\nq3RoTmvQ4TYiCTEwMNqWCeOFTWrxYKBDc1qLgoFIg6VSmbYJAgu0Sa31KWcgItumTWqtTzkDEdma\neBXUqIyFNqklm/YZiEhtRWcqk0qFBe9yuWXlLCSZFAxEpHZWnKkMhIXupqYWRwiSTNp0JtIAQVDi\n/PlJgqDNN1wVi+GIIK63N2yXlqXVRCI1UK3mUDvsH6iqylkIlMtLJ6hJS9LIQGSbVjukpm1HCOuc\nhSCtSSMDkW1aqjm0NIe+UHOo3fYTLBodhYMHL1pNtBVagZQMCgYi29TuNYdWlclsezSQfybP2LEx\nUt0pgvmA3KEco+9u0+m1hNNqIpEaaNeaQ/WkekaNp9pEInXWrjWH6kn1jJJFwUCkRtqx5lA9qZ5R\nsmg1kYg0heoZJYtyBiLSVFpN1DgqRyEizVWlqJ00nspRiEjz5PNhLaORkfBjXsdhJpFGBiJSPypq\nlygaGYhsUccUn6sXFbVrGVpaKrKKjio+Vy8qatcyNDIQqaLjis/Vi4ratQyNDESq6Mjic/VSw6J2\nUj8KBiJVdGzxuXqpQVE7qS9NE4lUkUplGBrK0dWVpru7n66uNENDOY0KEqj0eonJlyYpva4pvO2o\n29JSMzsM/DbwStT0X939ePS5+4C7gTngU+5+Imq/Hvg60Ac86u6fXuVra2mpNEQQlFR8LsFUAntz\nmrIDOQoGr7n7/1jRfg3wTeAG4ArgJPBOd3czexz4pLtPmtmjwJfc/TtVvraCgUiHUwnszWvmPoNq\n3/QQ8JC7z7l7ETgN7Dezy4Hd7j4ZXfcN4NY6909EWtRCCey4hRLYsnn1DgafNLMnzexrZnZJ1LYP\neDF2zUtR2z7gTKz9TNQmInIRlcCurW2tJjKz7wID8SbAgc8B9wN/EE3//CHwp8A92/l+cUeOHFl8\nPjw8zPDwcK2+tIi0gIUS2GOPjNHb3Ut5vqwS2CuMj48zPj6+oWsbUpvIzAaBv3H3f2Vm9wLu7p+P\nPnccOAxMAY+5+zVR++3Av3H3j1X5esoZiAigEtib0ZScQZQDWPAR4Nno+THgdjNLmdmVwFXAE+7+\nc+Ccme03MwPuAB6pV/9EpD1kdma4Yd8NCgTbVM9NZ18ws/cCFaAIfBTA3Z8zs4eB54Ay8PHYn/mf\nYPnS0uN17J+ItAKdhdAQKmEtIsmVz8PYWFj5NAjCukaj2kewVTrpTERaj85CqDmdZyAirUdnITSU\ngoF0BB1S04J0FkJDKRhI25uezjMxMchTT40wMTHI9LTO4G0JOguhoZQzkLYWBCUmJgapVJbmnbu6\n0hw4MKXCc61Cq4lqZq2cgc4zkLamQ2ragM5CaAhNE0lb0yE1IhujYCBtTYfUiGyMcgbSEXRIjYg2\nnYmICNp0JiKyrk4/S1nBQEQ6Xv6ZPINHBxl5cITBo4Pkn+28vSiaJhKRjtZJZylrmkhEZBU6Szmk\nYCAiHU1nKYcUDESkoy2cpZzuSdO/o590T7ojz1JWzkBEhM44S1n7DERERAlkERFZm4KBiHSGUgkm\nJ8OPchEFAxFpf/l8eJ7yyEj4Md95m8rWo5yBiLS3UikMADNLm8pIp2FqquPOSVDOQFqWzi6WbSsW\nIbV8Uxm9vWG7LNJJZ5JY09N5CoUxzFK4BwwN5RgYGG12t6TVZLMQLN9URrkctssijQwkkYKgRKEw\nRqUyw/z8OSqVGQqFMY0QZPMyGcjlwqmh/v7wYy7XcVNE69HIQBJJZxdLTY2OwsGD4dRQNqtAUIWC\ngSSSzi6WmstkFATWoGkiSSSdXSzSWFpaKomms4tFake1iURERPsMRERkbQoGIiKiYCAiIgoGIiKC\ngoGIiKBgICIiKBiIiAgKBiIigoKBiIigYCAiIigYiIgICgYiIoKCgYiIoGAgIiIoGIiICAoGIiLC\nNoOBmf17M3vWzObN7PoVn7vPzE6b2fNmdlOs/Xoze9rMfmRmR2PtKTN7KHrN35vZ27bTNxER2bjt\njgyeAf4t8H/jjWZ2DXAbcA1wC3C/mS2crvMVYMzdrwauNrMPRu1jwD+7+zuBo8AXttm3mhkfH292\nF+pK99e62vneQPfXSNsKBu5ecPfTwMpj1A4BD7n7nLsXgdPAfjO7HNjt7pPRdd8Abo295oHo+V8B\nN26nb7WUpDesHnR/raud7w10f41Ur5zBPuDF2L9fitr2AWdi7WeitmWvcfd54FUze3Od+iciIjE9\n611gZt8FBuJNgAOfc/e/qVfHuHi0ISIi9eLu234AjwHXx/59L/DZ2L+PA78KXA48H2u/HfhK/Jro\neTfwyhrfz/XQQw899Nj8Y7Xfq+uODDYh/pf8MeAvzeyLhNM/VwFPuLub2Tkz2w9MAncAfxZ7zZ3A\n48B/AL632jdyd40aRERqaFvBwMxuBf4n8Bbgb83sSXe/xd2fM7OHgeeAMvBxj/6kBz4BfB3oAx51\n9+NRew540MxOA78gHDWIiEgD2NLvaBER6VQdvwN5tY1zZjZoZm+Y2Q+jx/2xz7XMxrlO2hhoZofN\n7EzsPbs59rlN3WsrMLObzeyFqO+fbXZ/tsLMimb2lJmdMrMnora9ZnbCzApm9h0zuyR2fdX3MSnM\nLGdm02b2dKxt0/fTlJ/LWiSQW/kBDAHvJMxRxJPgg8DTq7zmceCG6PmjwAej5x8D7o+e/0fCvRZJ\nvb9rgFOEU4VZ4McsjRRb5v5W3Oth4L9Uad/0vSb9QfiH3I+jn9Ne4EngXc3u1xbu46fA3hVtnwc+\nEz3/LPDH0fNrV3sfk/IAfh14b/x3x1bupxk/lx0/MvDVN85Rra3VNs6tcX9ttTEwptr7uJV7Tbr9\nwGl3n3L3MvAQ4X22GuPiGYr4z9kDLL0nH6bK+9iITm6Uu38fOLuieVP306yfy44PBuvIRtMNj5nZ\nr0dt7bJxrl03Bn7SzJ40s6/FhuNbudekW3lPrdT3OAe+a2aTZnZP1Dbg7tMA7v5z4LKofbX3Meku\n2+T9NOXnspZLSxNrixvnfga8zd3PRnPt3zazazf7rTff283rpI2Ba90rcD/wB+7uZvaHwJ8C91z8\nVSRBfs3dXzazDHDCzAqE72dcu61ySeT9dEQwcPeRLbymTDTcc/cfmtlPgKsJo/dbY5deEbUR+9zP\nzKwb6Hf3f95O3zfY103fH6vfR+LuL24T9/rnwEIg3Mq9Jt1LQDyB30p9X+TuL0cfS2b2bcJpn2kz\nG3D36WjK5JXo8lZ9vzZ7P025T00TLbf4l66ZvcXMuqLnbyfcOPfTaJh3zsz2R5VY7wAeiV62sHEO\n1tk41yQrNwbeHq0QupKljYEte3/Rf7QFHwGejZ5v5V6TbhK4Klr1liLcl3OsyX3aFDN7k5ntip7v\nBG4irIR8DLgruuxOlv/8XfQ+NrTTG2Nc/H/truj5uvfTtJ/LZmffm/0gTMy8CMwALwP/J2pf+GXy\nQ+AHwIdir3kf4Q/taeBLsfYdwMNR+wSQTer9RZ+7j3AFw/PATa14fyvu9RvA04Qra75NOPe8pXtt\nhQdwM1CI+n5vs/uzhf5fGb1Xp6L34N6o/c3AyejeTgB71nsfk/IAvkk4xXwB+CfgPwF7N3s/zfi5\n1KYzERHRNJGIiCgYiIgICgYiIoKCgYiIoGAgIiIoGIiICAoGIiKCgoGIiAD/H8HLC7XSxlKfAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f3275d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGzZJREFUeJzt3WuQpNV93/Hvf2anZ9rLDlrHzZACqZsIMYJIKi0Umy3b\nVZkUDLpUSUuUBDN+IRGaFLYgAuUFAvvFUo5TllTlCOQUgpK7LKTyNsIJgXWKwGolxlVOZcRYWnGp\nRe11lB4taxieipfdBXr2mcs/L55nZntme+59ebr796nq6p7TT/ecUz3b/z3nfy7m7oiISHfraXUF\nRESk9RQMREREwUBERBQMREQEBQMREUHBQEREqFMwMLOCmU2b2ctVZbvN7LCZlczseTO7uOq5B8zs\nuJm9ZmY3VZVfa2Yvm9nfmtlD9aibiIisr149gz8DPrGi7H7giLsPAz8CHgAws2uAW4CrgU8Bj5iZ\nxa/5FpB396uAq8xs5XuKiEgD1CUYuPtfA6dWFO8HHo8fPw7cHD/+LPCEu8+5exk4Duw1s0uBXe4+\nGV/33arXiIhIAzUyZ3CJu08DuPubwCVx+WXAiarrTsZllwGvV5W/HpeJiEiDNTOBrH0vREQSakcD\n33vazIbcfToeAnorLj8JvL/qusvjstXKL2BmCiwiIlvg7larvJ49A4tviw4Bt8WPvwA8U1V+q5ml\nzOwK4ErgxXgo6bSZ7Y0Typ+ves0F3L1ptwMHDjT19zX7pva1762T26b21f+2lrr0DMzsIDAC/CMz\n+yVwAPgq8BdmdjswRTSDCHc/ZmZPAseAWeCLfr6WdwHfAQaAZ939uXrUT0RE1laXYODuv73KUzeu\ncv0fAX9Uo/wnwEfrUScREdk4rUDegJGRkVZXoaHUvvbVyW0Dta+ZbL1xpCQyM2/HeouItJKZ4U1I\nIIuISJtSMBAREQUDERFRMBARERQMREQEBQMREUHBQEREUDAQEREUDKTDBAFMTkb3IrJxCgbSMYpF\nyGZhdDS6LxZbXSOR9qHtKKQjBEEUACqV82XpNExNQSbTunqJJIm2o5COVy5DKrW8rK8vKheR9SkY\nSEfI5SAMl5fNzkblIrI+BQPpCJkMFArR0NDgYHRfKGiISGSjlDOQjhIE0dBQLqdAILLSWjkDBQNp\nC2EQMlOeYSA3QCqTWv8FInKBtYJBXY69FGmk6eI0pXwJSxkeOsOFYYbGhlpdLZGOop6BJFoYhExk\nJ1ioLCyV9aR72De1Tz0EkU3S1FJpWzPlGSy1/G/X+oyZ8kyLaiTSmRQMJNEGcgN4uLwX6LPOQG6g\nRTUS6UwKBpJoqUyK4cIwPekeegd76Un3MFwY1hCRSJ0pZyBtQbOJRLZPU0tFREQJZBERWZuCgYiI\nKBiIiIiCgYiIoGAgIiIoGIiICAoGIiKCgoHIkiCAycnoXqTbKBiIAMUiZLMwOhrdF4utrpFIc2kF\nsnS9IIgCQKVyviydhqkpnZYmnUUrkEXWUC5DasV2R319UblIt1AwkK6Xy0EYLi+bnY3KRbqFgoF0\njTAIOTN5hjBY/s2fyUChEA0NDQ5G94WChoikuyhnIF1hI+coB0E0NJTLKRBIZ9IW1tLVdI6ySEQJ\nZOlqOkdZZH0KBtLxdI6yyPoUDKTj6RxlkfUpZyBdQ+coS7dTAllERJRAFhGRtSkYiIhI44OBmZXN\n7CUzO2pmL8Zlu83ssJmVzOx5M7u46voHzOy4mb1mZjc1un4iItKcnsECMOLue9x9b1x2P3DE3YeB\nHwEPAJjZNcAtwNXAp4BHzKzm+JaIiNRPM4KB1fg9+4HH48ePAzfHjz8LPOHuc+5eBo4DexERkYZq\nRjBw4AdmNmlmd8RlQ+4+DeDubwKXxOWXASeqXnsyLhNJHJ2MJp1kRxN+x2+4+xtmlgEOm1mJKEBU\n2/Q80QcffHDp8cjICCMjI9upo8imFIuQz0fnIIRhtMvp2FirayWy3Pj4OOPj4xu6tqnrDMzsAPAO\ncAdRHmHazC4FXnD3q83sfsDd/Wvx9c8BB9z9xyveR+sMpGV0Mpq0q5atMzCzXzGzi+LHO4GbgFeA\nQ8Bt8WVfAJ6JHx8CbjWzlJldAVwJvNjIOopslk5Gk07U6GGiIeC/m5nHv+vP3f2wmf0N8KSZ3Q5M\nEc0gwt2PmdmTwDFgFviiugDSaiu3sdDJaNKJtB2FyBpWOxRnMWfQ1xcFAuUMpB1obyKRLVjvUByd\njCbtZq1g0IzZRCJtaelQnKpE8eKhOKlMikxGQUA6h/YmElmFDsWRbqJgILIKHYoj3UQ5A5F16FAc\n6RRKIIuIiA63ERGRtSkYiIiIgoGIiCgYiDSUtrmWdqFgINIgxWK0u+noaHRfLLa6RiKr02wikQbQ\nNteSRJpNJNIkYRByZvIM//doqG2upa1obyKROqne4XQhdPbNDfM8Q0vPa5trSTL1DETqIAxCSvkS\nC5UF5k/P45UF7rMSQwMhg4PREFGhoCEiSS71DETqoNYOp30Dxot/McP07pS2uZbEUzAQqYPVdji9\ndM8AH1AQkDagYSKROtAOp9LuNLVUpI60w6kkmXYtFRERrTMQEZG1KRiIiIiCgUgraSM7SQoFA5EW\n0UZ2kiRKIIs0WRiEvHl0hr37B5ieOT/jSBvZSaMpgSySENPFaSayE/zicy/x+MwE/4Lppee0kZ20\nkoKBSJNU71/Eu/P0s8B9lLiYENBGdtJaCgYiTbK0f1GVeYwP7pzRRnbSctqbSAAIwpDyzAy5gQEy\nKzfil7qotX/RRWnnW08NcMWe6OfJSbSpnbSEegZCcXqa7MQEoy+9RHZiguL09Povkk2rtX/RhwvD\nfHwP/NVjZ/joB0LNLJKW0WyiLheEIdmJCSoLC0tl6Z4epvbtUw+hQar3Lzp15BQ/z5d4p2L04nyd\nYV5gSDOLpCE0m0hWVZ6ZIWXL/zb6zCjPzLSoRp0vlUkxeP0gAKV8Ca8ssJN5BqoSyppZJM2mYNDl\ncgMDhCt6WbPu5AYGWlSj7lEroTyHcSkzmlkkTadg0OUyqRSF4WHSPT0M9vaS7umhMDysIaImqJVQ\n7sN5e2CAQiH6WVtVSLMoZyBA82YTadbSctPFaUr5EtZnLMw6fb83zBX/ajf/67/NcPd/GmCmP0UY\nRtNOx8ZaXVtpdzrPQBKhOD1NvlQiZUboTmF4mLGhoVZXq+WUUJZmUTCQltOspfWFQchEdiJaoRyb\noYdb2YcPpjhyBK6/voUVlLan2UTScpq1tD4llKWVFAykKRo9aykIQybPnCEIw7q8Xyusl1DWEJE0\nkoKBNEUjZy11ygrqlSuULd3Drv84zCu/TCl5LA2nnIE0Vb1nE3ViLqI6oZzKnG9DEEQL0bR3kWyV\ncgaSGJlUiusHB+v2RV3PXERShpoWVyhXB4LHHgt4//snueGGQHsXSUMoGEhbq1cuIslDTY89VuR3\nfifLuXOjnD2bpVIpks9rMZrUl4KBtLV65CKCMCRfKlFZWOD0/DyVhQXypdKmegiN6lUEQcA99+SB\nCnA6vs/T2xto7yKpK51nIG1vbGiIG3fv3nIuYnGoqVJVtjjUtJH3auRiunK5TCqV4ty56trt4Ny5\nZ7nook8DSh5IfahnIB1hO7mI7Qw1rder2G6PIZfLMTe38rVn2WF3c911WYpKHkidKBhI19vOUNNa\nCey18hAbDRKZTIZCoUA6neaii3YtlVfCd6hUKtx+2+0ESh5IHSRuaqmZfRJ4iChQFdz9azWu0dRS\nqbutTHtdbWrrT667jut+8pOaU16PnDq17rBSEASUy2VyuRyZTIYgCDj0xCG+9KUv8R7vLV23k50c\nfv4wv37Tr2+z9dIN2mZqqZn1AP8F+ATwT4ExM/twa2sl3WIrQ02r9SremZ+v2WM4evbsusnqYrFI\nNptldHSUbDYaCspkMnxs6GPMMrvsPeeY41Iu3V7DRUheAnkvcNzdpwDM7AlgP/DzltZKZA21EthB\nGNbMQwBrJquDICCfz1OpVKhUoqvy+Txnzpzhy1/+Mj3x/99SpDCM+/ru4/I9lzelndLZEtUzAC4D\nTlT9/HpcJpJoK3sVq/UY9uzatWayenH2ULUdO3Zwzz33UKlUOMc5ABzn2/3f5q7H71q2OE1kq5LW\nM9iwBx98cOnxyMgIIyMjLauLSC2rTXktDA+TL5XoM2M2zhksPpfL5QhXJJXDMIynl55bKkvvTPPB\npz7I0E1RriEIAo4ePQrAnj17yGi/CgHGx8cZHx/f2MXunpgbsA94rurn+4Gv1LjORdrZW+fO+Yun\nT/tb585d8NzBgwc9nU774OCgp9Npf/TRRz2dTjuwdOvv7/djx44tXd/X17f0XCqV8oMHDza7SdIG\n4u/Omt+/iZpNZGa9QAm4AXgDeBEYc/fXVlznSaq3SL2tnE1ULBbJ5/O4OzMzM6TTaQC+8Y1vcO+9\n9zKzYi+mdDrN1NSUegiyTFuddBZPLX2Y81NLv1rjGgUD6TqvvfYae/bsWTZc1N/fT29vL++9996y\na3fu3MlTTz3F7t27lwKKyFrBIHE5A3d/DhhudT1Ekuadd95hYGBgWTBYmUtYFIYhn/nMZ+jv72du\nbo5CocCYDkWQNSRtNpGI1BAEAadOnbrgi39ubo5vfvOb9PX1LZX19vYyOztLGIacPXuWSqVCPp/X\nSmVZk4KBSMItLkK75ZZbWFhYoK+vj8HBQdLpNIVCgTvvvJOTJ0/y/PPP8/3vf5/e3t4L3qO3t5ey\ntjmVNSQuZ7ARyhlItwiCgGw2u7QADaLk8NNPP11zCunk5CQ33HADZ8+eXVbe39/PiRMnlq4P3g0o\nv10m974cmZ3KJ3SLttmOQkSWq7UIra+vj927d9dMCke7nM5dUP7www8vXV98pUj2oSyj3xsl+1CW\n4qva+VQUDEQSrdYitNnZWXK5XM3rq3c53bVrF/39/Tz66KPceeedQNQjyB/KU5mrcPrcaSpzFfLP\n5AneVT6h2ykYiCRY9Zd7dZ5gramiY2NjTE1N8cMf/pATJ04sBQKCgGD8Wf5xZfkkwr7ePspvlxvY\nCmkHyhmIJEwYBszMlBkYyJFKxWP8KxahbVqxCPk8C307mKmc5fbPwvc/Gj2V3pFm6t4pAOUROlxb\nLTrbCAUD6VTT00VKpTxmKdxDhocLDA1tc31AEEA2C1VJ6Pf64CP3XcSb6XkK+wvgkD+UJ9WbIpwP\nKewvMPYRrUvoNAoGIm0gDAMmJrIsLJz/0u7pSbNv39RSD2FLJidhdBROn14qWhjcxc8P/gmZkU8D\nkH0oS2WuasZS3FtQD6GzaDaRSBuYmSljtnzmkFkfMzPl7b1xLgcrktA9s3Ncs/fTZN5DeQQBFAxE\nEmNgIIf78i9t91kGBnLbe+NMBgoFSKdhcDC6LxTgyBHIZvnwb/97Xvn6WX7rlfMvmZ2fJfe+HMG7\nAZMnJzXbqAtomEikxaoTxqdOHYlzBn24z9YnZ7AoCKBcjnoKoDxCF1LOQCShaiWMd+++8YLZRHWn\nPEJXUs5AJIHCMKBUyrOwUGF+/jQLCxVKpTwAg4PXNy4QwNp5hJ0Zym+XSfWuWPmsPEJHUzAQaZGG\nJYw3YrU8QryGIfe+HOH8ipXPcR5BOpOCgUiLNCxhvFFjYzA1FSWSp6ain2OZnRkK+wukd6QZ7B8k\nvSNNYX9BQ0QdTDkDkSZrWsK4DrS7aWdRAlkkIVqWMBZBwUAkERq2wlhkgzSbSCQBWpowbhItUmtf\nCgYiTdLyhHGD6dCc9qZhIpEmOp8zSGbCeKuCdwMtUmsDaw0T7ahVKCKNMTQ01pEJ48VFatXBYHGR\nmoJBe1AwEGmyVCrTMUFgkRaptT/lDERk27RIrf0pZyAiW1O9C2q8jYUWqSWb1hmISH3FZyqTSkUb\n3hUKy7azkGRSMBCR+qlxpjLpdLS/UUa9gSTTojORJgjDgDNnJgnDDl9wVS5HPYJqfX1RubQtzSYS\nqYNaew51wvqBmmqchcDs7PkT1KQtqWcgsk2rHVLTsT2Edc5CkPaknoHINp3fc+j8GPrinkOdtp5g\nydgY3HjjBbOJtkIzkJJBwUBkmzp9z6FVZTLb7g0UXymSP5Qn1ZsinA8p7C8w9pEOHV5LOM0mEqmD\nTt1zqJG0n1HzaW8ikQbr1D2HGkn7GSWLgoFInXTinkONpP2MkkWziUSkJbSfUbIoZyAiLaXZRM2j\n7ShEpLVqbGonzaftKESkdYrFaC+j0dHovqjjMJNIPQMRaRxtapco6hmIbFHXbD7XKNrUrm1oaqnI\nKrpq87lG0aZ2bUM9A5Eaum7zuUbRpnZtQz0DkRq6cvO5RqnjpnbSOAoGIjV07eZzjVKHTe2ksTRM\nJFJDKpVheLhAT0+a3t5BenrSDA8X1CtIoODdgMmTkwTvaghvOxo2tdTMDgD/DngrLvo9d38ufu4B\n4HZgDrjH3Q/H5dcC3wEGgGfd/d5V3ltTS6UpwjDQ5nMJpi2wN6clK5DjYHDW3f/zivKrgYPA9cDl\nwBHgQ+7uZvZj4G53nzSzZ4GH3f35Gu+tYCDS5bQF9ua1cp1BrV+6H3jC3efcvQwcB/aa2aXALnef\njK/7LnBzg+snIm1qcQvsaotbYMvmNToY3G1mPzOzPzWzi+Oyy4ATVdecjMsuA16vKn89LhMRuYC2\nwK6vbc0mMrMfAEPVRYADvw88AvxBPPzzh8AfA3ds5/dVe/DBB5cej4yMMDIyUq+3FpE2sLgFdv6Z\nPH29fczOz2oL7BXGx8cZHx/f0LVN2ZvIzLLAX7r7x8zsfsDd/Wvxc88BB4Ap4AV3vzouvxX45+7+\nuzXeTzkDEQG0BfZmtCRnEOcAFn0OeDV+fAi41cxSZnYFcCXworu/CZw2s71mZsDngWcaVT8R6QyZ\nnRmuv+x6BYJtauSis6+b2ceBBaAM3Ang7sfM7EngGDALfLHqv/l3sXxq6XMNrJ+ItAOdhdAU2sJa\nRJKrWIR8Ptr5NAyjfY3GtI5gq3TSmYi0H52FUHc6z0BE2o/OQmgqBQPpCjqkpg3pLISmUjCQjjc9\nXWRiIstLL40yMZFlelpn8LYFnYXQVMoZSEcLw4CJiSwLC+fHnXt60uzbN6WN59qFZhPVzVo5A51n\nIB1Nh9R0AJ2F0BQaJpKOpkNqRDZGwUA6mg6pEdkY5QykK+iQGhEtOhMREbToTERkXd1+lrKCgYh0\nveIrRbIPZRn93ijZh7IUX+2+tSgaJhKRrtZNZylrmEhEZBU6SzmiYCAiXU1nKUcUDESkqy2epZze\nkWawf5D0jnRXnqWsnIGICN1xlrLWGYiIiBLIIiKyNgUDEekOQQCTk9G9XEDBQEQ6X7EYnac8Ohrd\nF7tvUdl6lDMQkc4WBFEAqJxfVEY6DVNTXXdOgnIG0rZ0drFsW7kMqeWLyujri8pliU46k8Sani5S\nKuUxS+EeMjxcYGhorNXVknaTy0G4fFEZs7NRuSxRz0ASKQwDSqU8CwsV5udPs7BQoVTKq4cgm5fJ\nQKEQDQ0NDkb3hULXDRGtRz0DSSSdXSx1NTYGN94YDQ3lcgoENSgYSCLp7GKpu0xGQWANGiaSRNLZ\nxSLNpamlkmg6u1ikfrQ3kYiIaJ2BiIisTcFAREQUDERERMFARERQMBARERQMREQEBQMREUHBQERE\nUDAQEREUDEREBAUDERFBwUBERFAwEBERFAxERAQFAxERQcFARETYZjAws39tZq+a2byZXbviuQfM\n7LiZvWZmN1WVX2tmL5vZ35rZQ1XlKTN7In7N/zazD2ynbiIisnHb7Rm8AvxL4K+qC83sauAW4Grg\nU8AjZrZ4us63gLy7XwVcZWafiMvzwD+4+4eAh4Cvb7NudTM+Pt7qKjSU2te+OrltoPY107aCgbuX\n3P04sPIYtf3AE+4+5+5l4Diw18wuBXa5+2R83XeBm6te83j8+L8CN2ynbvWUpA+sEdS+9tXJbQO1\nr5kalTO4DDhR9fPJuOwy4PWq8tfjsmWvcfd54G0z+9UG1U9ERKrsWO8CM/sBMFRdBDjw++7+l42q\nGBf2NkREpFHcfds34AXg2qqf7we+UvXzc8A/Ay4FXqsqvxX4VvU18eNe4K01fp/rpptuuum2+dtq\n36vr9gw2ofp/8oeAPzezbxAN/1wJvOjubmanzWwvMAl8Hvhm1Wu+APwY+DfAj1b7Re6uXoOISB1t\nKxiY2c3AnwC/BvwPM/uZu3/K3Y+Z2ZPAMWAW+KLH/6UH7gK+AwwAz7r7c3F5AfiemR0H/h9Rr0FE\nRJrAzn9Hi4hIt+r6FcirLZwzs6yZvWdmP41vj1Q91zYL57ppYaCZHTCz16s+s09WPbeptrYDM/uk\nmf08rvtXWl2frTCzspm9ZGZHzezFuGy3mR02s5KZPW9mF1ddX/NzTAozK5jZtJm9XFW26fa05O+y\nHgnkdr4Bw8CHiHIU1UnwLPDyKq/5MXB9/PhZ4BPx498FHokf/xbRWouktu9q4CjRUGEO+DvO9xTb\npn0r2noA+A81yjfd1qTfiP4j93fx32kf8DPgw62u1xba8Qtg94qyrwH3xY+/Anw1fnzNap9jUm7A\nbwIfr/7u2Ep7WvF32fU9A1994Ry1ytpt4dwa7euohYFVan2OW2lr0u0Fjrv7lLvPAk8QtbPdGBeO\nUFT/nT3O+c/ks9T4HJtRyY1y978GTq0o3lR7WvV32fXBYB25eLjhBTP7zbisUxbOderCwLvN7Gdm\n9qdV3fGttDXpVrapnepezYEfmNmkmd0Rlw25+zSAu78JXBKXr/Y5Jt0lm2xPS/4u6zm1NLG2uHDu\n74EPuPupeKz9aTO7ZrO/evO13bxuWhi4VluBR4A/cHc3sz8E/hi448J3kQT5DXd/w8wywGEzKxF9\nntU6bZZLItvTFcHA3Ue38JpZ4u6eu//UzP4PcBVR9H5/1aWXx2VUPff3ZtYLDLr7P2yn7hus66bb\nx+rtSFz7qm2ird8GFgPhVtqadCeB6gR+O9V9ibu/Ed8HZvY00bDPtJkNuft0PGTyVnx5u35em21P\nS9qpYaLllv6na2a/ZmY98eN/QrRw7hdxN++0me2Nd2L9PPBM/LLFhXOwzsK5Flm5MPDWeIbQFZxf\nGNi27Yv/oS36HPBq/HgrbU26SeDKeNZbimhdzqEW12lTzOxXzOyi+PFO4CainZAPAbfFl32B5X9/\nF3yOTa30xhgX/lu7LX68bnta9nfZ6ux7q29EiZkTQAV4A/ifcfnil8lPgb8BPl31muuI/miPAw9X\nlfcDT8blE0Auqe2Ln3uAaAbDa8BN7di+FW39LvAy0cyap4nGnrfU1na4AZ8ESnHd7291fbZQ/yvi\nz+po/BncH5f/KnAkbtth4H3rfY5JuQEHiYaYzwG/BP4tsHuz7WnF36UWnYmIiIaJREREwUBERFAw\nEBERFAxERAQFAxERQcFARERQMBARERQMREQE+P+YM1x37sMS3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1135c8f50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFvlJREFUeJzt3X2MZWV9wPHvz52Z3XGXtWsdlgRxhgbRtbZxMLttqkkn\ngUXwD6GkpUz/EMrQ+Bpf+keF9o/dVBNfEotYgyHtWLDRWbGhiA0FnMA0MWllxFU0IK6aGQFxval0\nxbjuzLJP/7jn7p65e+f9vpxz7/eT3HDmmXPnnoc7e3/zPL/n95xIKSFJ6m0v6fQFSJI6z2AgSTIY\nSJIMBpIkDAaSJAwGkiSaFAwiYjIijkbE47m2XRHxUEQ8FREPRsTLct+7JSKORMSTEXF5rv2SiHg8\nIn4QEZ9qxrVJklbXrJHBvwBvqWu7GZhOKb0GeBi4BSAiXgdcC+wBrgRuj4jInvNZYCKldDFwcUTU\n/0xJUgs0JRiklL4OPF/XfBVwV3Z8F3B1dvw24FBK6WRKaQ44AuyLiPOAc1JKs9l5n889R5LUQq3M\nGZybUjoKkFL6GXBu1n4+8HTuvGeztvOBZ3Ltz2RtkqQWa2cC2X0vJKmg+lr4s49GxO6U0tFsCujn\nWfuzwAW5816ZtS3XfpaIMLBI0gaklKJRezNHBpE9au4DbsiOrwe+kmu/LiIGIuJC4CLg0Wwq6VhE\n7MsSym/PPecsKaW2PQ4cONDW12v3w/6V99HNfbN/zX+spCkjg4j4IjAG/HZE/AQ4AHwM+HJE3AjM\nU11BRErpiYi4G3gCWATenc5c5XuAO4FtwP0ppQeacX2SpJU1JRiklP5imW9dtsz5HwU+2qD9MeD3\nmnFNkqS1swJ5DcbGxjp9CS1l/8qrm/sG9q+dYrV5pCKKiFTG65akTooIUhsSyJKkkjIYSJIMBpIk\ng4EkCYOBJAmDgSQJg4EkCYOBJAmDgbpMpVJhdnaWSqXS6UuRSsVgoK4xNTXF8PAw+/fvZ3h4mKmp\nqU5fklQabkehrlCpVBgeHub48eOn2wYHB5mfn2doaKiDVyYVh9tRqOvNzc0xMDCwpK2/v5+5ubnO\nXJBUMgYDdYWRkREWFhaWtC0uLjIyMtKZC5JKxmCgrjA0NMTk5CSDg4Ps3LmTwcFBJicnnSKS1sic\ngbpKpVJhbm6OkZERA4FUZ6WcgcFApeCHvLR5JpBVai4ZlVrPkYEKzSWjUvM4MlBpuWRUag+DgQrN\nJaNSexgMVGguGZXaw5yBSsHVRNLmubRUkmQCWZK0MoOBJMlgIEkyGEiSMBhIkjAYSJIwGEiSMBhI\np1UqFWZnZ6lUKp2+FKntDAYSbpMtWYGsnuc22eoVViBLK3CbbMlgILlNtoTBQD1kuQSx22RL5gzU\nI6amppiYmGBgYICFhQUmJycZHx9fco7bZKvbuYW1epoJYqnKBLJ6mgliaXUGA3U9E8TS6gwG6nom\niKXVmTNQzzBBrF5nAlmSZAJZkrQyg4EkqfXBICLmIuI7EXE4Ih7N2nZFxEMR8VREPBgRL8udf0tE\nHImIJyPi8lZfnySpPSODU8BYSmk0pbQva7sZmE4pvQZ4GLgFICJeB1wL7AGuBG6PiIbzW5Kk5mlH\nMIgGr3MVcFd2fBdwdXb8NuBQSulkSmkOOALsQ5LUUu0IBgn4WkTMRsRNWdvulNJRgJTSz4Bzs/bz\ngadzz302a5MKxzujqZv0teE13pRSei4ihoCHIuIpqgEib93rRA8ePHj6eGxsjLGxsc1co7Qua9n4\nTuq0mZkZZmZm1nRuW+sMIuIA8CvgJqp5hKMRcR7wSEppT0TcDKSU0sez8x8ADqSUvlH3c6wzUMe4\n8Z3KqmN1BhHx0ojYkR1vBy4HvgvcB9yQnXY98JXs+D7guogYiIgLgYuAR1t5jdJ6ufGdulGrp4l2\nA/8eESl7rS+klB6KiG8Cd0fEjcA81RVEpJSeiIi7gSeAReDdDgHUafXbWLjxnbqR21FIK1guN1Br\n7+/vZ3Fx0ZyBSsG9iaQNWC034MZ3KpuVgkE7VhNJpVTLDeSDQS03MDQ0dPohdQP3JpKWYW5AvcRg\nIC3Dm+Kol5gzkFZhbkDdwgSyJMmb20iSVmYwkCQZDCRJBgOppdzmWmVhMJBaZGpqiuHhYfbv38/w\n8DBTU1OdviRpWa4mklrAba5VRG5HIbVJrSbh+eefX3ErC6loDAZSk+R3OD1x4gSnTp1a8n23slCR\nOU0kNUGjaaH+/n76+vrc5lqF4TSR1GKNdjgdHBzky1/+Mrt27XIrCxWewUBqguV2OB0dHTUIqBRc\nWio1gTucquzMGUhN5A6nKjJ3LZUkuWupJGllBgNJksFA6iQ3slNRGAykDnEjOxWJCWSpzSqVCocP\nH+bqq692Izu1lQlkqSBqo4FrrrlmSSCAMxvZSZ3gyEBqk0b7F+U5MlCrOTKQCqC2f1G97du3W7Gs\njnNvIgFWzrZDo/2Ltm3bxj333MPo6CgAs7OzvgfqCEcGclVLmzTav+hzn/sco6Oj3HHHHb4H6ihz\nBj3O2zO2X34UNj09zY033shvfvObJef4HqgVzBloWY3msV3V0lpDQ0Ps3bsXgImJibMCAfgeqP0M\nBj1uuX34vT1j6y2XUAbfA7WfwaDHuQ9/5zQKxMDp9wBwqwq1jcFAjI+PMz8/z/T0NPPz8y29T697\n8ZxRH4i3bdvGhz/8YR577DF+9KMfmVBWW5lAVttMTU0xMTHBwMAACwsL3iA+Y0JZ7eLNbdRxrlpa\n3UoVyjt37mR6evp04lnaCFcTqeNctbQ6E8rqJIOB2qLVq5a6IRexWkLZEZRayWCgtmjlqqVuqaBe\nLqHc6qS+BOYM1GbN3gOpG3MRy/0/cv8obZY5AxVGrfq2WR9mzcxFFGWqqdH/o3vuuINrLriA6y69\ntNSjHxWXwUCl1qxcRJGnmn55xx1c8c538tUTJ/juCy/wtuPHmZiY6HjQUndxmkilV6tf6O/vZ3Fx\ncd31C82YamrZFE6lwosXXMCWEydON/0aeP2OHXzp4Yddaqp1cZpIXW2zFdSbnWpq6ahibo6ou7aT\nwKUnTnDhjh3Nex31PEcG6nmbGRms9txNjxgqFRgehtzPT8DCwAADW7YQk5PgSiOtkSMDaQWbWfa6\n0qhipRHDmpPVQ0MwOQmDg5zavp0EBLB1YYE4fpyTN9xQDRjSJhVuZBARVwCfohqoJlNKH29wjiMD\nNd1G/opfbmTw2GOP8cY3vrHhiGF6enr1PZoqFZibg5GRakCoVPjloUPwvvexM3faMeDUgw+y6/LL\nN9xv9Y7SjAwi4iXAZ4C3AL8LjEfEazt7VeoVG1n2utyo4le/+lXDEcPhw4eZmJjg+PHjHDt2jOON\nVgZNTVWnhvbvr/53agqGhpjfvZv6zSr6gbkN91g6o6/TF1BnH3AkpTQPEBGHgKuA73f0qqQVjI+P\nc9llly0ZVVQqlYZLXgEGBgaWjBhq00pD2QiAiYlqjqB2zsQE/PKXvP6DH+Q41ZxB7dnv7O/nk6Oj\nre+kul6hRgbA+cDTua+fydqkQqsfVSw3YhgdHV25LmJuDuo3q+vrg/e/nzh+nJdSzRm8BPijrVu5\n8q67rEZWUxRtZLBmBw8ePH08NjbG2NhYx65FaqTRiAFgcnLyrLqI0x/oIyNQv1ndwkI1QORqDfq2\nb+eRe+45kyuoVODw4erx6Gg1z6CeNzMzw8zMzJrOLVQCOSL+EDiYUroi+/pmINUnkU0gq+xWTFZP\nTVWnhvr7YXERbr0VPvjBJctL2bq1+uG/Z0/1/Ouvr54L1cBx550uOdVZSnNzm4jYAjwFXAo8BzwK\njKeUnqw7z2Cg7la/mqgWIKAaFAYHq8e33gof+ADU3RmNwUGYn3eEoCVKEwzg9NLS2ziztPRjDc4x\nGKj3PPlkdQooN13E1q2wZQv8+tdLz92+He65B3btOhNQ1PNKFQzWwmCgnjQ7W11ueuzYmbYdO6o5\nhfo8Q38/RFSDxcmT1cI1p416XmnqDCQto1KB559fOioAePFF+PSnqx/+NVu2VPMHCwvwwgvVaaWJ\nCSuVtSKDgVR0tSK0a6+FU6eqH/w7d1bzApOT8I53wLPPwoMPwpe+VA0G9fr6qjkIaRlOE0lF1mCj\nOgYH4d57Gy8hnZ2FSy+tjgjytm6Fp58+fb53TetNThNJZdWoCK2/v5oYbvQhPjJSzRHUu+220+cX\n+UY+6hxHBlKRLTcyWGnZaG0Zai13cNtt1akkuvOe0Vo7RwZSWeW2sF6SJ1jpg3t8vBosHn64OjWU\nCwT3338/fX1LNx7Y6D2j1V0cGUgF03A+v74IbZ1qtwbt6+vjhbp8Qm1kAJhH6HLWGUglUfvQXvFe\nB+vUaGoI4JxzzuHkyZNMTk4CNP11VTwGA6kEWjWfPzs7y/79+zmWK1bbsWMHn/nMZ3jrW98KYB6h\nR5gzkEpgpVtobsbIyMhZ22a/+OKLpwOBeQSBIwOpMFq50qc2/ZTfNhswj9BjnCaSCiyfMK7dHzn/\nod2sufv868DZU0NgHqHbGQykgmqUMG50Q5xmM4/QmwwGUgF1sgBstdduFCx27tzJ9PQ0e/fubem1\nqXVMIEsF1KqE8Vosd4/mWhBqlHRecq9mdR1HBlKHFGFriJU2rGuUdDZnUG5OE0kF0q6EcTO4u2l3\nWSkY9DVqlNQajRLG8/Pzhf3AHRoaKtw1qTUcGUhtUoRpIfU2E8hSAXQyYdwulUqF2dlZKt5is3QM\nBlKbdPsKHW+aU25OE0lt1K0rdJwCKwcTyFJBjI+Pt6XCuN1qU2D5YFCbAuuWPnY7g4HUZt24Qqfb\np8B6gTkDSZu2WkWzis+cgaQNaVSQZpFasVmBLKmpWnF7TrWewUBS07hyqLwsOpPaoFcKrnqheK4X\nGQykJuilgitXDnUnp4mkTerFaZNuLZ7rdhadSS3UiwVXzSyecwVSMThNJG1Sr06bDA0NsXfv3k19\ngPfS9FrROU0kNYHTJuvXi9NrneY0kdRi3brnUCv14vRakRkMpCbpxj2HWqlXp9eKypyBpI5wP6Ni\nMWcgqaNcTdQ+bkchqaP8wC8Gt6OQ1DEuHy0HRwaSWsblo8XiyEDaoF7ZfK5V3NSuPAwG0jKc3tg8\nl4+Wh9NEUgNObzSP1dnFYQWytE5WxzaP1dnlYDCQGnB6o7mszi4+cwZSA1bHlodJ/uZoWc4gIg4A\nfwX8PGv625TSA9n3bgFuBE4C708pPZS1XwLcCWwD7k8pfWCZn23OQG1hsVSx1fIRAwMDLCwsmI9Y\nRUcqkLNg8EJK6R/q2vcAXwT2Aq8EpoFXp5RSRHwDeG9KaTYi7gduSyk92OBnGwykHmeSf/06WWfQ\n6EWvAg6llE6mlOaAI8C+iDgPOCelNJud93ng6hZfn6SSsoahuVodDN4bEd+OiH+OiJdlbecDT+fO\neTZrOx94Jtf+TNYmSWcxyd9cm1pNFBFfA3bnm4AE/B1wO/D32fTPR4BPAjdt5vXyDh48ePp4bGyM\nsbGxZv1oSSVQS/LX1zA4RXTGzMwMMzMzazq3LUVnETEMfDWl9PsRcTOQUkofz773AHAAmAceSSnt\nydqvA/44pfSuBj/PnIEkwCT/enQkZ5DlAGquAb6XHd8HXBcRAxFxIXAR8GhK6WfAsYjYFxEBvB34\nSquuT1J3GBoaYu/evQaCTWpl0dknIuINwClgDngHQErpiYi4G3gCWATenfsz/z0sXVr6QAuvT1IJ\n+Jd/e7g3kaTCso6gubzTmaTSsY6g+byfgaTSsY6gvQwG6gnuX1M+1hG0l8FAXc+b1JSTmwW2lzkD\ndTXnncvP1UTN481t1LO8SU35eS+E9nCaSF3NeWdpbQwG6mrOO0trY85APcF5Z8miM0kSFp1J0qp6\nvRbFYCCp51mL4jSRpB7XS7UoThNJ0jLcA6nKYCCpp1mLUmUwkNTTrEWpMmcgSfRGLYp1BpIkE8iS\npJUZDCT1hF4vKluNwUBS17OobHXmDCR1tV4qKluNOQOVlkN7bZZFZWtjMFBhObRXM1hUtjZOE6mQ\nHNqrmaamppiYmKC/v5/FxUUmJycZHx/v9GW1nfdAVul472I10/j4OJdddlnXF5VthsFAheTQXs02\nNDRkEFiBOQMVkvvFSO1lzkCF1gv7xUjt4t5EkiTrDCRJKzMYSJIMBpIkg4EkCYOBJAmDgSQJg4Ek\nCYOBJAmDgSQJg4EkCYOBJAmDgSQJg4EkCYOBJAmDgSQJg4EkiU0Gg4j404j4XkS8GBGX1H3vlog4\nEhFPRsTlufZLIuLxiPhBRHwq1z4QEYey5/x3RLxqM9cmSVq7zY4Mvgv8CfBf+caI2ANcC+wBrgRu\nj4ja3XU+C0yklC4GLo6It2TtE8AvUkqvBj4FfGKT19Y0MzMznb6ElrJ/5dXNfQP7106bCgYppadS\nSkeA+tuoXQUcSimdTCnNAUeAfRFxHnBOSmk2O+/zwNW559yVHf8bcOlmrq2ZivSGtYL9K69u7hvY\nv3ZqVc7gfODp3NfPZm3nA8/k2p/J2pY8J6X0IvB/EfHyFl2fJCmnb7UTIuJrwO58E5CAv0spfbVV\nF8bZow1JUquklDb9AB4BLsl9fTPwodzXDwB/AJwHPJlrvw74bP6c7HgL8PMVXi/58OHDh4/1P5b7\nXF11ZLAO+b/k7wO+EBG3Up3+uQh4NKWUIuJYROwDZoG3A5/OPed64BvAnwEPL/dCKSVHDZLURJsK\nBhFxNfCPwCuA/4iIb6eUrkwpPRERdwNPAIvAu1P2Jz3wHuBOYBtwf0rpgax9EvjXiDgC/C/VUYMk\nqQ3izGe0JKlX9XwF8nKFcxExHBG/johvZY/bc98rTeFcLxUGRsSBiHgm955dkfveuvpaBhFxRUR8\nP7v2D3X6ejYiIuYi4jsRcTgiHs3adkXEQxHxVEQ8GBEvy53f8H0sioiYjIijEfF4rm3d/enI72Uz\nEshlfgCvAV5NNUeRT4IPA48v85xvAHuz4/uBt2TH7wJuz47/nGqtRVH7twc4THWqcAT4IWdGiqXp\nX11fDwB/3aB93X0t+oPqH3I/zH5P+4FvA6/t9HVtoB8/BnbVtX0c+Jvs+EPAx7Lj1y33PhblAbwZ\neEP+s2Mj/enE72XPjwzS8oVzNGorW+HcCv3rqsLAnEbv40b6WnT7gCMppfmU0iJwiGo/yyY4e4Yi\n/3t2F2fek7fR4H1sx0WuVUrp68Dzdc3r6k+nfi97PhisYiSbbngkIt6ctXVL4Vy3Fga+NyK+HRH/\nnBuOb6SvRVffpzJde14CvhYRsxFxU9a2O6V0FCCl9DPg3Kx9ufex6M5dZ3868nvZzKWlhbXBwrmf\nAq9KKT2fzbXfGxGvW+9Lr/9q16+XCgNX6itwO/D3KaUUER8BPgncdPZPUYG8KaX0XEQMAQ9FxFNU\n38+8blvlUsj+9EQwSCnt38BzFsmGeymlb0XEj4CLqUbvC3KnvjJrI/e9n0bEFmBnSukXm7n2NV7r\nuvvH8v0oXP/y1tHXfwJqgXAjfS26Z4F8Ar9M135aSum57L+ViLiX6rTP0YjYnVI6mk2Z/Dw7vazv\n13r705F+Ok201Om/dCPiFRHxkuz4d6gWzv04G+Ydi4h92U6sbwe+kj2tVjgHqxTOdUh9YeB12Qqh\nCzlTGFja/mX/0GquAb6XHW+kr0U3C1yUrXoboFqXc1+Hr2ldIuKlEbEjO94OXE51J+T7gBuy065n\n6e/fWe9jWy96bYKz/63dkB2v2p+O/V52Ovve6QfVxMzTwHHgOeA/s/bah8m3gG8Cb809541Uf2mP\nALfl2rcCd2ft/wOMFLV/2fduobqC4Ung8jL2r66vnwcep7qy5l6qc88b6msZHsAVwFPZtd/c6evZ\nwPVfmL1Xh7P34Oas/eXAdNa3h4DfWu19LMoD+CLVKeYTwE+AvwR2rbc/nfi9tOhMkuQ0kSTJYCBJ\nwmAgScJgIEnCYCBJwmAgScJgIEnCYCBJAv4fG/ROEQ3X6wQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e30e710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "53\n",
      "43\n",
      "0.811320754717\n"
     ]
    }
   ],
   "source": [
    "#matches=[4,3,5,0,1,2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2]\n",
    "matches=[3,1,0,5,2,4,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2,-2]\n",
    "hitmatches_sorted=[]\n",
    "\n",
    "sorting = np.argsort(particles)\n",
    "particles_sorted = particles[sorting]\n",
    "\n",
    "print sorting\n",
    "\n",
    "tracks_sorted = tracks[sorting]\n",
    "pixelx_sorted = pixelx[sorting]\n",
    "pixely_sorted = pixely[sorting]\n",
    "\n",
    "print particles_sorted\n",
    "print tracks_sorted\n",
    "\n",
    "cmap = ['b','m','c','y','g','r','k']\n",
    "for iparticle in np.unique(particles_sorted):\n",
    "#for itrack in [7]:\n",
    "    hits_particle = (particles_sorted == iparticle)\n",
    "#    print hitmatches\n",
    "    print iparticle\n",
    "    hitmatches_sorted = hitmatches_sorted + sum(hits_particle == 1)*[matches[iparticle]]\n",
    "    plt.scatter(pixelx_sorted[hits_particle],pixely_sorted[hits_particle],color=cmap[int(matches[iparticle]) % 7])\n",
    "axes = plt.gca()\n",
    "axes.set_xlim([-200,200])\n",
    "axes.set_ylim([-200,200])\n",
    "plt.axis('equal')\n",
    "plt.show()\n",
    "\n",
    "cmap = ['b','m','c','y','g','r','k','k','k','k','k']\n",
    "for itrack in np.unique(tracks_sorted):\n",
    "    if itrack >= 0:\n",
    "        hits_track = (tracks_sorted == itrack)\n",
    "        plt.scatter(pixelx_sorted[hits_track],pixely_sorted[hits_track],color=cmap[int(itrack)])\n",
    "axes = plt.gca()\n",
    "axes.set_xlim([-200,200])\n",
    "axes.set_ylim([-200,200])\n",
    "plt.axis('equal')\n",
    "plt.show()\n",
    "\n",
    "        \n",
    "\n",
    "\n",
    "plt.scatter(pixelx_sorted[hitmatches_sorted == tracks_sorted],pixely_sorted[hitmatches_sorted == tracks_sorted],color='k')\n",
    "plt.scatter(pixelx_sorted[hitmatches_sorted != tracks_sorted],pixely_sorted[hitmatches_sorted != tracks_sorted],color='r')\n",
    "\n",
    "axes = plt.gca()\n",
    "axes.set_xlim([-200,200])\n",
    "axes.set_ylim([-200,200])\n",
    "plt.axis('equal')\n",
    "plt.show()\n",
    "\n",
    "print len(tracks_sorted)\n",
    "print sum(hitmatches_sorted == tracks_sorted)\n",
    "print 43./53."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see which hits were not assigned to any cluster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cmap = ['k','b','y','g','r']\n",
    "itrack = -1\n",
    "hits_track = (tracks == itrack)\n",
    "plt.scatter(pixelx[hits_track],pixely[hits_track],color=cmap[int(itrack) % 5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scoring the assignment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The score function measures the fraction of hits with good track assignment. More precisely, it is the overlap fraction between the true vs predicted cluster_id, averaged over all events. However, notice that the reconstruction algorithm has the freedom to label the indices in a way different than the true cluster_ids, therefore the \"overlap\" should be more clearly defined.\n",
    "\n",
    "For this, we first determine a one-to-one matching between the true cluster_ids and predicted cluster_ids, particle-by-particle for a given event. For a given particle (defined by a single true cluster_id) the corresponding predicted cluster_id is the most frequent one among its hits. \n",
    "\n",
    "Once every particle is associated to a track in this way, we loop over to remove the duplicates, i.e. the tracks which contain the maximum hits of more than one particle. Among those particles, the particle with the most hits inside this track is considered to be the right match, and the other particles are considered to be lost.\n",
    "\n",
    "After the above procedure, the well-associated hits in the event are counted, and their ratio to the total number of hits is determined as the efficiency of this event.\n",
    "\n",
    "In the end, every event, regardless of its hit multiplicity, is considered equivalent and the efficiency is averaged over the events in the full sample.\n",
    "\n",
    "The input to the score function is the true and predicted arrays of [event_id, cluster_id]. Note that in order for the scoring to function, the order of the hits in the predicted array must match exactly the same order in the input. If your algorithm reshuffles the order of the hits, make sure you save the original indices before doing that, and re-sort the result array corresponding to those indices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def score_function(y_true, y_pred):\n",
    "    '''Compute a clustering score.\n",
    "\n",
    "    Cluster ids should be nonnegative integers. A negative integer\n",
    "    will mean that the corresponding point does not belong to any\n",
    "    cluster.\n",
    "\n",
    "    We first identify assigned clusters by taking the max count of \n",
    "    unique assigned ids for each true cluster. We remove all unassigned\n",
    "    clusters (all assigned ids are -1) and all duplicates (the same\n",
    "    assigned id has majority in several true clusters) except the one\n",
    "    with the largest count. We add the counts, then divide by the number\n",
    "    of events. The score should be between 0 and 1. \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    y_true : np.array, shape = (n, 2)\n",
    "        The ground truth.\n",
    "        first column: event_id\n",
    "        second column: cluster_id\n",
    "    y_pred : np.array, shape = n\n",
    "        The predicted cluster assignment (predicted cluster_id)\n",
    "    \"\"\"\n",
    "    '''\n",
    "    score = 0.\n",
    "    event_ids = y_true[:, 0]\n",
    "    y_true_cluster_ids = y_true[:, 1]\n",
    "    y_pred_cluster_ids = y_pred\n",
    "\n",
    "    unique_event_ids = np.unique(event_ids)\n",
    "    for event_id in unique_event_ids:\n",
    "        event_indices = (event_ids==event_id)\n",
    "        cluster_ids_true = y_true_cluster_ids[event_indices]\n",
    "        cluster_ids_pred = y_pred_cluster_ids[event_indices]\n",
    "\n",
    "        unique_cluster_ids = np.unique(cluster_ids_true)\n",
    "        n_cluster = len(unique_cluster_ids)\n",
    "        n_sample = len(cluster_ids_true)\n",
    "\n",
    "        # assigned_clusters[i] will be the predicted cluster id\n",
    "        # we assign (by majority) to true cluster i \n",
    "        assigned_clusters = np.empty(n_cluster, dtype='int64')\n",
    "        # true_positives[i] will be the number of points in \n",
    "        # predicted cluster[assigned_clusters[i]]\n",
    "        true_positives = np.full(n_cluster, fill_value=0, dtype='int64')\n",
    "        for i, cluster_id in enumerate(unique_cluster_ids):\n",
    "            # true points belonging to a cluster\n",
    "            true_points = cluster_ids_true[cluster_ids_true == cluster_id]\n",
    "            # predicted points belonging to a cluster\n",
    "            found_points = cluster_ids_pred[cluster_ids_true == cluster_id]\n",
    "            # nonnegative cluster_ids (negative ones are unassigned)\n",
    "            assigned_points = found_points[found_points >= 0]\n",
    "            # the unique nonnegative predicted cluster ids on true_cluster[i]\n",
    "            n_sub_cluster = len(np.unique(assigned_points))\n",
    "            # We find the largest predicted cluster in the true cluster.\n",
    "            if(n_sub_cluster > 0):\n",
    "                # sizes of predicted assigned cluster in true cluster[i]\n",
    "                predicted_cluster_sizes = np.bincount(\n",
    "                    assigned_points.astype(dtype='int64'))\n",
    "                # If there are ties, we assign the tre cluster to the predicted\n",
    "                # cluster with the smallest id (combined behavior of np.unique\n",
    "                # which sorts the ids and np.argmax which returns the first \n",
    "                # occurence of a tie).\n",
    "                assigned_clusters[i] = np.argmax(predicted_cluster_sizes)\n",
    "                true_positives[i] = len(\n",
    "                    found_points[found_points == assigned_clusters[i]])\n",
    "            # If none of the assigned ids are positive, the cluster is unassigned\n",
    "            # and true_positive = 0\n",
    "            else:\n",
    "                assigned_clusters[i] = -1\n",
    "                true_positives[i] = 0\n",
    "\n",
    "        # resolve duplicates and count good assignments\n",
    "        sorted = np.argsort(true_positives)\n",
    "        true_positives_sorted = true_positives[sorted]\n",
    "        assigned_clusters_sorted = assigned_clusters[sorted]\n",
    "        good_clusters = assigned_clusters_sorted >= 0\n",
    "        for i in range(len(assigned_clusters_sorted) - 1):\n",
    "            assigned_cluster = assigned_clusters_sorted[i]\n",
    "            # duplicates: only keep the last count (which is the largest\n",
    "            # because of sorting)\n",
    "            if assigned_cluster in assigned_clusters_sorted[i+1:]:\n",
    "                good_clusters[i] = False\n",
    "        n_good = np.sum(true_positives_sorted[good_clusters])\n",
    "        score += 1. * n_good / n_sample\n",
    "    score /= len(unique_event_ids)\n",
    "    return score\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Check the score of a single event:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y_pred = np.zeros((len(y_predicted),2))\n",
    "\n",
    "y_pred[:,1] = y_predicted\n",
    "y_pred[:,0] = [ievent] * len(y_predicted)\n",
    "\n",
    "y_test_event = y[y[:,0]==ievent]\n",
    "\n",
    "print y_test_event\n",
    "print particles\n",
    "print y_pred\n",
    "print tracks\n",
    "score_function(y_test_event,y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-validation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to be aware that when we split the dataset for cross-validation, we shuffle and split by the events (which correspond to many rows), not by rows. The function provided helps doing that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# We do a single fold because blending would not work anyway:\n",
    "# mean of cluster_ids make no sense\n",
    "def get_cv(y_train_array):\n",
    "    unique_event_ids = np.unique(y_train_array[:, 0])\n",
    "    event_cv = ShuffleSplit(\n",
    "        n_splits=1, test_size=0.5, random_state=57)\n",
    "    for train_event_is, test_event_is in event_cv.split(unique_event_ids):\n",
    "        train_is = np.where(\n",
    "            np.in1d(y_train_array[:, 0], unique_event_ids[train_event_is]))[0]\n",
    "        test_is = np.where(\n",
    "            np.in1d(y_train_array[:, 0], unique_event_ids[test_event_is]))[0]\n",
    "        yield train_is, test_is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unit testing\n",
    "\n",
    "It is <b><span style=\"color:red\">important that you test your submission files before submitting them</span></b>.\n",
    "For this we provide a unit test. \n",
    "Place the python file <code>clusterer.py</code> and <code>user_test_model.py</code></a> in a directory, set the paths to the data files in <code>user_test_model.py</code>, and run:\n",
    "\n",
    "<code>python user_test_model.py</code>\n",
    "\n",
    "**If it runs and prints**   \n",
    "<code>\n",
    "Reading file ...\n",
    "Training ...\n",
    "Testing ...\n",
    "score =  0.846205321725\n",
    "</code>\n",
    "**you can submit the code.**\n",
    "\n",
    "You can also run the unit test from the notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading file ...\n",
      "Training ...\n",
      "p[0]+x+p[1]*x+p[2]*x+p[3]*x\n",
      "Testing ...\n",
      "score =  0.477611940299\n"
     ]
    }
   ],
   "source": [
    "!python user_test_submission.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You're now ready to implement your own algorithm. Enjoy tracking!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}