{"cells":[{"metadata":{"_uuid":"3550939f8248c7d15b3cb87695b4d6a713da66da"},"cell_type":"markdown","source":"# TrackML Problem Explanation and Data Exploration"},{"metadata":{"_uuid":"52784d4de9ec81506da32c01a2a2d3975f23d4c8"},"cell_type":"markdown","source":"Code along of [Wesam Elshamy's Kaggle kernel](https://www.kaggle.com/wesamelshamy/trackml-problem-explanation-and-data-exploration)"},{"metadata":{"_uuid":"c6d3096d0d123d5c04b21a56ba2f05965dcf5241"},"cell_type":"markdown","source":"## 0. Problem Description"},{"metadata":{"_uuid":"1433355257a60759fcf0258e944cb92edc5d75f2"},"cell_type":"markdown","source":"Link every **track** to one **hit**.\n\nEvery particle leaves a track behind it. We want to link every track to a unique (max 1 per detector) set of hits.\n\nIn every **event**, a large number of **particles** are released. They move along a path leaving behind their **tracks**. They eventually **hit** a particle detector surface on the other end.\n\nIn the training data we have the following information on each **event**:\n\n- **Hits**: $x, y, z$ coords of each hit on the particle detector\n- **Particles**: Each particle's initial position ($v_x, v_y, v_y$), momentum ($p_x, p_y, p_z$), charge ($q$), and number of hits.\n- **Truth**: Mapping between hits and generating particles, particle trajectory, momentum, and hit weight.\n- **Cells**: Precise location of where each particle hit the detector and how much energy it deposited."},{"metadata":{"_uuid":"bd356c10b36b588fed3f31d8ae158ace7dd8e447"},"cell_type":"markdown","source":"## 1. Data Exploration:"},{"metadata":{"_uuid":"af83498fdc2c5ebfab10e82f23741a8ac4a8976e"},"cell_type":"markdown","source":"1. settings –> add a custom package –>GitHub user/repo (LAL/trackml-library)\n2. restart kernel"},{"metadata":{"trusted":true,"_uuid":"863328fab34d25b61afa1e8c42129fbf88ed061a"},"cell_type":"code","source":"%matplotlib inline\n\nimport os\nimport numpy as np\nimport pandas as pd\n\nfrom trackml.dataset import load_event\nfrom trackml.randomize import shuffle_hits\nfrom trackml.score import score_event\n\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nimport seaborn as sns\n\nfrom pathlib import Path","execution_count":5,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"be94863da0855d6e5e4338e56ea895e09a9a541b"},"cell_type":"code","source":"PATH = Path('../input/train_1')\nevent_prefix = 'event000001000'\nhits,cells,particles,truth = load_event(PATH/event_prefix)\n\nmem_bytes = (hits.memory_usage(index=True).sum()\n             + cells.memory_usage(index=True).sum()\n             + particles.memory_usage(index=True).sum()\n             + truth.memory_usage(index=True).sum())\nprint(f'{event_prefix} memory usage {mem_bytes/2**20:.2f} MB')","execution_count":10,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"c2668fcecc8739e72d0948aa9548ef6f58abc01c"},"cell_type":"markdown","source":"### 2.1. Hits Data"},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"997913bf5f0e0041cad8f483b34d70645ce240f5"},"cell_type":"markdown","source":"#### Where Did it Hit?\n\nHere we have the $x,y,z$ gloval coords [mm] of where the particles hit the detector surface."},{"metadata":{"trusted":true,"_uuid":"c41327e1fcd609f892be44400c654477816a4ea8"},"cell_type":"code","source":"hits.head()","execution_count":11,"outputs":[]},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"626a18929b8392e394107e8ad660b70ff290155a"},"cell_type":"markdown","source":"Here's the distribution of $x,y,z$ hit locations in event 1000. This is only for one of 8,850 events."},{"metadata":{"trusted":true,"collapsed":true,"_uuid":"a0d92d9246b3e260bb7c43ad6036dad12e3f9cb8"},"cell_type":"markdown","source":"#### Vertical Intersection ($x, y$) in Detection Layers\n\nAs shown in the figure below, the hits are semi-evenly distributed on the detector surface $x, y$. The white circle in the center of the plot is where the beam pipe lies. [Clarification](https://www.kaggle.com/wesamelshamy/trackml-problem-explanation-and-data-exploration/comments#323803): [agerom](https://www.kaggle.com/artemiosgeromitsos).\n\nThe colors represent different detector volumes. See [Joshua Bonatt's notebook](https://www.kaggle.com/jbonatt/trackml-eda-etc)."},{"metadata":{"trusted":true,"_uuid":"470778ae61803db17258c68da6b22d30a954eea8"},"cell_type":"code","source":"g = sns.jointplot(hits.x, hits.y, s=1, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n    vol = hits[hits.volume_id==volume]\n    plt.scatter(vol.x, vol.y, s=3, label='volume {}'.format(volume))\n    \nplt.xlabel('X [mm]'); plt.ylabel('Y [mm]'); plt.legend(); plt.show()","execution_count":14,"outputs":[]},{"metadata":{"_uuid":"961c822ccdb0ea3ab3a1c792e77c338112842449"},"cell_type":"markdown","source":"#### Horizontal Intersection ($y, z$) in Detection Layers\n\nYou can think of the chart below as a horizontal intersection in the detection surface, where every dot is a hit. Notice the relationship between the different activity levels in this char and the one above for $x, y$.\n\nAgain, the colors represent different volumes in the detector surface."},{"metadata":{"trusted":true,"_uuid":"8e0f66da3c3b7691f7d1cc767d78ad5f992bc513"},"cell_type":"code","source":"g = sns.jointplot(hits.z, hits.y, s=1, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nvolumes = hits.volume_id.unique()\nfor volume in volumes:\n    vol = hits[hits.volume_id==volume]\n    plt.scatter(vol.z, vol.y, s=3, label=f'volume {volume}')\n\nplt.xlabel('Z [mm]');plt.ylabel('Y [mm]');plt.legend();plt.show()","execution_count":17,"outputs":[]},{"metadata":{"_uuid":"128a631f6904935d960bc2d7a93fda95cc291e21"},"cell_type":"markdown","source":"And here is how the hits in this event look in 3D. Again, a sample from 1 event. This combines the previous 2 charts in 3D.\n\nNotice how th eparticles penetrate the detector surface along the $z$ coordinate:"},{"metadata":{"trusted":true,"_uuid":"d8040e052c4f67d81b1836813c1db54677d260b3"},"cell_type":"code","source":"fig = plt.figure(figsize=(12,12))\nax  = fig.add_subplot(111, projection='3d')\nfor volume in volumes:\n    vol = hits[hits.volume_id==volume]\n    ax.scatter(vol.z,vol.x,vol.y, s=1, label=f'volume {volume}', alpha=0.5)\nax.set_title('Hit Locations');ax.set_xlabel('Z [mm]');ax.set_ylabel('X [mm]')\nax.set_zlabel('Y [mm]'); plt.show()","execution_count":20,"outputs":[]},{"metadata":{"_uuid":"a96b5531c9f5aad723135e7fc437bb2525aa40c2"},"cell_type":"markdown","source":"#### Affected Surface Object\n\nThe **volume**, **layer**, and **module** are nested parts on the detector surface. The volume is made of layers, which in turn have modules. Analyzing their response could help us understand if some of them are dead/defective, and therefore we may need to account for the bias they cause.\n\nThe figure betlow shows a plot of every combination of `x`, `y`, `volume`, `layer`, and `module`. The colors identify different *volumes*. Along the main diagonal we have the variables' histograms.\n\nThe (`hit_id`, `x`) and (`hit_id`, `y`) pairs show us how different volumes are layered."},{"metadata":{"trusted":true,"_uuid":"dfdc33204bdbb26c261f995543e12ca3e4a027d6"},"cell_type":"code","source":"hits_sample = hits.sample(8000)\nsns.pairplot(hits_sample, hue='volume_id', size=8); plt.show()","execution_count":21,"outputs":[]},{"metadata":{"_uuid":"ea763e447b504231249764e61352585875e94cb8"},"cell_type":"markdown","source":"### 2.2 Particle Data"},{"metadata":{"_uuid":"7c834ab061f0fc8f53eac07f6f3ade896c9ccbbc"},"cell_type":"markdown","source":"The particle data help us undrestand each particle's intitial position, momentum, and charge, which we can join with the event truth dataset to get the particle's final position and momentum. This is needed to identify the tracks that each particle generated.\n\nThe data looks like this:"},{"metadata":{"trusted":true,"_uuid":"3a8c4b096c4437ff02a7d698bf6c548264607976"},"cell_type":"code","source":"particles.head()","execution_count":22,"outputs":[]},{"metadata":{"_uuid":"3b849a836d045d7836fb02647d46d892192ccb04"},"cell_type":"markdown","source":"#### Hit Rate and Charge Distribution\n\nLet's see the distribution of the number of hits per particle, show below. A significant number of particles had no attributed hits, and most of them have positive charge in this event:"},{"metadata":{"trusted":true,"_uuid":"04bf31682cfa562f952ace1640e301131f0bf56a"},"cell_type":"code","source":"plt.figure(figsize=(15,5))\nplt.subplot(1,2,1)\nsns.distplot(particles.nhits.values, axlabel='Hits/Particle', bins=50)\nplt.title('Distribution of number of hits per particle for event 1000.')\nplt.subplot(1,2,2)\nplt.pie(particles.groupby('q')['vx'].count(),\n        labels=['negative', 'positive'],\n        autopct='%.0f%%', shadow=True, radius=0.8)\nplt.title('Distribution of particle charges.'); plt.show()","execution_count":24,"outputs":[]},{"metadata":{"_uuid":"05d892e01acf914edb3dfeec8ab07723767b54f4"},"cell_type":"markdown","source":"#### Initial Position and Momentum\n\nLet's now take a look at the initial position of the particles around the global coordinates' origin $(x, y) = (0, 0)$, as shown in the figure below.\n\nThe initial position distribution is more concentrated around the origin (less variance) than its hit position (shown above under the Hits Data section). As the particles hit th edetection surface, they tend to scatter as shown in the particle trafjectory plot at the end of teh notebook.\n\nThe colors here show the number of hits for each particle."},{"metadata":{"trusted":true,"_uuid":"bb02d2496a2cb0a9a6a51a1c2fee1984f088388b"},"cell_type":"code","source":"g = sns.jointplot(particles.vx, particles.vy, s=3, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nn_hits = particles.nhits.unique()\nfor n_hit in n_hits:\n    p = particles[particles.nhits==n_hit]\n    plt.scatter(p.vx, p.vy, s=3, label=f'Hits {n_hit}')\n\nplt.xlabel('X [mm]'); plt.ylabel('Y [mm]'); plt.legend(); plt.show()","execution_count":27,"outputs":[]},{"metadata":{"_uuid":"9b3650ab0dfd90e0e175115e3a51b8ca11c92883"},"cell_type":"markdown","source":"And here's the initial position of the particles in a $z, y$ view. Colors show number of hits."},{"metadata":{"trusted":true,"_uuid":"fac0d59831dca05404c38e9c4d70a460097d39fc"},"cell_type":"code","source":"g = sns.jointplot(particles.vz, particles.vy, s=3, size=12)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nn_hits = particles.nhits.unique()\nfor n_hit in n_hits:\n    p = particles[particles.nhits == n_hit]\n    plt.scatter(p.vz, p.vy, s=3, label=f'Hits {n_hit}')\n\nplt.xlabel('Z [mm]'); plt.ylabel('Y [mm]'); plt.legend(); plt.show()","execution_count":28,"outputs":[]},{"metadata":{"_uuid":"9de87012d99f44d14de820046c4b7b74b02405cf"},"cell_type":"markdown","source":"And this is what they look like in 3D:"},{"metadata":{"trusted":true,"_uuid":"1801de3d72c14c8feafcf644ade788b58575cbf3"},"cell_type":"code","source":"fig = plt.figure(figsize=(12,12))\nax = fig.add_subplot(111, projection='3d')\nfor charge in [-1,1]:\n    q = particles[particles.q==charge]\n    ax.scatter(q.vz, q.vx, q.vy, s=1, label=f'Charge {charge}', alpha=0.5)\n\nax.set_title('Sample of 1000 Particle initial locations')\nax.set_xlabel('Z [mm]'); ax.set_ylabel('X [mm]'); ax.set_zlabel('Y [mm]')\nax.legend(); plt.show()","execution_count":29,"outputs":[]},{"metadata":{"_uuid":"4a0391fc825cb0cc5f068e987b3dc43ad0fa34e2"},"cell_type":"markdown","source":"#### Pair plot\n\nLet's now take a look at the relationship between different combinations of the particle variables. Again, the colors represent the number of hits.\n\nThere's no large skew in the distribution of the number of hits over other variables. It looks like the particles are targetted towards the global origin $(x, y) = (0, 0)$ and are evently distributed aroudn it."},{"metadata":{"trusted":true,"_uuid":"45b4a076cf24534f5692376f43fa3ebfdbf85542"},"cell_type":"code","source":"p_sample = particles.sample(8000)\nsns.pairplot(p_sample, vars=['particle_id', 'vx', 'vy', 'vz', 'px', 'py', 'pz',\n                             'nhits'], hue='nhits', size=8)\nplt.show()","execution_count":30,"outputs":[]},{"metadata":{"_uuid":"e72c5b0ab0101c38bb2dcc74e85712b19b4478ae"},"cell_type":"markdown","source":"#### Particle Trajectory\n\nWe can reconstruct the trajectories for a few particles given their intersection points with the detector layers. As explained in the [competition evaluation page](https://www.kaggle.com/c/trackml-particle-identification#evaluation), hits from straight tracksshave larger weights, and random tracks or hits with very short tracks have weights of zero. The figure below shows 2 such exampmles.\n\nSee [Makahana's notebook for trajectory plotting](https://www.kaggle.com/makahana/quick-trajectory-plot)."},{"metadata":{"trusted":true,"_uuid":"ba3def904f4dc85727159b140894c63ed6c03516"},"cell_type":"code","source":"# Get particle ID with max number of hits in this event\nparticle0 = particles.loc[particles.nhits==particles.nhits.max()].iloc[0]\nparticle1 = particles.loc[particles.nhits==particles.nhits.max()].iloc[1]\n\n# Get points where the same particle intersected subsequent layers of the observation material\np_traj_surface0 = truth[truth.particle_id==particle0.particle_id][['tx','ty','tz']]\np_traj_surface1 = truth[truth.particle_id==particle1.particle_id][['tx','ty','tz']]\n\np_traj0 = (p_traj_surface0.append({'tx':particle0.vx, \n                                   'ty':particle0.vy,\n                                   'tz':particle0.vz}, \n                                  ignore_index=True).sort_values(by='tz'))\np_traj1 = (p_traj_surface1.append({'tx':particle1.vx,\n                                   'ty':particle1.vy,\n                                   'tz':particle1.vz},\n                                  ignore_index=True).sort_values(by='tz'))\n\nfig = plt.figure(figsize=(10,10))\nax  = fig.add_subplot(111, projection='3d')\n\nax.plot(xs=p_traj0.tx, ys=p_traj0.ty, zs=p_traj0.tz, marker='o')\nax.plot(xs=p_traj1.tx, ys=p_traj1.ty, zs=p_traj1.tz, marker='o')\n\nax.set_xlabel('X [mm]'); ax.set_ylabel('Y [mm]'); ax.set_zlabel('Z [mm] –– Detection layers')\nplt.title('Trajectories of 2 particles as they cross the detection surface ($Z$ axis).')\nplt.show()","execution_count":33,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.5","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}