{"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":"markdown","source":"![header.jpg](https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/header.jpg)\n## Goal of the Competition\nThe goal of this competition is to predict a neutrino particle’s direction. *Participants* will develop models based on data from the \"IceCube\" detector, which observes the cosmos from deep within the South Pole ice.<br>\n*This* work could help scientists better understand exploding stars, gamma-ray bursts, and cataclysmic phenomena involving black holes, neutron stars and the fundamental properties of the neutrino itself.","metadata":{}},{"cell_type":"markdown","source":"## Sections\n* [The Observatory in brief](#brief_intro)\n    * [Observatory Architecture](#s1ss1)\n    * [Neutrino Events](#s1ss2)\n    * [Sources of Noise](#s1ss3)\n    * [A catch in the assumption about Ice](#s1ss4)\n* [Data Overview](#data_overview)\n* [Visualizing Events](#event_visualization)\n    * [Some cool events](#event_visualization)\n    * [Event Animation](#event_animation)","metadata":{}},{"cell_type":"markdown","source":"## The Observatory in Brief <a class=\"anchor\"  id=\"brief_intro\"></a>","metadata":{}},{"cell_type":"markdown","source":"### Observatory Architecture <a class=\"anchor\"  id=\"s1ss1\"></a>\nIceCube, the South Pole neutrino observatory, is a cubic-kilometer particle detector made of Antarctic ice and located near the Amundsen-Scott South Pole Station. It is buried beneath the surface, extending to a depth of about 2,500 meters. A surface array, IceTop, and a denser inner subdetector, DeepCore, significantly enhance the capabilities of the observatory, making it a multipurpose facility.<br>\nThe **`in-ice component of IceCube consists of 5,160 digital optical modules (DOMs)`**, each with a ten-inch photomultiplier tube and associated electronics. The **`DOMs are attached to vertical “strings,” frozen into 86 boreholes`**, and arrayed over a cubic kilometer from 1,450 meters to 2,450 meters depth. The **`strings are deployed on a hexagonal grid with 125 meters spacing`** and hold 60 DOMs each. The **`vertical separation of the DOMs is 17 meters`**.<br>\nIceTop consists of 81 stations located on top of the same number of IceCube strings. Each station has two tanks, each equipped with **`two downward facing DOMs. IceTop, built as a veto`** and calibration detector for IceCube, also detects air showers from primary cosmic rays in the 300 TeV to 1 EeV energy range. The **`surface array measures the cosmic-ray arrival directions in the Southern Hemisphere`** as well as the flux and composition of cosmic rays.<br>\n![Icecube-architecture-diagram2009.png](https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/Icecube-architecture-diagram2009.png)","metadata":{"_kg_hide-input":false}},{"cell_type":"markdown","source":"### Neutrino Events <a class=\"anchor\"  id=\"s1ss2\"></a>\nThere are primarity three types of signatures depending on the flavour of the neutrino that had some interaction.\n* The interaction may occur outside the detector especially in case of muon track type detections.\n* Electron neutrino events decay in energy over short distances so are more likely to be confined within the detector.","metadata":{}},{"cell_type":"markdown","source":"<img src=\"https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/signatures.PNG\" width=\"600\"/><br>\nRef: <a href=\"https://indico.cern.ch/event/472838/contributions/1150248/attachments/1296103/1932596/CAP_2016_IceCube.pdf\"> IceCube-DeepCore-PINGU, CAP Congress 2016, Darren R Grant </a>","metadata":{}},{"cell_type":"markdown","source":"### Sources of Noise <a class=\"anchor\"  id=\"s1ss3\"></a>\n\n> *There is a **large background of muons created** not by neutrinos from astrophysical sources but **by cosmic rays impacting the atmosphere above the detector**. There are about 106 times more cosmic ray muons than neutrino-induced muons observed in IceCube. Most of these can be rejected using the fact that they are traveling downwards. Most of the remaining (up-going) events are from neutrinos, but most of these neutrinos are from cosmic rays hitting the far side of the Earth; some unknown fraction may come from astronomical sources, and these neutrinos are the key to IceCube point source searches.* <a href=\"https://en.wikipedia.org/wiki/IceCube_Neutrino_Observatory\">- Wikipedia</a>","metadata":{}},{"cell_type":"markdown","source":"<div style=\"align:center; margin:auto\">\n    <div style=\"width:400px;\">\n        <img src=\"https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/event_downgoing_neutrino.PNG\" \n             width=\"400px\" \n             style=\"float: left\" \n             />\n    </div>\n    <div style=\"width:100px; color:#FFFFFF; float: left;\">words</div>\n    <div style=\"width:800px; vertical-align: middle; display: table-cell; min-height: 10em;\"><br><br><br><br><br>\n        The IceCube Observatory is most sensitive to muon track events. Atmospheric <strong>\"downgoing\" muon tracks such as the one on the left are hence most likely to be noise</strong> in the data. \n        Necessary filtering mechanisms have to be implemented for the model.\n    </div>\n</div>","metadata":{}},{"cell_type":"markdown","source":"### A catch in the assumption about Ice <a class=\"anchor\"  id=\"s1ss4\"></a>","metadata":{}},{"cell_type":"markdown","source":"<div style=\"align:center; margin:auto\">\n    <div style=\"width:700px;\">\n        <img src=\"https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/ice_properties.PNG\" \n             width=\"600px\" \n             style=\"float: left\" \n             />\n    </div>\n    <div style=\"width:100px; color:#FFFFFF; float: left;\">words</div>\n    <div style=\"width:800px; vertical-align: middle; display: table-cell; min-height: 10em;\"><br><br><br>\n        While the South Pole ice at the depth of instrumentation is the clearest ice available, <strong>it is not uniform in its properties. There are wide variations in the optical properties of the ice.</strong> This may be taken into account when creating the model.\n    </div>\n</div>","metadata":{}},{"cell_type":"markdown","source":"## Data Overview <a class=\"anchor\"  id=\"data_overview\"></a>","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport plotly.graph_objects as go\nimport plotly.offline as po\nfrom plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\ninit_notebook_mode(connected=True)\nimport os","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:21:54.458840Z","iopub.execute_input":"2023-02-05T07:21:54.459304Z","iopub.status.idle":"2023-02-05T07:21:54.553920Z","shell.execute_reply.started":"2023-02-05T07:21:54.459211Z","shell.execute_reply":"2023-02-05T07:21:54.552812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"PATH = '/kaggle/input/icecube-neutrinos-in-deep-ice/'\ndoms = pd.read_csv(f'{PATH}sensor_geometry.csv')\ndata_files = os.listdir(f'{PATH}train')\ntrain_meta = pd.read_parquet(f'{PATH}train_meta.parquet')\ntrain_events = train_meta.event_id.to_list()","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:21:54.555928Z","iopub.execute_input":"2023-02-05T07:21:54.556280Z","iopub.status.idle":"2023-02-05T07:23:17.365348Z","shell.execute_reply.started":"2023-02-05T07:21:54.556248Z","shell.execute_reply":"2023-02-05T07:23:17.362677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"event = 244748\nevent_meta = train_meta[train_meta.event_id == event].iloc[0].to_dict()\ncurrent_event = pd.read_parquet(f'{PATH}train/batch_'+str(train_meta[train_meta.event_id == event].batch_id.iloc[0])+'.parquet')\ncurrent_event = current_event[current_event.index == event].reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:17.367997Z","iopub.execute_input":"2023-02-05T07:23:17.368601Z","iopub.status.idle":"2023-02-05T07:23:24.695926Z","shell.execute_reply.started":"2023-02-05T07:23:17.368527Z","shell.execute_reply":"2023-02-05T07:23:24.694838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(event_meta, index = [0])","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:24.698464Z","iopub.execute_input":"2023-02-05T07:23:24.699042Z","iopub.status.idle":"2023-02-05T07:23:29.772778Z","shell.execute_reply.started":"2023-02-05T07:23:24.699002Z","shell.execute_reply":"2023-02-05T07:23:29.770849Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The meta data of an event has the ***target variables:*** `azimuth` and `zenith` angles. Together **the two angles are defined in a coordinate system with origin at** roughly the **center of the detector**. They **specify a direction vector** for observers **to point in to identify the source** of the neutrino.  <font color='red'><b>Since the target is a direction vector, it can be expected that in most cases the actual position of the neutrino events will be in a different position in the detector.</b></font>","metadata":{}},{"cell_type":"markdown","source":"<div style=\"align:center; margin:auto\">\n    <img src=\"https://raw.githubusercontent.com/Chahbaz-Aman/datastore/main/IceCube%20-%20Neutrinos%20in%20deep%20ice/az_zen.PNG\" width=\"500\"/><br>\n    Ref: <a href=\"https://docs.icecube.aq/icetray/main/projects/dataclasses/coordinates.html\"> IceCube Docs »icetray (8f4d09b9) »Dataclasses documentation »Coordinate Systems </a>\n</div>","metadata":{}},{"cell_type":"code","source":"current_event.iloc[12:12+10]","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:29.777962Z","iopub.execute_input":"2023-02-05T07:23:29.778325Z","iopub.status.idle":"2023-02-05T07:23:29.794662Z","shell.execute_reply.started":"2023-02-05T07:23:29.778292Z","shell.execute_reply":"2023-02-05T07:23:29.793317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The event data records each sensor activation, the relative time instant & the intensity of the light (name: charge). The field `auxiliary` is set to true/false by the detector software. The idea is if a DOM detects photons, it checks if any neighboring DOMs are also detecting the same event. If yes, auxiliary is set to False indicating the detection is most likely not random noise. However, the `auxiliary=False` may be due to a downward going muon which would still be noise for our purpose. Auxiliary detections are highly likely to be carrying valuable information.","metadata":{}},{"cell_type":"code","source":"del event, event_meta, current_event","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:29.796092Z","iopub.execute_input":"2023-02-05T07:23:29.796400Z","iopub.status.idle":"2023-02-05T07:23:29.806731Z","shell.execute_reply.started":"2023-02-05T07:23:29.796371Z","shell.execute_reply":"2023-02-05T07:23:29.805711Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Visualizing Events <a class=\"anchor\"  id=\"event_visualization\"></a>","metadata":{}},{"cell_type":"markdown","source":"### Some Cool Events","metadata":{}},{"cell_type":"code","source":"cool_events = [\n    1128590,    # Finder's Credit: Peter https://www.kaggle.com/competitions/icecube-neutrinos-in-deep-ice/discussion/381935\n    244748,     # muon track starting from the top end of the observatory - this is noise labeled as non-auxiliary\n    1567684848, # downward going muon - this is noise labeled as non-auxiliary\n    1636247611, # possible electron neutrino / neutral current interaction\n] ","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:29.808907Z","iopub.execute_input":"2023-02-05T07:23:29.810018Z","iopub.status.idle":"2023-02-05T07:23:29.820895Z","shell.execute_reply.started":"2023-02-05T07:23:29.809974Z","shell.execute_reply":"2023-02-05T07:23:29.819724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def view_event_by_id(event_id = None):\n\n    PICK_RANDOM_EVENT = event_id == None\n    AUX_LABEL = 'auxiliary = True'\n    NON_AUX_LABEL = 'auxiliary = False'\n\n    if PICK_RANDOM_EVENT:\n        print('Picking random event...')\n        event = np.random.choice(train_meta.event_id) # random event for preliminary EDA\n    else:\n        event = event_id\n\n    print('1. Gathering event data...')\n    event_meta = train_meta[train_meta.event_id == event].iloc[0].to_dict()\n    current_event = pd.read_parquet(f'{PATH}train/batch_'+str(train_meta[train_meta.event_id == event].batch_id.iloc[0])+'.parquet')\n    current_event = current_event[current_event.index == event]\n    current_event = pd.merge(current_event, doms, on = 'sensor_id').sort_values(by = 'time').reset_index()\n    non_aux = current_event[current_event[\"auxiliary\"] == False].reset_index()\n\n    ANIMATION_FRAME_DURATION = int(5000/current_event.shape[0]+1)\n\n    # add DOMs\n    print('2. Generating plotly animation...')\n    dom_fig = go.Scatter3d(x = doms[\"x\"], y = doms[\"y\"], z = doms[\"z\"], \n                        opacity = 0.3,\n                        marker = dict(size = 2, \n                                      color = 'lightblue', \n                                     ),\n                        mode = 'markers',\n                        showlegend = False\n                       )\n\n    # add path\n    lx = np.cos(event_meta['azimuth'])*np.sin(event_meta['zenith'])\n    ly = np.sin(event_meta['azimuth'])*np.sin(event_meta['zenith'])\n    lz = np.cos(event_meta['zenith'])\n    path_fig = go.Scatter3d(name = 'neutrino direction',\n                            showlegend = False,\n                            x = [-500*lx, 500*lx],\n                            y = [-500*ly, 500*ly],\n                            z = [-500*lz, 500*lz],\n                            mode = 'lines',\n                            line = dict(width = 5, \n                                        color = 'red'\n                                       ),\n                           )\n\n    x,y,z = current_event['x'], current_event['y'], current_event['z']\n    xna,yna,zna = non_aux['x'], non_aux['y'], non_aux['z']\n\n\n    # Create figure\n    size_delta = current_event['charge'].max() - current_event['charge'].min()\n\n    fig = go.Figure(data=[go.Scatter3d(name = AUX_LABEL,\n                                       x=[], y=[], z=[],\n                                       mode=\"markers\",\n                                       marker=dict(color = current_event['time'], \n                                                   size  = [(int(v)*20/size_delta + 10) for v in current_event['charge'].to_list()],\n                                                   colorscale = 'aggrnyl', \n                                                  ),\n                                      ),\n\n                          go.Scatter3d(name = 'auxiliary = False',\n                                       x=[], y=[], z=[],\n                                       mode=\"markers\",\n                                       marker=dict(color = non_aux['time'], \n                                                   size  = [(int(v)*20/size_delta + 10) for v in non_aux['charge'].to_list()],\n                                                   line = dict(width = 3, color = 'red'),\n                                                   colorscale = 'agsunset',\n                                                  ),\n                                      ),\n                         ]\n                   )\n\n    frames = []\n    i = 0\n    d = int(current_event.shape[0]//100) + 1\n    for k in range(0,current_event.shape[0]-1,d):\n        time = current_event['time'][k]\n\n        data = [go.Scatter3d(x=x[:k+1], \n                             y=y[:k+1], \n                             z=z[:k+1],\n                             name = AUX_LABEL\n                            ),\n                go.Scatter3d(x=xna[:i], \n                             y=yna[:i], \n                             z=zna[:i],\n                             name = NON_AUX_LABEL\n                            ),\n               ]\n\n        if time in non_aux['time'].to_list():\n            data[1] = go.Scatter3d(x=xna[:i+1], \n                                   y=yna[:i+1], \n                                   z=zna[:i+1],\n                                   name = NON_AUX_LABEL\n                                  )\n            i+=d\n        frames.append(go.Frame(data=data,\n                               traces = [0,1],\n                               name = f'frame_{time}'\n                              ))\n    \n    # LAST FRAME -----------------------------------\n    time = current_event['time'].iloc[-1]\n\n    data = [go.Scatter3d(x=x[:], \n                         y=y[:], \n                         z=z[:],\n                         name = AUX_LABEL\n                        ),\n            go.Scatter3d(x=xna[:], \n                         y=yna[:], \n                         z=zna[:],\n                         name = NON_AUX_LABEL\n                        ),\n           ]\n\n    if time in non_aux['time'].to_list():\n        data[1] = go.Scatter3d(x=xna[:], \n                               y=yna[:], \n                               z=zna[:],\n                               name = NON_AUX_LABEL\n                              )\n        i+=1\n    frames.append(go.Frame(data=data,\n                           traces = [0,1],\n                           name = f'frame_{time}'\n                              ))\n    # LAST FRAME -----------------------------------\n\n    fig.update(frames=frames)\n    #print(\"frames added\")\n    fig.update_layout(updatemenus=[dict(type = \"buttons\",\n                                        buttons = [dict(label=\"Play\", \n                                                        method=\"animate\",\n                                                        args=[None, \n                                                              dict(frame=dict(redraw=True,\n                                                                              duration = ANIMATION_FRAME_DURATION,\n                                                                             ),\n                                                                   fromcurrent = True,\n                                                                   transition=dict(duration = ANIMATION_FRAME_DURATION,\n                                                                                   easing = \"quadratic-in-out\",\n                                                                                  )\n                                                                  )      \n                                                             ]\n                                                       ),\n                                                   dict(label=\"Pause\", \n                                                        method=\"animate\",\n                                                        args=[[None], \n                                                              dict(frame=dict(redraw=False,\n                                                                              duration = 0,\n                                                                             ),\n                                                                   mode=\"immediate\",\n                                                                   transition=dict(duration = 0)\n                                                                  )      \n                                                             ]\n                                                       ),\n                                                  ],\n                                        pad = {\"b\":0,\n                                               \"r\":0,\n                                               \"l\":0,\n                                               \"t\":0\n                                              },\n                                        showactive = False,\n                                        direction = \"right\",\n                                        xanchor = \"auto\",\n                                        yanchor = \"auto\",\n                                        y = -0.20,\n                                        x = -0.05,\n                                       )\n                                  ]\n                     )\n\n    sliders_dict = {\n        \"active\": 0,\n        \"yanchor\": \"top\",\n        \"xanchor\": \"left\",\n        \"currentvalue\": {\n            \"font\": {\"size\": 20},\n            \"prefix\": \"Time(ns):\",\n            \"visible\": True,\n            \"xanchor\": \"right\"\n        },\n        \"transition\": {\"duration\": ANIMATION_FRAME_DURATION, \"easing\": \"cubic-in-out\"},\n        \"pad\": {\"b\": 10, \"t\": 50},\n        \"len\": 0.9,\n        \"x\": 0.1,\n        \"y\": 0,\n        \"steps\": [{\"args\": [[f'frame_{time}'],\n                            {\"frame\": {\"duration\": ANIMATION_FRAME_DURATION, \n                                       \"redraw\": True\n                                      },\n                             \"mode\": \"immediate\",\n                             \"transition\": {\"duration\": ANIMATION_FRAME_DURATION}\n                            }\n                           ],\n                   \"label\": time,\n                   \"method\": \"animate\"\n                  } \n                  for time in current_event['time'] \n                 ]\n    }\n\n    fig.update_layout(sliders=[sliders_dict])\n\n    fig.update_scenes(xaxis_visible=False, \n                      yaxis_visible=False,\n                      zaxis_visible=True\n                     )\n\n    fig.update_layout(autosize=True,\n                      height=800,\n                     )\n\n    fig.update_layout(\n        title={\n            'text': \"<b>Event \"+str(event)+\"</b>\",\n            'y':0.9,\n            'x':0.5,\n            'xanchor': 'center',\n            'yanchor': 'top'},\n        template = \"simple_white\",\n        legend=dict(yanchor=\"top\", y=0.95, xanchor=\"left\", x=0.2),\n    )\n\n    # add DOMs to event\n    fig.add_trace(dom_fig)\n\n    # add path to event\n    fig.add_trace(path_fig)\n\n    iplot(fig)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-02-05T07:23:29.822948Z","iopub.execute_input":"2023-02-05T07:23:29.823667Z","iopub.status.idle":"2023-02-05T07:23:31.210470Z","shell.execute_reply.started":"2023-02-05T07:23:29.823610Z","shell.execute_reply":"2023-02-05T07:23:31.208302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<a class=\"anchor\"  id=\"event_animation\"></a>","metadata":{}},{"cell_type":"code","source":"view_event_by_id(cool_events[0])","metadata":{"execution":{"iopub.status.busy":"2023-02-05T07:23:31.215842Z","iopub.execute_input":"2023-02-05T07:23:31.216218Z","iopub.status.idle":"2023-02-05T07:23:40.477302Z","shell.execute_reply.started":"2023-02-05T07:23:31.216184Z","shell.execute_reply":"2023-02-05T07:23:40.474699Z"},"trusted":true},"execution_count":null,"outputs":[]}]}