{"cells":[{"metadata":{"_cell_guid":"2e6866a2-fa92-46c9-a5bc-2785e1b17d2b","_uuid":"62bb7399cb3194a89d59ca77f6883d159f646eb1"},"cell_type":"markdown","source":"**Introduction**\n\nIn this kernel I'm investigating how the geometry is defined and trying to correctly calculate the cell coordinates.\nI'm trying to validate the calculated coordinates by restoring the hit coordinates from them and comparing those with the ones provided in the dataset."},{"metadata":{"_cell_guid":"61dc29c4-1206-41ed-84f8-8910e0433591","_uuid":"c19c88cb6bedd5538a496b629efac9cc7f6b7b8c"},"cell_type":"markdown","source":"Ok, so let's import everything we need..."},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","collapsed":true,"trusted":true},"cell_type":"code","source":"import gc\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom mpl_toolkits.mplot3d import Axes3D","execution_count":1,"outputs":[]},{"metadata":{"_cell_guid":"b00dbdc6-e558-49cd-b6b4-a490b53cc7cb","_uuid":"c5365f8e308ad1f50e297413ade859546a5cc8e1"},"cell_type":"markdown","source":"...and load some data:"},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","collapsed":true,"trusted":true},"cell_type":"code","source":"det_descr = pd.read_csv(\"../input/detectors.csv\")\ndet_descr.set_index(['volume_id', 'layer_id', 'module_id'], inplace=True)\n\n#hits  = pd.read_csv(\"../input/train_1/event000001160-hits.csv\")\n#cells = pd.read_csv(\"../input/train_1/event000001160-cells.csv\")\nhits  = pd.read_csv(\"../input/train_1/event000001100-hits.csv\")\ncells = pd.read_csv(\"../input/train_1/event000001100-cells.csv\")\n\n# Fetch module descriptions for each of the hits\nhit_descrs = det_descr.loc[[tuple(x) for x in hits[['volume_id', 'layer_id', 'module_id']].values]]\nhits_with_descr = pd.concat([hits, hit_descrs.set_index(hits.index)], axis=1)\nhits_with_descr.set_index('hit_id', inplace=True)\ndel hit_descrs; gc.collect()\n\n# Append hit and module info to each of the cells entries\nhits_aug = pd.concat([hits_with_descr.loc[cells.hit_id.values].set_index(cells.index), cells], axis=1)","execution_count":2,"outputs":[]},{"metadata":{"_cell_guid":"02655d11-d93c-4cb2-b44c-39e6d080d630","_uuid":"0640f53b43546e1208d164a2f86ac51370805de3"},"cell_type":"markdown","source":"Then, this bit of code below calculates cell center coordinates in the global frame, then averages them weighted by *value* to reproduce the hit coordinates and compares those with the provided hit coordinates."},{"metadata":{"_cell_guid":"2a75a490-c6a7-4333-a272-1e13bbb84754","_uuid":"3430a16e3dc79fbb8679779fd4a8b59309a6dc26","trusted":true},"cell_type":"code","source":"# Define columns for shorter formulas later\nx, y, z = hits_aug.x, hits_aug.y, hits_aug.z\n\ncell_iv = hits_aug.ch1\ncell_iu = hits_aug.ch0\ncx, cy, cz = hits_aug.cx, hits_aug.cy, hits_aug.cz\npitch_u, pitch_v = hits_aug.pitch_u, hits_aug.pitch_v\nmodule_hv = hits_aug.module_hv\nmodule_hu = hits_aug.module_maxhu\n\nrot_xu, rot_xv, rot_xw = hits_aug.rot_xu, hits_aug.rot_xv, hits_aug.rot_xw\nrot_yu, rot_yv, rot_yw = hits_aug.rot_yu, hits_aug.rot_yv, hits_aug.rot_yw\nrot_zu, rot_zv, rot_zw = hits_aug.rot_zu, hits_aug.rot_zv, hits_aug.rot_zw\n\n# Calculate the number of cells across each dimension\nnu = (module_hu * 2 / pitch_u).round()\nnv = (module_hv * 2 / pitch_v).round()\n\n# Some checks for the cell indexes\nassert (cell_iv >= 0).all()\nassert (cell_iv < nv).all()\nassert (cell_iu >= 0).all()\nassert (cell_iu < nu).all()\n\n# Calculating locall cell coordinates\nhit_v = -module_hv + (cell_iv + 0.5) * pitch_v\nhit_u = -module_hu + (cell_iu + 0.5) * pitch_u\n\n# Transforming to global (w = 0, i.e. we're interested in the depth center of the cell)\nhit_x = rot_xu * hit_u + rot_xv * hit_v + cx\nhit_y = rot_yu * hit_u + rot_yv * hit_v + cy\nhit_z = rot_zu * hit_u + rot_zv * hit_v + cz\n\n# Distance between the calculated cell coordinates and the provided hit coordinates\ndist_x = hit_x - x\ndist_y = hit_y - y\ndist_z = hit_z - z\n\ndist_df = pd.DataFrame({'dist_x' : dist_x         ,\n                        'dist_y' : dist_y         ,\n                        'dist_z' : dist_z         ,\n                        'value'  : hits_aug.value ,\n                        'hit_id' : hits_aug.hit_id})\ndist_df['dist_x_times_value'] = dist_df.dist_x * dist_df.value\ndist_df['dist_y_times_value'] = dist_df.dist_y * dist_df.value\ndist_df['dist_z_times_value'] = dist_df.dist_z * dist_df.value\ng = dist_df.groupby('hit_id')\n\nmean_dist = ((g.dist_x_times_value.sum() / g.value.sum())**2 + \\\n             (g.dist_y_times_value.sum() / g.value.sum())**2 + \\\n             (g.dist_z_times_value.sum() / g.value.sum())**2)**0.5\nplt.hist(mean_dist, log=True, bins=100);\nplt.show()","execution_count":3,"outputs":[]},{"metadata":{"_cell_guid":"512dfaa0-a8b1-4f18-988d-5f81a9a93a0a","_uuid":"3c5b08b9b94edc498dd977162afae7032d8c29f9"},"cell_type":"markdown","source":"Hmm, that looks weird... So we cannot reproduce the provided hit coordinates exactly (*any suggestions why is that so?*)...\n\nLets then transform global hit coordinates to local (u, v, w) and see where exactly that falls compared to our calculated average:"},{"metadata":{"_cell_guid":"28e7bb34-3ae8-4edb-84c0-76c1fa410ec9","_uuid":"14b97da56445d2fb4efdaa7c6cde5e0673aa34aa","trusted":true},"cell_type":"code","source":"# Let's also calculate local hit coordinates and compare those with the provided ones\n# These are the provided coordinates:\nu = (x - cx) * rot_xu + (y - cy) * rot_yu + (z - cz) * rot_zu\nv = (x - cx) * rot_xv + (y - cy) * rot_yv + (z - cz) * rot_zv\nw = (x - cx) * rot_xw + (y - cy) * rot_yw + (z - cz) * rot_zw\n\nlocal_df = pd.DataFrame({'u' : u,\n                         'v' : v,\n                         'w' : w,\n                         'cell_u_times_value' : hit_u * hits_aug.value,\n                         'cell_v_times_value' : hit_v * hits_aug.value,\n                         'value' : hits_aug.value,\n                         'hit_id' : hits_aug.hit_id,\n                         'tr' : (hits_aug.module_maxhu != hits_aug.module_minhu)})\n\n# And here are the coords calculated from the cell coords:\ng2 = local_df.groupby('hit_id')\n\nhit_u = g2.cell_u_times_value.sum() / g2.value.sum()\nhit_v = g2.cell_v_times_value.sum() / g2.value.sum()\n\n# Mean over same values to reduce the df\nu = g2.u.mean()\nv = g2.v.mean()\nw = g2.w.mean()\ntr = g2.tr.mean().astype(bool)\n\nfig = plt.figure(figsize=(18,5))\nax = fig.add_subplot(141, projection='3d')\nax.scatter((u - hit_u).loc[tr],\n           (v - hit_v).loc[tr],\n           w          .loc[tr], c='b', marker='o')\nax.scatter((u - hit_u).loc[~tr],\n           (v - hit_v).loc[~tr],\n           w          .loc[~tr], c='r', marker='^')\nax.set_xlabel('U')\nax.set_ylabel('V')\nax.set_zlabel('W')\n\nax = fig.add_subplot(142)\nax.scatter((u - hit_u).loc[tr],\n           (v - hit_v).loc[tr], c='b', marker='o')\nax.scatter((u - hit_u).loc[~tr],\n           (v - hit_v).loc[~tr], c='r', marker='^')\nax.set_xlabel('U')\nax.set_ylabel('V')\n\nax = fig.add_subplot(143)\nax.scatter((u - hit_u).loc[tr],\n           w          .loc[tr], c='b', marker='o')\nax.scatter((u - hit_u).loc[~tr],\n           w          .loc[~tr], c='r', marker='^')\nax.set_xlabel('U')\nax.set_ylabel('W')\n\nax = fig.add_subplot(144)\nax.scatter((v - hit_v).loc[tr],\n           w          .loc[tr], c='b', marker='o')\nax.scatter((v - hit_v).loc[~tr],\n           w          .loc[~tr], c='r', marker='^')\nax.set_xlabel('V')\nax.set_ylabel('W')\n\nplt.tight_layout()\nplt.show()","execution_count":4,"outputs":[]},{"metadata":{"_cell_guid":"b3dc0b15-8386-4b45-87c7-bb24ea1cff53","_uuid":"6fad2eef51b55fd33cfd314b6ebc9ebe45f00f73"},"cell_type":"markdown","source":"Below are some other checks I did earlier."},{"metadata":{"_cell_guid":"a7dcfdec-5d50-4983-b8ce-567b6ea8234f","_uuid":"46ff8bcd004c15252b83161a6824150db0e78a81","trusted":false,"collapsed":true},"cell_type":"code","source":"# let's check the rotations are orthogonal:\nrotations = [[hits_aug.rot_xu, hits_aug.rot_xv, hits_aug.rot_xw],\n             [hits_aug.rot_yu, hits_aug.rot_yv, hits_aug.rot_yw],\n             [hits_aug.rot_zu, hits_aug.rot_zv, hits_aug.rot_zw],]\n\nsum2_x = rotations[0][0]**2 + rotations[0][1]**2 + rotations[0][2]**2\nsum2_y = rotations[1][0]**2 + rotations[1][1]**2 + rotations[1][2]**2\nsum2_z = rotations[2][0]**2 + rotations[2][1]**2 + rotations[2][2]**2\n\nsum_xy = rotations[0][0]*rotations[1][0] + rotations[0][1]*rotations[1][1] + rotations[0][2]*rotations[1][2]\nsum_yz = rotations[1][0]*rotations[2][0] + rotations[1][1]*rotations[2][1] + rotations[1][2]*rotations[2][2]\nsum_zx = rotations[2][0]*rotations[0][0] + rotations[2][1]*rotations[0][1] + rotations[2][2]*rotations[0][2]\n\nplt.figure(figsize=(18,8))\nplt.subplot(231)\nplt.hist(sum2_x, log=True, bins=100);\nplt.subplot(232)\nplt.hist(sum2_y, log=True, bins=100);\nplt.subplot(233)\nplt.hist(sum2_z, log=True, bins=100);\nplt.subplot(234)\nplt.hist(sum_xy, log=True, bins=100);\nplt.subplot(235)\nplt.hist(sum_yz, log=True, bins=100);\nplt.subplot(236)\nplt.hist(sum_zx, log=True, bins=100);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"24d4bca2-0505-400d-bb3b-edb814054def","_uuid":"752e1c22eaba6695210c90fba7ac03611668a5f6","trusted":false,"collapsed":true},"cell_type":"code","source":"# Norm to module plane in local coords is u = 0, v = 0, w = 1.\n# Hence in global coordinates it is (rot_xw, rot_yw, rot_zw)\n# Let's plot its dot product with two vectors: (cx, cy, 0) and (0, 0, cz):\n\nnorm_xy = (det_descr.cx**2 + det_descr.cy**2)**0.5\nsel1 = det_descr.module_minhu != det_descr.module_maxhu\nsel2 = det_descr.module_minhu == det_descr.module_maxhu\n\nplt.figure(figsize=(15,5))\nplt.subplot(121)\nplt.scatter(((det_descr.rot_xw * det_descr.cx + det_descr.rot_yw * det_descr.cy) / norm_xy).loc[sel1],\n             (det_descr.rot_zw * np.sign(det_descr.cz)).loc[sel1]);\nplt.subplot(122)\nplt.scatter(((det_descr.rot_xw * det_descr.cx + det_descr.rot_yw * det_descr.cy) / norm_xy).loc[sel2],\n             (det_descr.rot_zw * np.sign(det_descr.cz)).loc[sel2]);","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b79ad6ee-44d8-40a2-91c7-c4fab6f9892e","_uuid":"9a59f0c274c6be39d75002ed5822ecf73888725c","collapsed":true,"trusted":false},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"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}