{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport plotly.express as px\nimport plotly.graph_objects as go\n\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-02-12T14:45:16.380410Z","iopub.execute_input":"2023-02-12T14:45:16.380843Z","iopub.status.idle":"2023-02-12T14:45:17.581489Z","shell.execute_reply.started":"2023-02-12T14:45:16.380746Z","shell.execute_reply":"2023-02-12T14:45:17.580248Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Loading and checking the data","metadata":{}},{"cell_type":"code","source":"PATH_DATASETS = '/kaggle/input/icecube-neutrinos-in-deep-ice/'","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:17.583134Z","iopub.execute_input":"2023-02-12T14:45:17.583505Z","iopub.status.idle":"2023-02-12T14:45:17.589506Z","shell.execute_reply.started":"2023-02-12T14:45:17.583471Z","shell.execute_reply":"2023-02-12T14:45:17.587839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta = pd.read_parquet(os.path.join(PATH_DATASETS,'train_meta.parquet'))","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:17.591824Z","iopub.execute_input":"2023-02-12T14:45:17.592335Z","iopub.status.idle":"2023-02-12T14:45:29.763725Z","shell.execute_reply.started":"2023-02-12T14:45:17.592294Z","shell.execute_reply":"2023-02-12T14:45:29.762077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_batch_125 = train_meta[train_meta.batch_id == 125]\nprint(f\"Shape: {train_meta_batch_125.shape}\")\ntrain_meta_batch_125.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:29.766949Z","iopub.execute_input":"2023-02-12T14:45:29.767386Z","iopub.status.idle":"2023-02-12T14:45:30.296622Z","shell.execute_reply.started":"2023-02-12T14:45:29.767349Z","shell.execute_reply":"2023-02-12T14:45:30.295547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_geometry = pd.read_csv(os.path.join(PATH_DATASETS, 'sensor_geometry.csv'))\nprint(f\"Shape: {sensor_geometry.shape}\")\nsensor_geometry.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:30.298039Z","iopub.execute_input":"2023-02-12T14:45:30.298962Z","iopub.status.idle":"2023-02-12T14:45:30.322712Z","shell.execute_reply.started":"2023-02-12T14:45:30.298919Z","shell.execute_reply":"2023-02-12T14:45:30.321685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_batch_125 = pd.read_parquet( \\\n    os.path.join(PATH_DATASETS, 'train/batch_125.parquet')).reset_index()\nprint(f\"Shape: {train_batch_125.shape}\")\nprint(f\"The number of unique events: {train_batch_125.event_id.nunique()}\")\ntrain_batch_125.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:30.324113Z","iopub.execute_input":"2023-02-12T14:45:30.324717Z","iopub.status.idle":"2023-02-12T14:45:32.089291Z","shell.execute_reply.started":"2023-02-12T14:45:30.324677Z","shell.execute_reply":"2023-02-12T14:45:32.087818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Adding and exploring some new features for events","metadata":{}},{"cell_type":"code","source":"train_meta_batch_125['num_of_pulses'] = train_meta_batch_125['last_pulse_index'] - \\\n                train_meta_batch_125['first_pulse_index'] + 1","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:32.091233Z","iopub.execute_input":"2023-02-12T14:45:32.092157Z","iopub.status.idle":"2023-02-12T14:45:32.103643Z","shell.execute_reply.started":"2023-02-12T14:45:32.092101Z","shell.execute_reply":"2023-02-12T14:45:32.101500Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(8, 6))\nsns.distplot(train_meta_batch_125.num_of_pulses)\nplt.xlabel('Number of pulses per event', size=11)\nplt.ylabel('Frequency', size=11)\nplt.title('Distribution of the number of pulses', size=17)\nplt.show()\n\nprint(f'The mean number of pulses per event: \\\n{train_meta_batch_125.num_of_pulses.mean():.1f} pulses')\nprint(f'The median number of pulses per event: \\\n{train_meta_batch_125.num_of_pulses.median():.1f} pulses')\nprint(f'The number of pulses of 95% of all events in this batch is less than \\\n{train_meta_batch_125.num_of_pulses.quantile(0.95):.1f} pulses')","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:32.106301Z","iopub.execute_input":"2023-02-12T14:45:32.106936Z","iopub.status.idle":"2023-02-12T14:45:33.257214Z","shell.execute_reply.started":"2023-02-12T14:45:32.106867Z","shell.execute_reply":"2023-02-12T14:45:33.255596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"total_event_sensors = train_batch_125.groupby('event_id').sensor_id.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:33.259193Z","iopub.execute_input":"2023-02-12T14:45:33.259765Z","iopub.status.idle":"2023-02-12T14:45:46.772963Z","shell.execute_reply.started":"2023-02-12T14:45:33.259713Z","shell.execute_reply":"2023-02-12T14:45:46.771697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_batch_125 = train_meta_batch_125 \\\n                .merge(total_event_sensors, on='event_id') \\\n                .rename(columns={'sensor_id': 'num_of_unique_sensors'})","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:46.774872Z","iopub.execute_input":"2023-02-12T14:45:46.775774Z","iopub.status.idle":"2023-02-12T14:45:46.872101Z","shell.execute_reply.started":"2023-02-12T14:45:46.775717Z","shell.execute_reply":"2023-02-12T14:45:46.870929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(8, 6))\nsns.distplot(train_meta_batch_125.num_of_unique_sensors)\nplt.xlabel('Number of unique sensors', size=11)\nplt.ylabel('Frequency', size=11)\nplt.title('Distribution of the number of unique sensors', size=17)\nplt.show()\n\nprint(f'The average number of unique sensors per event: \\\n{train_meta_batch_125.num_of_unique_sensors.mean():.1f}')\nprint(f'95% of all events in this batch have less than \\\n{train_meta_batch_125.num_of_unique_sensors.quantile(0.95)} unique sensors')","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:46.873495Z","iopub.execute_input":"2023-02-12T14:45:46.874015Z","iopub.status.idle":"2023-02-12T14:45:47.929865Z","shell.execute_reply.started":"2023-02-12T14:45:46.873945Z","shell.execute_reply":"2023-02-12T14:45:47.928546Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sum_of_charges_per_event = train_batch_125 \\\n                    .groupby('event_id').charge.sum().round(2)","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:47.931269Z","iopub.execute_input":"2023-02-12T14:45:47.931663Z","iopub.status.idle":"2023-02-12T14:45:48.714333Z","shell.execute_reply.started":"2023-02-12T14:45:47.931627Z","shell.execute_reply":"2023-02-12T14:45:48.713254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_meta_batch_125 = train_meta_batch_125 \\\n                    .merge(sum_of_charges_per_event, on='event_id') \\\n                    .rename(columns={'charge': 'sum_of_charges'})","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:48.720552Z","iopub.execute_input":"2023-02-12T14:45:48.720964Z","iopub.status.idle":"2023-02-12T14:45:48.823607Z","shell.execute_reply.started":"2023-02-12T14:45:48.720928Z","shell.execute_reply":"2023-02-12T14:45:48.822647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(8, 6))\nsns.distplot(train_meta_batch_125.sum_of_charges)\nplt.xlabel('Total charge of light produced per event', size=11)\nplt.ylabel('Frequency', size=11)\nplt.title('Distribution of the total charge of light produced per event', size=12)\nplt.show()\n\nprint(f'The mean total charge of light produced per event: \\\n{train_meta_batch_125.sum_of_charges.mean():.1f} p.e.')\nprint(f'The median total charge of light produced per event: \\\n{train_meta_batch_125.sum_of_charges.median():.1f} p.e.')\nprint(f'The total light charge of 95% of all events in this batch is less than \\\n{train_meta_batch_125.sum_of_charges.quantile(0.95):.1f} p.e.')","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:48.824809Z","iopub.execute_input":"2023-02-12T14:45:48.826110Z","iopub.status.idle":"2023-02-12T14:45:49.898572Z","shell.execute_reply.started":"2023-02-12T14:45:48.826068Z","shell.execute_reply":"2023-02-12T14:45:49.897144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(8, 6))\n\nsns.scatterplot(train_meta_batch_125.num_of_pulses, \\\n                train_meta_batch_125.sum_of_charges, alpha=0.5)\nplt.xlabel('Number of pulses', size=11)\nplt.ylabel('Total charge of light', size=11)\nplt.title('Scatterplot of the number of pulses vs the total charge of light', \\\n          size=12)\nplt.show()\n\ncorr_coef = train_meta_batch_125[['num_of_pulses', 'sum_of_charges']] \\\n                .corr(method='spearman').values[0][1]\nprint(f'Spearman’s correlation coefficient : {corr_coef:.2f}')","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:49.900603Z","iopub.execute_input":"2023-02-12T14:45:49.901183Z","iopub.status.idle":"2023-02-12T14:45:50.638568Z","shell.execute_reply.started":"2023-02-12T14:45:49.901122Z","shell.execute_reply":"2023-02-12T14:45:50.636964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Spearman’s correlation coefficient shows that there is a very strong monotonic relationship between this two variables. This means the less light produced per event, the fewer number of sensors hit and record pulses.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(8, 6))\n\nsns.scatterplot(train_meta_batch_125.num_of_unique_sensors, \\\n                train_meta_batch_125.sum_of_charges, alpha=0.5)\nplt.xlabel('Number of unique sensors', size=11)\nplt.ylabel('Total charge of light', size=11)\nplt.title('Scatterplot of the number of unique sensors vs the total charge of light', \\\n          size=11)\nplt.show()\n\ncorr_coef = train_meta_batch_125[['num_of_unique_sensors', 'sum_of_charges']] \\\n                .corr(method='spearman').values[0][1]\nprint(f'Spearman’s correlation coefficient : {corr_coef:.2f}')","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:50.640711Z","iopub.execute_input":"2023-02-12T14:45:50.641334Z","iopub.status.idle":"2023-02-12T14:45:51.204866Z","shell.execute_reply.started":"2023-02-12T14:45:50.641262Z","shell.execute_reply":"2023-02-12T14:45:51.203252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Spearman’s correlation coefficient shows that there is a strong monotonic relationship between this two variables. This means the less light produced per event, the fewer number of unique sensors record pulses.","metadata":{}},{"cell_type":"markdown","source":"Thus, on average events have around 63 pulses that have 60 p.e. of total charge of light and that are recorded by 61 unique sensors. \\\n But it should be noted that distributions of pulses, unique sensors and total charge of light have very long tails on the right side of them. It means there are very \"bright\" events: with a huge number of pulses that have a large total charge of light and which were detected by a large number of unique sensors. \\\n Let's take a look at one of these \"bright\" events.","metadata":{}},{"cell_type":"markdown","source":"### Visualizing one of the \"brightest\" events","metadata":{}},{"cell_type":"markdown","source":"Let's visualize the event with the largest number of unique triggered sensors.","metadata":{}},{"cell_type":"code","source":"event_id_1 =  train_meta_batch_125.loc[train_meta_batch_125.num_of_unique_sensors == \\\n                    train_meta_batch_125.num_of_unique_sensors.max(), 'event_id'].values[0]","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:51.207215Z","iopub.execute_input":"2023-02-12T14:45:51.207659Z","iopub.status.idle":"2023-02-12T14:45:51.216765Z","shell.execute_reply.started":"2023-02-12T14:45:51.207619Z","shell.execute_reply":"2023-02-12T14:45:51.215058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses_1 = train_batch_125[train_batch_125.event_id == event_id_1].copy()\nprint(event_pulses_1.shape)\nevent_pulses_1.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:51.218826Z","iopub.execute_input":"2023-02-12T14:45:51.219250Z","iopub.status.idle":"2023-02-12T14:45:51.812924Z","shell.execute_reply.started":"2023-02-12T14:45:51.219210Z","shell.execute_reply":"2023-02-12T14:45:51.811362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses_1 = event_pulses_1.reset_index()\nevent_pulses_1[['x', 'y', 'z']] = sensor_geometry.loc[event_pulses_1.sensor_id] \\\n             .reset_index()[['x', 'y', 'z']]\n\nevent_pulses_1.x = event_pulses_1.x - event_pulses_1.x.mean()\nevent_pulses_1.y = event_pulses_1.y - event_pulses_1.y.mean()\nevent_pulses_1.z = event_pulses_1.z - event_pulses_1.z.mean()\nevent_pulses_1.head()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:51.814165Z","iopub.execute_input":"2023-02-12T14:45:51.814561Z","iopub.status.idle":"2023-02-12T14:45:51.856557Z","shell.execute_reply.started":"2023-02-12T14:45:51.814526Z","shell.execute_reply":"2023-02-12T14:45:51.855474Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"zenith = train_meta_batch_125.loc[ \\\n    train_meta_batch_125.event_id == event_id_1, 'zenith'].values[0]\nazimuth = train_meta_batch_125.loc[ \\\n    train_meta_batch_125.event_id == event_id_1, 'azimuth'].values[0]\n\nvector = [np.cos(azimuth) * np.sin(zenith), \\\n          np.sin(azimuth) * np.sin(zenith), \\\n          np.cos(zenith)]\n\nvector_base = np.array([-500, 500])\n\nx = vector_base * vector[0]\ny = vector_base * vector[1]\nz = vector_base * vector[2]\n\nvector_df = pd.DataFrame({'x': x, 'y': y, 'z': z})\nvector_df","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:51.857963Z","iopub.execute_input":"2023-02-12T14:45:51.858568Z","iopub.status.idle":"2023-02-12T14:45:51.878601Z","shell.execute_reply.started":"2023-02-12T14:45:51.858531Z","shell.execute_reply":"2023-02-12T14:45:51.877334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_3d(event_pulses_1, x='x', y='y', z='z', opacity=0.5, color='time')\nfig.update_traces(marker_size=2)\nfig_line = px.line_3d(vector_df, x='x', y='y', z='z')\nfig_line.update_traces(line=dict(color=\"Black\", width=5.5))\nfig = go.Figure(data = fig.data + fig_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:51.880575Z","iopub.execute_input":"2023-02-12T14:45:51.881015Z","iopub.status.idle":"2023-02-12T14:45:52.405058Z","shell.execute_reply.started":"2023-02-12T14:45:51.880951Z","shell.execute_reply":"2023-02-12T14:45:52.399638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Because of the fact that so many Digital Optical Modules (DOMs) have detected the Cherenkov light it has become more difficult to identify the true diraction of the neutrino. Let's try to filter those triggered sensors by finding only the most active sensors in the event.","metadata":{}},{"cell_type":"markdown","source":"### Finding the most active sensors","metadata":{}},{"cell_type":"markdown","source":"Let's calculate the total light charge and the number of pulses that each sensor registers in the event. This will allow us to prioritize sensors and highlight the most active sensors for those events that have a very large number of unique triggered sensors.","metadata":{}},{"cell_type":"code","source":"sensors_info = train_batch_125.groupby(['event_id', 'sensor_id'], as_index=False) \\\n                .agg(charge_per_sensor=('charge', 'sum'), \\\n                     n_pulses_per_sensor=('charge', 'count'))","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:45:52.406616Z","iopub.execute_input":"2023-02-12T14:45:52.407258Z","iopub.status.idle":"2023-02-12T14:46:03.132272Z","shell.execute_reply.started":"2023-02-12T14:45:52.407155Z","shell.execute_reply":"2023-02-12T14:46:03.131213Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensors_info.head()\n","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:03.135546Z","iopub.execute_input":"2023-02-12T14:46:03.136176Z","iopub.status.idle":"2023-02-12T14:46:03.150420Z","shell.execute_reply.started":"2023-02-12T14:46:03.136117Z","shell.execute_reply":"2023-02-12T14:46:03.149123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensors_info[['charge_per_sensor', 'n_pulses_per_sensor']].describe().round(2)","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:03.152391Z","iopub.execute_input":"2023-02-12T14:46:03.154020Z","iopub.status.idle":"2023-02-12T14:46:04.611858Z","shell.execute_reply.started":"2023-02-12T14:46:03.153946Z","shell.execute_reply":"2023-02-12T14:46:04.610409Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's find the most active sensors for our event.","metadata":{}},{"cell_type":"code","source":"sensors_info_per_event = sensors_info[sensors_info.event_id == event_id_1]","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:04.614857Z","iopub.execute_input":"2023-02-12T14:46:04.615398Z","iopub.status.idle":"2023-02-12T14:46:04.664989Z","shell.execute_reply.started":"2023-02-12T14:46:04.615356Z","shell.execute_reply":"2023-02-12T14:46:04.663507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sensor_pulses_threshold = sensors_info_per_event.n_pulses_per_sensor.quantile(0.7)\nsensor_charge_threshold = sensors_info_per_event.charge_per_sensor.quantile(0.7)","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:04.666670Z","iopub.execute_input":"2023-02-12T14:46:04.667168Z","iopub.status.idle":"2023-02-12T14:46:04.679117Z","shell.execute_reply.started":"2023-02-12T14:46:04.667125Z","shell.execute_reply":"2023-02-12T14:46:04.677340Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"more_active_sensors = sensors_info_per_event[ \\\n        (sensors_info_per_event.n_pulses_per_sensor > sensor_pulses_threshold) & \\\n        (sensors_info_per_event.charge_per_sensor > sensor_charge_threshold)] \\\n        .sensor_id.to_list()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:04.682246Z","iopub.execute_input":"2023-02-12T14:46:04.683277Z","iopub.status.idle":"2023-02-12T14:46:04.697131Z","shell.execute_reply.started":"2023-02-12T14:46:04.683215Z","shell.execute_reply":"2023-02-12T14:46:04.694737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event_pulses_relevant = event_pulses_1[event_pulses_1.sensor_id.isin(more_active_sensors)]\nevent_pulses_relevant.shape","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:04.699857Z","iopub.execute_input":"2023-02-12T14:46:04.700634Z","iopub.status.idle":"2023-02-12T14:46:04.728221Z","shell.execute_reply.started":"2023-02-12T14:46:04.700573Z","shell.execute_reply":"2023-02-12T14:46:04.727116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_3d(event_pulses_relevant, x='x', y='y', z='z', opacity=0.5, color='time')\nfig.update_traces(marker_size=2)\nfig_line = px.line_3d(vector_df, x=\"x\", y=\"y\", z=\"z\")\nfig_line.update_traces(line=dict(color=\"Black\", width=5.5))\nfig = go.Figure(data = fig.data + fig_line.data)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-02-12T14:46:04.729480Z","iopub.execute_input":"2023-02-12T14:46:04.730716Z","iopub.status.idle":"2023-02-12T14:46:04.956178Z","shell.execute_reply.started":"2023-02-12T14:46:04.730662Z","shell.execute_reply":"2023-02-12T14:46:04.953568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks much better. Thus, for those events that have a very large number of pulses and, therefore, the number of unique sensors that recorded these pulses we can filter their sensors by finding the most active sensors: with more charge of light and more number of pulses recorded by them.","metadata":{}}]}