{"cells":[{"metadata":{"_uuid":"c2b7ca346409732793bdfdb75e477e97214247e3","_cell_guid":"053c84e4-ee69-462f-85a8-cbb79cae1e0c"},"cell_type":"markdown","source":"# Introduction and table of contents\n___\n\nIn this notebook, I made some exploratory analysis in order to better understand what we are dealing with in this competition. The trajectories (both in position and momentum space) of some particles are visualized using the Truth dataset. There are still lots of possible explorations and I'll be adding more analyses as the competetion goes on. \n\n## Contents\n\n-[Hits dataset](#Hits-dataset)\n\n-[Truth dataset](Truth-dataset)\n\n-[Particles dataset](Particles-dataset)"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"trusted":true},"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D\nimport matplotlib.gridspec as gridspec\n%matplotlib inline\n\nfrom trackml.dataset import load_event\nfrom trackml.randomize import shuffle_hits\nfrom trackml.score import score_event","execution_count":1,"outputs":[]},{"metadata":{"_uuid":"ccb292153fb5e8308ba35ce1327e762585a93a40"},"cell_type":"markdown","source":"## Load event 1000"},{"metadata":{"_cell_guid":"674c96c1-3444-4c35-a395-d5bfd17082d9","_uuid":"da7672d80d228666d234190d7bda38e9e7f53304","collapsed":true,"trusted":true},"cell_type":"code","source":"hits, cells, particles, truth = load_event('../input/train_1/event000001000')","execution_count":2,"outputs":[]},{"metadata":{"_uuid":"c94c21f5da9f2c978511b292174e8cd7e6f609b0","_cell_guid":"6421bc22-d7b2-40a7-8dad-1ae459061f70"},"cell_type":"markdown","source":"## Exploring the datasets\n\n### Hits dataset"},{"metadata":{"_uuid":"21a80ee9718be3344f17c69f664f89d66115b06e","_cell_guid":"6292d015-4642-4f76-b184-1fdd4d50d297","trusted":false,"collapsed":true},"cell_type":"code","source":"sns.jointplot(x = \"x\", y= \"y\", data = hits, alpha = 0.05, size = 8);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7f35c0d32c6376bf126680bea63e3fc459ee22da","_cell_guid":"16dcc781-1d89-4546-bfd3-f9b09b2ab570"},"cell_type":"markdown","source":"Here I set alpha = 0.05 so that each blue point represents 20 actual points. With this alpha number, it is possible to see where the particles are mostly concentrated. They concentric rings, with a few particles outside them (few compared to the actual number of particles, of course)."},{"metadata":{"_uuid":"df71088eea77ceb422a42dc1e0404c7e5c9d7fa9","_cell_guid":"72f96a1f-a9b9-4c91-b9dd-0b0110a58cec","trusted":true},"cell_type":"code","source":"sns.jointplot(x = \"x\", y= \"z\", data = hits, alpha = 0.05, size = 8);","execution_count":3,"outputs":[]},{"metadata":{"_uuid":"1d628d9a236d50e6c7cd4041172cc4701da49010","_cell_guid":"0d5bedc6-3bd6-4896-89ce-97709abb8468","trusted":false,"collapsed":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8, 8))\nax = fig.add_subplot(111, projection='3d')\nax.scatter(\n    xs=hits.x.values,\n    ys=hits.y.values,\n    zs=hits.z.values,\n    alpha = 0.05\n)\nax.set_title('Hit Locations of event 1000')\nax.set_xlabel('X (millimeters)')\nax.set_ylabel('Y (millimeters)')\nax.set_zlabel('Z (millimeters)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a43213fc7c29ab5d01b580d7fc71b3665026b30e","_cell_guid":"6eef6086-c59e-4a19-8fa6-1f627e0c391c","trusted":false,"collapsed":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8, 8))\nax = fig.add_subplot(111, projection='3d')\nax.scatter(\n    xs=truth.tpx.values,\n    ys=truth.tpy.values,\n    zs=truth.tpz.values,\n    alpha = 0.01\n)\nax.set_title('Momentum space of event 1000')\nax.set_xlabel('Px (GeV/c)')\nax.set_ylabel('Py (GeV/c)')\nax.set_zlabel('Pz (GeV/c)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"27366f5453ac09a4251aa4c627b71600b7490a38","_cell_guid":"073aed63-ed51-48f7-bd17-8b79d822aab5"},"cell_type":"markdown","source":"When alpha = 0.01, we can see that lots of particles have high momenta pointing towards z-direction (there are two clear points on the north and south pole of the \"momentum sphere\")."},{"metadata":{"_uuid":"6a3804605d97ca70e6fbfe1358e9a8a7c8343bab","_cell_guid":"dafe7e54-386c-4873-8a23-b7d5c99c7c9e","trusted":false,"collapsed":true},"cell_type":"code","source":"hits_per_layer = pd.DataFrame(hits.groupby(\"layer_id\", as_index=False)[\"hit_id\"].count())\nhits_per_layer = hits_per_layer.rename(columns = {\"hit_id\": \"num_of_hits\"})\n\nsns.barplot(x = hits_per_layer.layer_id, y = hits_per_layer.num_of_hits);","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"5cefe795d60c87e90142742718e6ee99364692bb","_cell_guid":"52325b55-597c-4e97-aabd-57bb38685650"},"cell_type":"markdown","source":"Some layers have way less number of hits than others. Let's see how the cross section is looking at the different layers (based on [Joshua Bonatt's kernel](https://www.kaggle.com/jbonatt/trackml-eda-etc))"},{"metadata":{"_uuid":"d45c83d73f3340c9675a77e45cbce01c71e48971","_cell_guid":"e78f36a3-f6e3-4b12-af3e-9b315cf0194c","trusted":false,"collapsed":true},"cell_type":"code","source":"g = sns.jointplot(hits.x, hits.y,  s=1, size=10)\ng.ax_joint.cla()\nplt.sca(g.ax_joint)\n\nlayers = hits.layer_id.unique()\nfor layer in layers:\n    lay_hit = hits[hits.layer_id == layer]\n    plt.scatter(lay_hit.x, lay_hit.y, s=1, label='layer {}'.format(layer))\n\nplt.xlabel('x (mm)')\nplt.ylabel('y (mm)')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"2387ccdb61a6e1f88947c58d12dfb9237d971978","_cell_guid":"5293e536-df10-42cd-b555-0822259cb833"},"cell_type":"markdown","source":"Note that the layers doesn't have a specific order, as there are lots of colored points all over the cross section. However, there are concentric circles with apparently only one type of layer."},{"metadata":{"_kg_hide-input":true,"_uuid":"381ef8a1752137513b6901ea860fbc2384204693","_cell_guid":"c6e56a97-5687-418b-8dc1-a2ec1086be96","trusted":false,"collapsed":true},"cell_type":"code","source":"fig = plt.figure(figsize=(8, 8))\nax = fig.add_subplot(111, projection='3d')\n\n#get the unique layer values\nlayers = hits.layer_id.unique()\n\n#for each layer\nfor layer in layers:\n    #get the data belonging to that layer alone\n    lay_hit = hits[hits.layer_id == layer]\n    #make a scatterplot of only points of a specific layer\n    #and give them colors (using label)\n    ax.scatter(lay_hit.x, lay_hit.y, lay_hit.z, s=1, label='layer {}'.format(layer))\n\n#set axes names\nax.set_xlabel('x (mm)')\nax.set_ylabel('y (mm)')\nax.set_zlabel('z (mm)')\nax.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_uuid":"46e099d36f8aec6c1dc4bbf3c8ba105efa44e407"},"cell_type":"code","source":"hits_per_vol = pd.DataFrame(hits.groupby(\"volume_id\", as_index=False)[\"hit_id\"].count())\nhits_per_vol = hits_per_vol.rename(columns = {\"hit_id\": \"num_of_hits\"})\n\nsns.barplot(x = hits_per_vol.volume_id, y = hits_per_vol.num_of_hits);","execution_count":44,"outputs":[]},{"metadata":{"_uuid":"b3c2a25b5a93282f2b465776959ee05e5b732dc0"},"cell_type":"markdown","source":"Similar to the layers, some volumes were more hit than others. Let's make a scatterplot to actually \"see\" the different volumes."},{"metadata":{"trusted":true,"_uuid":"8e29164aa8b36de091d4e9413bce886af00cfa60"},"cell_type":"code","source":"fig = plt.figure(figsize=(25, 8))\ngs = gridspec.GridSpec(nrows=2, ncols=3, left=0.05, right=0.48, wspace=0.5, hspace = 0.3)\nax = fig.add_subplot(gs[0:2,0:2], projection = '3d')\nax2 =fig.add_subplot(gs[0,2])\nax3 =fig.add_subplot(gs[1,2])\n\n#get the unique volume values\nvolumes = hits.volume_id.unique()\n\n#for each volume\nfor volume in volumes:\n    #get the data belonging to that volume alone\n    vol_hit = hits[hits.volume_id == volume]\n    #make a scatterplot of only points of a specific volume\n    #and give them colors (using label)\n    ax.scatter(\n        vol_hit.x, \n        vol_hit.y, \n        vol_hit.z, \n        s=0.1, \n        label='volume {}'.format(volume))\n\n    ax2.scatter( \n    x = vol_hit.x,\n    y = vol_hit.y,\n    s = 0.1,\n    label='volume {}'.format(volume))\n        \n    ax3.scatter(\n    x = vol_hit.x,\n    y = vol_hit.z,\n    s = 0.1,\n    label='volume {}'.format(volume))\n    \n#set axes names\nax.set_xlabel('x (mm)')\nax.set_ylabel('y (mm)')\nax.set_zlabel('z (mm)')\nax.set_title('Colored volumes')\nax.legend(loc = \"upper left\")\n\nax2.set_xlabel('x (mm)')\nax2.set_ylabel('y (mm)')\nax2.set_title('Colored volumes x-y cross section')\n\nax3.set_xlabel('x (mm)')\nax3.set_ylabel('z (mm)')\nax3.set_title('Colored volumes x-z cross section')\n\nplt.show()","execution_count":57,"outputs":[]},{"metadata":{"_uuid":"bdcb864eef4efdf3424921e1ea98e7fb35380921"},"cell_type":"markdown","source":"Here I made the points very tiny so it becomes better to see where each volume is located. "},{"metadata":{"_uuid":"69723263f0d6a1a6dfd52ce9886cff5ee2d488b7","_cell_guid":"0ae9a2d9-f179-4645-aaa6-29b0f0d9d6f7"},"cell_type":"markdown","source":"### Truth dataset\n\nFor starters, let's see some trajectories from the truth dataset."},{"metadata":{"_kg_hide-input":false,"_uuid":"b41f7e2f10a5115bb1b9e3f2cfd264a9eed703e3","_cell_guid":"f0680d76-b6e9-4e79-b018-60b480025da0","trusted":false,"collapsed":true},"cell_type":"code","source":"#get the information for some particles\ntruth_0 = truth[truth.particle_id == particles.iloc[20,0]]\ntruth_1 = truth[truth.particle_id == particles.iloc[10,0]]\ntruth_2 = truth[truth.particle_id == particles.iloc[5,0]]\n\n#create figure instance\nfig = plt.figure(figsize=(8, 8))\nax = fig.add_subplot(111, projection='3d')\n\n#plot each particle's path\nax.plot(\n    xs=truth_0.tx,\n    ys=truth_0.ty,\n    zs=truth_0.tz, marker='o')\nax.plot(\n    xs=truth_1.tx,\n    ys=truth_1.ty,\n    zs=truth_1.tz, marker='o')\nax.plot(\n    xs=truth_2.tx,\n    ys=truth_2.ty,\n    zs=truth_2.tz, marker='o')\n\nax.set_title('Trajectories of 3 different particle_id')\nax.set_xlabel('x (mm)')\nax.set_ylabel('y (mm)')\nax.set_zlabel('z (mm)')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"7c9b04fd4d98d2145acacf61ba03d12ab264176d","_cell_guid":"5fbb3b97-8cc0-4e2a-a116-5184e563c9b8"},"cell_type":"markdown","source":"Now let's automate this a bit and see a larger number of trajectories. Not too large so we can still see some distinct paths."},{"metadata":{"_kg_hide-input":true,"_uuid":"d90295cb4c66e6117a2767da886c9255394dbff8","_cell_guid":"50f670a9-8e1f-411f-a680-2b40dea50a66","trusted":true},"cell_type":"code","source":"#get the information of a given particle from the truth dataframe\ninformation = []\n\n#append the information of each desired particle on a list\nfor i in np.arange(0,20,1):\n    #select the true values for a given particle_id\n    particle_information = truth[truth.particle_id == particles.iloc[i,0]]\n    #append on the list\n    information.append(particle_information)\n\n#create figure instance\nfig = plt.figure(figsize=(25, 8))\ngs = gridspec.GridSpec(nrows=2, ncols=3, left=0.05, right=0.48, wspace=0.5, hspace = 0.3)\nax = fig.add_subplot(gs[0:2,0:2], projection = '3d')\nax2 =fig.add_subplot(gs[0,2])\nax3 =fig.add_subplot(gs[1,2])\n\n#plot the trajectory for each particle on the information list\nfor trajectory in information:\n    \n    ax.plot(\n    xs=trajectory.tx,\n    ys=trajectory.ty,\n    zs=trajectory.tz, marker='o')\n    \n    ax2.scatter( \n    x = trajectory.tx,\n    y = trajectory.ty)\n    \n    \n    ax3.scatter(\n    x= trajectory.tx,\n    y = trajectory.tz)\n    \n#labels\nax.set_xlabel('x (mm)')\nax.set_ylabel('y (mm)')\nax.set_zlabel('z (mm)')\nax.set_title('20 different trajectories')\n\nax2.set_xlabel('x (mm)')\nax2.set_ylabel('y (mm)')\nax2.set_title('Detector x-y cross section')\n\nax3.set_xlabel('x (mm)')\nax3.set_ylabel('z (mm)')\nax3.set_title('Detector x-z cross section')\nplt.show()","execution_count":39,"outputs":[]},{"metadata":{"_uuid":"0c846d53d2010883b67bc93e0950790e829af78b","_cell_guid":"157da532-941b-4ae1-ae55-21a85a6962df"},"cell_type":"markdown","source":"I think this image is really beautiful. All particles are moving from the center towards all directions, changing their paths according to their charges and other properties on the instant they were created, such as initial energy and momentum.\n\nNow a similar plot but in the momentum space (I decided to leave the first two blocks of code - the ones that store information about the particles we will plot - because one might want to see a different number of \"momentum trajectories\" than what they saw in position space."},{"metadata":{"_kg_hide-input":true,"_uuid":"3e7d081654c8f554841a0a05c8e2dbf8767a3b88","_cell_guid":"5881e012-3bc5-440d-a2b2-60fb7cce7277","trusted":true},"cell_type":"code","source":"#get the information of a given particle from the truth dataframe\ninformation = []\n\n#append the information of each desired particle on a list\nfor i in np.arange(0,20,1):\n    #select the true values for a given particle_id\n    particle_information = truth[truth.particle_id == particles.iloc[i,0]]\n    #append on the list\n    information.append(particle_information)\n\n#create figure instance\nimport matplotlib.gridspec as gridspec\nfig = plt.figure(figsize=(25, 8))\ngs = gridspec.GridSpec(nrows=2, ncols=3, left=0.05, right=0.48, wspace=0.3, hspace = 0.3)\nax = fig.add_subplot(gs[0:2,0:2], projection = '3d')\nax2 =fig.add_subplot(gs[0,2])\nax3 =fig.add_subplot(gs[1,2])\n\n#plot the trajectory for each particle on the information list\nfor trajectory in information:\n    \n    ax.plot(\n    xs=trajectory.tpx,\n    ys=trajectory.tpy,\n    zs=trajectory.tpz, marker='o')\n    \n    ax2.scatter( \n    x = trajectory.tpx,\n    y = trajectory.tpy)\n    \n    ax3.scatter(\n    x= trajectory.tpx,\n    y = trajectory.tpz)\n\n#labels\nax.set_xlabel('Px (GeV/c)')\nax.set_ylabel('Py (GeV/c)')\nax.set_zlabel('Pz (GeV/c)')\nax.set_title(\"Trajectories in momentum space\")\n\nax2.set_xlabel('Px (GeV/c)')\nax2.set_ylabel('Py (GeV/c)')\nax2.set_title (\"Trajectories in momentum space x-y cross section\")\n\nax3.set_xlabel('Px (GeV/c)')\nax3.set_ylabel('Pz (GeV/c)')\nax3.set_title (\"Trajectories in momentum space x-z cross section\")\nplt.show()","execution_count":42,"outputs":[]},{"metadata":{"_uuid":"fc13d3033beaee4fa6c2c3a487afb688bb31db1a","_cell_guid":"9d280e66-03cc-442b-9a90-1fd704f4f963"},"cell_type":"markdown","source":"The momentum plot didn't look as pretty as the position one due to one outlier. This particle had too much momentum in all three directions. Over 12 GeV/c of magnitude in both x and y directions and about 40 GeV/c in z direction."},{"metadata":{"_cell_guid":"f234da63-f384-40b3-8583-731c412f66ff","_kg_hide-input":false,"_uuid":"7e98cdac2d2ce8878cba2adfdf3a196bbb42ec9f","_kg_hide-output":false,"trusted":false,"collapsed":true},"cell_type":"code","source":"plt.hist(truth.weight, bins = 100)\nplt.xlabel(\"Weight\")\nplt.ylabel(\"Number of observations\")\nplt.title(\"Weight distribution\");","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d1dfe60c257f6a97f2eb4864e436af91d6721fe2","_cell_guid":"1f814eeb-b2dc-4fc6-9a15-d6b74cecec7e"},"cell_type":"markdown","source":"The weight properties are the following, according to the project description:\n - the few first (starting from the center of the detector) and last hits have a larger weight\n - hits from the more straight tracks (more rare, but more interesting) have a larger weight\n - random hits or hits from very short tracks have weight zero\n\nFrom these points and the above histogram, we can conclude that:\n- there are indeed very few straight tracks - with higher weights -, as explained on the description\n- there is a considerable amount of random hits or hits with very short tracks."},{"metadata":{"_uuid":"170fffe9552b8990757246625c1125a25ff5e8c7","_cell_guid":"0619cb47-fc9b-4e4b-ba8a-ae6911aef162"},"cell_type":"markdown","source":"### Particles dataset\n\nFirst, how many hits do particles of different charges have?"},{"metadata":{"_uuid":"196b7ac9c9ad9b82e57d292751c4b9f591b65356","_cell_guid":"a7dab669-5465-47d5-9b86-db19d99c88d6","trusted":true},"cell_type":"code","source":"g = sns.FacetGrid(data = particles, col = \"q\", hue = \"q\", size = 5)\ng.map(sns.distplot, \"nhits\", kde=False)\ng.add_legend();","execution_count":4,"outputs":[]},{"metadata":{"_uuid":"f09087996dfde830bf3064f820d9bdc520eceeaa"},"cell_type":"markdown","source":"Let's compute the absolute momentum of a particle and see if its somehow correlated with the number of hits."},{"metadata":{"trusted":true,"_uuid":"9cae1c33d585660e308bd3c0324fe35f002c60f9"},"cell_type":"code","source":"#computing absolute velocity\nparticles[\"abs_p\"] = np.sqrt(particles.pz**2 + particles.px**2 + particles.py**2)\n\n#making a scatterplot\nplt.scatter(particles.nhits, particles.abs_p)\nplt.title(\"Number of hits x Absolute momentum\")\nplt.xlabel(\"Number of hits\")\nplt.ylabel(\"Absolute momentum (GeV/c)\");\n\n#computing correlation\nparticles.abs_p.corr(particles.nhits)","execution_count":23,"outputs":[]},{"metadata":{"_uuid":"2b2fee17a669d248f2627ea41687050a5dd15178"},"cell_type":"markdown","source":"Both scatterplot and correlation coefficient indicate that there is no linear relationship between absolute momentum and the number of hits. "},{"metadata":{"trusted":true,"_uuid":"d2f57e2702b20c7f66f2d98541b67fc131dc5109"},"cell_type":"code","source":"#create figure instance\nfig = plt.figure(figsize=(25, 8))\ngs = gridspec.GridSpec(nrows=2, ncols=3, left=0.05, right=0.48, wspace=0.5, hspace = 0.3)\nax = fig.add_subplot(gs[0:2,0:2], projection = '3d')\nax2 =fig.add_subplot(gs[0,2])\nax3 =fig.add_subplot(gs[1,2])\n\n#plot the initial position for each particle on the information list\nax.scatter(\nxs=particles.vx,\nys=particles.vy,\nzs=particles.vz,\nalpha = 0.1)\n    \nax2.scatter( \nx = particles.vx,\ny = particles.vy,\nalpha = 0.1)\n    \n    \nax3.scatter(\nx= particles.vx,\ny = particles.vz,\nalpha = 0.1)\n    \n#labels\nax.set_xlabel('x (mm)')\nax.set_ylabel('y (mm)')\nax.set_zlabel('z (mm)')\nax.set_title('Initial position 3d view')\n\nax2.set_xlabel('x (mm)')\nax2.set_ylabel('y (mm)')\nax2.set_title('Initial position x-y cross section')\n\nax3.set_xlabel('x (mm)')\nax3.set_ylabel('z (mm)')\nax3.set_title('Initial position x-z cross section')\nplt.show()","execution_count":38,"outputs":[]},{"metadata":{"_uuid":"85442f4a661ebbeb7c8cd6fdf195371f24c3d38e"},"cell_type":"markdown","source":"Particles were generated very close to the origin of the detector, according to the global coordinates system, with some created in other parts of the apparatus."}],"metadata":{"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"},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"}},"nbformat":4,"nbformat_minor":1}