```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, fig.align = "center", warning = FALSE, message = FALSE) ``` ```{r echo=FALSE} library(arrow) library(corrplot) library(data.table) library(DT) library(dplyr) library(ggplot2) library(ggraph) library(knitr) library(igraph) library(plotly) library(png) ``` ## 1) The IceCube Project ### 1.1) The Neutrino Neutrinos are the only particles showing deviation from the Standard Model of Fundamental Interactions (SM). It states that they are massless but they are not [REF4]. There are many natural and human-made neutrino sources. Some are known as our Sun (and most probably the rest of the stars), supernovae, cosmic rays interacting with the atmosphere, nuclear reactors, neutrino accelerators and radioactive decays of unstable isotopes, and some others are only conjectured like black holes, dark matter and gamma ray bursts. In addition, there is the widely accepted cosmic neutrino background generated after the very first second from the Big Bang. This diversity in origin also determines their characteristics, such as flavour or energy, providing information about a broad range of physics phenomena and properties [REF4]. Although neutrinos are the most abundant of the known massive particles in our universe, adding up to 0.3% of its total mass, they are extremely hard to detect due to their tiny interaction cross sections, only through nuclear weak force. For that reason, in order to detect them and study their properties, huge experiments with very large amount of active matter and very low background are needed. One of the most successful technology employed in neutrino experiments is the so called water-Čerenkov. They measure the radiation emitted by charged particles traveling in the water (for instance those originated in a neutrino interaction) with momentum larger than their Čerenkov threshold [REF4].
FIGURE 1: The Standard Model of particle physics (source: Wikipedia) FIGURE 1: The Standard Model of particle physics (source: Wikipedia)
### 1.2) The IceCube Observatory At present, several experiments are in operation: ANTARES (Astronomy with a Neutrino Telescope and Abyss environmental RESearch), BAIKAL and IceCube. The first two are located in the Northern Hemisphere, and in both cases, natural sea or lake water is used as a material means to perform neutrino detection. ANTARES is located in the Mediterranean Sea, about 40 km from the French coast. The project consists of 12 columns anchored to the seabed, at a depth of about 2500m, each containing 75 photomultipliers that collect the signals produced by the Cherenkov effect. The BAIKAL project, located in the lake of the same name in Russia, has a similar structure. The largest current detector is being developed by the IceCube project, which involves about 300 researchers from 45 institutions and research centers in 12 different countries. IceCube , which is located in Antarctica, differs from the other two projects by using ice as the neutrino detection process. At a depth of one kilometer, the pressure is so high that the ice no longer contains air bubbles, which makes it a material well suited to detect the light signals produced during the interaction of neutrinos. The IceCube project occupies a volume of 1km3 and consists of 86 columns containing 5160 photomultipliers (digital optical modules, or DOMs for short), immersed at a depth between 1450m and 2450m. IceCube points all its detectors downwards. The goal is to eliminate most of the muons produced by cosmic rays, using the planet as a filter.
FIGURE 2: The IceCube observatory architecture (source: [REF1])
FIGURE 2: The IceCube observatory architecture (source: [REF1])
Thus, the objective of IceCube is to produce a map of the neutrino flux of the Northern Hemisphere. This information will complement that provided by the two other experiments, ANTARES and BAIKAL, which detect neutrinos in the Southern Hemisphere. ### 1.3) Instrumentation Each DOM includes 12 LEDs on a “flasher board” that produce pulsed light detectable by other DOMs located up to 0.5 km away. The primary purpose of the measurements with these flashers is calibration of the detector. These calibration studies include determining the detector geometry, verifying the calibration of time offsets and the time resolution, verifying the linearity of photon intensity measurement, and extracting the optical properties of the detector ice [REF2].
FIGURE 3: A schematic view of an IceCube Digital Optical Module (source: [REF2])
FIGURE 3: A schematic view of an IceCube Digital Optical Module (source: [REF1])
### 1.4) Influence of Ice on the Model To realize the full potential of the detector, the properties of light propagation in the ice in and around the detector must be well understood. As the photons propagate from the point of emission to the receiving sensor, they are affected by absorption and scattering in the ice. These propagation effects must be considered for both simulation and reconstruction of IceCube data and thus need to be carefully modeled. The important parameters to describe photon propagation in a transparent medium are: the average distance to absorption, the average distance between successive scatters of photons, and the angular distribution of the new direction of a photon at each given scattering point. To determine the ice parameters, dedicated measurements are performed with the IceCube detector. Photons are emitted by the LEDs in DOMs and recorded by other DOMs. A global fit of these data was performed, and the result is a set of scattering and absorption parameters that best describes the full data set. Several dust loggers were used during the deployment of six of the IceCube strings resulting in a survey of the structure of ice dust layers with an effective resolution of approximately 2 mm. The thickness of the ice layers was somewhat arbitrarily chosen to be 10 m. The scattering and absorption coefficients of each ice layer are best interpreted as the average of their true values over the thickness of the ice layer. The chosen thickness of 10 m is smaller than the vertical DOM spacing of 17 m. Due to small depth offsets between the DOMs on different strings, we retain at least 1 receiving DOM per layer. Here is the table of scattering length and absorption length vs. depth, unit is meter (source: [REF2]): ```{r echo=FALSE, cache=FALSE} # Cleaning environment. rm(list = ls()) iceTransparency <- data.table(depth = c(-1398.4, -1408.4, -1418.4, -1428.4, -1438.4, -1448.4, -1458.4, -1468.4, -1478.4, -1488.4, -1498.4, -1508.5, -1518.6, -1528.7, -1538.8, -1548.7, -1558.7, -1568.5, -1578.5, -1588.5, -1598.5, -1608.5, -1618.5, -1628.5, -1638.5, -1648.4, -1658.4, -1668.4, -1678.5, -1688.5, -1698.5, -1708.5, -1718.5, -1728.5, -1738.5, -1748.5, -1758.5, -1768.5, -1778.5, -1788.5, -1798.5, -1808.5, -1818.4, -1828.4, -1838.4, -1848.4, -1858.4, -1868.5, -1878.5, -1888.5, -1898.5, -1908.5, -1918.5, -1928.5, -1938.5, -1948.5, -1958.5, -1968.5, -1978.5, -1988.4, -1998.4, -2008.4, -2018.5, -2028.5, -2038.5, -2048.5, -2058.5, -2068.5, -2078.5, -2088.5, -2098.5, -2108.5, -2118.4, -2128.4, -2138.4, -2148.4, -2158.3, -2168.3, -2178.3, -2188.2, -2198.2, -2208.2, -2218.2, -2228.2, -2238.3, -2248.3, -2258.3, -2268.2, -2278.2, -2288.1, -2298.0, -2308.0, -2318.0, -2328.0, -2338.0, -2348.0, -2357.9, -2367.8, -2377.8, -2387.8, -2397.9, -2408.0, -2418.0, -2428.1, -2438.1, -2448.2, -2458.2, -2468.3, -2478.4, -2488.4), scatteringLength = c(13.2, 14.0, 14.7, 17.0, 16.0, 14.4, 16.0, 20.8, 26.7, 34.7, 39.7, 38.7, 27.8, 16.6, 13.7, 13.5, 15.7, 15.7, 14.7, 17.6, 21.6, 24.0, 20.0, 17.8, 28.9, 36.9, 42.1, 46.5, 45.4, 39.1, 30.6, 26.5, 19.3, 20.8, 20.1, 20.3, 24.5, 33.5, 36.2, 35.4, 32.3, 40.2, 44.7, 34.5, 30.6, 27.5, 19.7, 21.4, 28.8, 38.3, 38.4, 44.2, 50.5, 46.6, 36.8, 26.7, 20.3, 17.4, 16.1, 9.4, 10.6, 13.2, 10.9, 6.8, 5.5, 5.0, 7.2, 9.8, 12.2, 21.1, 54.3, 50.5, 33.5, 34.6, 48.4, 53.2, 46.3, 32.9, 27.4, 30.5, 28.9, 35.1, 39.9, 48.0, 53.3, 54.8, 57.9, 61.1, 76.8, 79.0, 75.6, 75.3, 78.0, 59.4, 51.8, 32.9, 23.9, 28.6, 32.5, 44.5, 56.9, 57.5, 54.3, 61.3, 68.8, 77.6, 79.8, 89.4, 80.7, 56.7), absorptionLength = c(45.1, 48.6, 53.2, 57.6, 57.6, 52.2, 60.1, 74.6, 96.6, 110.5, 135.6, 134.7, 98.2, 64.7, 48.5, 44.3, 54.4, 56.7, 52.1, 60.7, 72.7, 78.9, 68.7, 66.6, 100.0, 128.6, 148.2, 165.7, 156.0, 138.5, 113.9, 90.2, 73.5, 75.9, 67.8, 68.6, 83.8, 119.5, 121.6, 108.3, 113.4, 139.1, 148.1, 122.8, 113.8, 89.9, 71.7, 70.6, 95.9, 116.5, 143.6, 169.4, 178.0, 156.5, 135.3, 103.9, 75.2, 66.2, 53.7, 33.6, 36.2, 44.0, 40.4, 24.9, 20.1, 17.9, 28.4, 34.4, 41.6, 84.4, 173.1, 180.8, 116.7, 120.4, 164.4, 172.8, 149.2, 108.4, 91.1, 98.9, 94.0, 113.1, 134.8, 154.1, 157.6, 180.5, 179.7, 185.2, 227.2, 220.8, 223.9, 256.6, 264.4, 193.7, 159.1, 118.7, 86.2, 104.0, 119.7, 140.6, 203.5, 201.8, 178.2, 206.0, 205.2, 232.1, 259.4, 276.1, 244.3, 185.2), key = "depth") setorder(iceTransparency, -depth) datatable(iceTransparency, rownames = FALSE, options = list(pageLength = 10 ) ) %>% formatRound("depth", digits = 1) %>% formatRound("scatteringLength", digits = 1) %>% formatRound("absorptionLength", digits = 1) ``` These two variables will be associated with the different sensors according to their depth. ### 1.5) The Detection Results The shape of the signals produced in the detector by the light, as well as the amount of energy measured by the DOMs, allow researchers to estimate the energy, direction and sometimes even the flavor of the interacting neutrino.
FIGURE 4: The three distinct neutrino signatures (source: [REF3])

cascade signature

track signature

double bang signature
In the image above, three distinct neutrino signatures are shown: - The "Cascade" signature, on the left, is a typical signature of an electron neutrino, which interacts in the detector by producing an electromagnetic cascade of particles. - The "Track" signature, in the center, is characteristic of a muon neutrino, which interacts by producing a muon as the only secondary particle. This muon then crosses the whole detector, leaving a light line in its wake. - Finally, the "Double Bang" signature, on the right, is characteristic of a tauonic neutrino. The latter interacts in the detector by producing first a hadronic cascade (the first reddish cascade on the image) and then a tau, which decays almost immediately, thus triggering a second cascade (in green on the figure). Unfortunately, neutrinos are not the only particles capable of reaching IceCube\'s sensors. Every day, **millions of muons created by the interaction of cosmic rays with the Earth\'s atmosphere** leave a trace in the detector. A muon is a charged particle similar to the electron, but with a mass 200 times greater than the electron. In addition to muons, a large amount of neutrinos is also produced during these interactions, making it difficult to identify neutrinos of cosmic origin. **What does IceCube see?** - **250 million signals** are detected every day, but only a few hundred of them are neutrinos. - Most of the observed neutrinos were produced in the Earth\'s atmosphere; **only a few dozen neutrinos per year** were produced beyond our solar system. In July 2018, the IceCube Neutrino Observatory announced that they have traced an extremely-high-energy neutrino that hit their detector in September 2017 back to its point of origin in the blazar TXS 0506 +056 located 5.7 billion light-years away in the direction of the constellation Orion. This was the first time that a neutrino detector had been used to locate an object in space, and indicated that a source of cosmic rays had been identified. In 2020, evidence of the Glashow resonance at 2.3σ (formation of the W boson in antineutrino-electron collisions) was announced. In February 2021, a possible detection of a tidal disruption event AT2019dsg was reported and a second candidate AT2019fdr in June 2022. In November 2022, IceCube announced the detection of a neutrino source emitted by the active galactic nucleus of Messier 77. It is the second detection by IceCube after TXS 0506+056, and only the fourth known source including SN1987A and solar neutrinos. OKS 1424+240 and GB9 are others possible candidates. (cf. Wikipedia) ## 2) Explorating Data ```{r echo=FALSE} # Defining data files location. dataDirectory <- "/kaggle/input/icecube-neutrinos-in-deep-ice/" ``` ### 2.1) Locating Data What\'s in the input directory? ```{r echo=FALSE, cache=FALSE} list.files(path = dataDirectory) ``` ### 2.2) The sample_submission.parquet File An example submission with the correct columns and properly ordered event IDs. The sample submission is provided in the parquet format so it can be read quickly but **your final submission must be a csv**. Here\'s the content of this file: ```{r echo=FALSE, cache=FALSE} sampleSubmissionFilename <- file.path(dataDirectory, "sample_submission.parquet") sampleSubmission = read_parquet(sampleSubmissionFilename) datatable(sampleSubmission, rownames = FALSE, options = list(pageLength = 10 ) ) ``` ### 2.3) The sensors_geometry.csv File First of all, let\'s check out sensors information, located in _sensors_geometry.csv_ file. This file contains the **x**, **y**, and **z** positions for each of the 5160 IceCube sensors. The row index corresponds to the sensor_idx feature of pulses. The **x**, **y**, and **z** coordinates are in units of meters, with the origin at the center of the IceCube detector. The coordinate system is right-handed, and the z-axis points upwards when standing at the South Pole. You can convert from these coordinates to azimuth and zenith with the following formulas (here the vector (x,y,z) is normalized): ``` x = cos(azimuth) * sin(zenith) y = sin(azimuth) * sin(zenith) z = cos(zenith) ``` ```{r echo=FALSE, cache=FALSE} sensorGeometryFilename <- file.path(dataDirectory, "sensor_geometry.csv") sensorGeometry = fread(file = sensorGeometryFilename, header = TRUE, sep = ",", na.strings = c("", "NA"), verbose = FALSE) setkey(sensorGeometry, sensor_id) ``` Here are the 10 first records of the _sensors_geometry.csv_ file: ```{r echo=FALSE, cache=FALSE} head(sensorGeometry, n = 10) ``` and the 10 last records: ```{r echo=FALSE, cache=FALSE} tail(sensorGeometry, n = 10) ``` Min/Max values for X coordinates are `r min(sensorGeometry$x)` and `r max(sensorGeometry$x)`. Min/Max values for Y coordinates are `r min(sensorGeometry$y)` and `r max(sensorGeometry$y)`. Min/Max values for Z coordinates are `r min(sensorGeometry$z)` and `r max(sensorGeometry$z)`. Here\'s the content of this file: ```{r echo=FALSE, cache=FALSE} datatable(sensorGeometry, rownames = FALSE, options = list(pageLength = 10 ) ) ``` More details on Digital Optical modules (DOMs): 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. The center (x=0, y=0, z=0) of IceCube observatory is located at a depth of 1950 meters. I will add **depth** and **string** property for all sensors. Eight of these strings at the center of the array were deployed more compactly, with a horizontal separation of about 70 meters and a vertical DOM spacing of 7 meters. This denser configuration forms the DeepCore subdetector, which lowers the neutrino energy threshold to about 10 GeV, creating the opportunity to study neutrino oscillations (source: University of Wisconsin–Madison). The (X, Y) coordinates of the DeepCore sensors are as follows: ``` x = -10.97, y = 6.72 x = -9.68, y = -79.5 x = 31.25, y = -72.93 x = 41.6, y = 35.49 x = 57.2, y = -105.52 x = 72.37, y = -66.6 x = 106.94, y = 27.09 x = 113.19, y = -60.47 ``` ```{r echo=FALSE, cache=FALSE} # Depth of the sensors. # https://arxiv.org/pdf/1301.5361.pdf # Center of IceCube observatory : -1950 meters centerDepth <- 1950 sensorGeometry[, depth := z - centerDepth] # String of the sensors. sensorGeometry[, string := sensor_id %/% 60 + 1] # Adding a column dedicated to the color associated with the sensor. # Coordinates of Deepcore DOMs: # x = -10.97, y = 6.72 # x = -9.68, y = -79.5 # x = 31.25, y = -72.93 # x = 41.6, y = 35.49 # x = 57.2, y = -105.52 # x = 72.37, y = -66.6 # x = 106.94, y = 27.09 # x = 113.19, y = -60.47 colorsSensors = c("orange", "darkred") sensorGeometry[, sensorSection := "IceCube"] sensorGeometry[x==-10.97 & y==6.72, sensorSection := "DeepCore"] sensorGeometry[x==-9.68 & y==-79.5, sensorSection := "DeepCore"] sensorGeometry[x==31.25 & y==-72.93, sensorSection := "DeepCore"] sensorGeometry[x==41.6 & y==35.49, sensorSection := "DeepCore"] sensorGeometry[x==57.2 & y==-105.52, sensorSection := "DeepCore"] # Bug with -66.6 value? the devil is in the details! sensorGeometry[x==72.37, sensorSection := "DeepCore"] sensorGeometry[x==106.94 & y==27.09, sensorSection := "DeepCore"] sensorGeometry[x==113.19 & y==-60.47, sensorSection := "DeepCore"] ``` So, here\'s the sensors data set, including all additional properties: ```{r echo=FALSE, cache=FALSE} # Merging sensorGeometry and iceTransparency data.tables. augmentedSensorGeometry <- iceTransparency[sensorGeometry, roll = "nearest", on = "depth"] setkey(augmentedSensorGeometry, sensor_id) setcolorder(augmentedSensorGeometry, c("sensor_id", "x", "y", "z", "depth", "string", "sensorSection", "scatteringLength", "absorptionLength")) datatable(augmentedSensorGeometry, rownames = FALSE, options = list(pageLength = 10 ) ) ``` Here\'s the 3D representation of the sensors network: ```{r figs, echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} # See https://plotly.com/r/3d-scatter-plots/ IceCubeArchitecture <- plot_ly(augmentedSensorGeometry, x=~x, y=~y, z=~z, color=~sensorSection, size=3, colors=colorsSensors, text = ~paste('Depth:', depth, '
String:', string, '
Scattering:', scatteringLength, '
Absorption:', absorptionLength, '
Sensor:',sensor_id)) IceCubeArchitecture <- IceCubeArchitecture %>% add_markers() IceCubeArchitecture <- IceCubeArchitecture %>% layout(title="FIGURE 5: The 3d representation of the sensors", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') IceCubeArchitecture ``` ### 2.4) The train_meta.parquet File. ```{r echo=FALSE, cache=FALSE} explainTrain <- c("batch_id (int)", "The ID of the batch the event was placed into.", "event_id (int)", "The event ID.", "first_pulse_index (int)", "Index of the first row in the features dataframe belonging to this event.", "last_pulse_index (int)", "Index of the last row in the features dataframe belonging to this event.", "azimuth (float32)", "The azimuth angle in radians of the neutrino. A value between 0 and 2π.", "zenith (float32)", "The zenith angle in radians of the neutrino. A value between 0 and π.") explainTrainTable <- data.frame(matrix(explainTrain, ncol = 2, byrow = T)) names(explainTrainTable) <- c("Field", "Description") knitr::kable(explainTrainTable) ``` The direction vector represented by **azimuth** and **zenith** points to where the neutrino came from. The two columns **azimuth** and **zenith** are target columns, that's why they're only in the _train_meta.parquet_ file. ```{r echo=FALSE, cache=FALSE} trainMetaFilename <- file.path(dataDirectory, "train_meta.parquet") trainMetaSet <- read_parquet(trainMetaFilename) setDT(trainMetaSet) trainMetaSet[, batch_id := as.integer(batch_id)] trainMetaSet[, event_id := as.integer(event_id)] trainKeycols = c("batch_id", "event_id") setkeyv(trainMetaSet, trainKeycols) ``` Here are the 10 first records of the _train_meta.parquet_ file: ```{r echo=FALSE, cache=FALSE} head(trainMetaSet, n = 10) ``` and the 10 last records: ```{r echo=FALSE, cache=FALSE} tail(trainMetaSet, n = 10) ``` The train data set contains `r nrow(trainMetaSet)` rows (`r uniqueN(trainMetaSet$batch_id)` batches and `r uniqueN(trainMetaSet$event_id)` distinct events). **batch_id** values are numbered from `r min(trainMetaSet$batch_id)` to `r max(trainMetaSet$batch_id)`. **event_id** Min/Max values are `r min(trainMetaSet$event_id)` and `r max(trainMetaSet$event_id)`. **first_pulse_index** Min/Max values are `r min(trainMetaSet$first_pulse_index)` and `r max(trainMetaSet$first_pulse_index)`. **last_pulse_index** Min/Max values are `r min(trainMetaSet$last_pulse_index)` and `r max(trainMetaSet$last_pulse_index)`. What do **first_pulse_index** and **last_pulse_index** variables represent exactly? For exemple, for **event_id** = `r trainMetaSet$event_id[1]`, **first_pulse_index** = `r trainMetaSet$first_pulse_index[1]` and **last_pulse_index** = `r trainMetaSet$last_pulse_index[1]`; that means in _batch_1.parquet_ file associated with this event, records `r trainMetaSet$first_pulse_index[1] + 1` to `r trainMetaSet$last_pulse_index[1] + 1` contain data associated with this event. I choosed to display **azimuth** and **zenith** Min/max values with 15 decimals. Min/Max values for **azimuth** are `r sprintf("%.15f", min(trainMetaSet$azimuth))` and `r sprintf("%.15f", max(trainMetaSet$azimuth))`. Min/Max values for **zenith** are `r sprintf("%.15f", min(trainMetaSet$zenith))` and `r sprintf("%.15f", max(trainMetaSet$zenith))`. You can verify all values of **azimuth** are between 0 and 2$\pi$ and all values of **zenith** are between 0 and $\pi$. ```{r echo=FALSE, cache=FALSE} displayedBatchId <- 3 ``` Here\'s the content of _train_meta.parquet_ file for batch_id = `r displayedBatchId`: ```{r echo=FALSE, cache=FALSE} datatable(subset(trainMetaSet, batch_id == displayedBatchId), rownames = FALSE, options = list(pageLength = 10)) ``` #### 2.4.1) Number of events per batch Number of events for **batch_id** = 1: `r nrow(trainMetaSet[batch_id == 1,])`. Number of events for **batch_id** = 2: `r nrow(trainMetaSet[batch_id == 2,])`. Number of events for **batch_id** = 3: `r nrow(trainMetaSet[batch_id == 3,])`. ```{r echo=FALSE, cache=FALSE} commonNumberOfEventsPerBatch <- nrow(trainMetaSet[batch_id == 1,]) ``` Considering these first results, we can look for batches that are not associated with `r commonNumberOfEventsPerBatch` events. ```{r echo=FALSE, cache=FALSE} tmpEventsPerBatch <- trainMetaSet[, .(numberOfEvents = sum(.N)), by = batch_id] head(tmpEventsPerBatch[numberOfEvents != commonNumberOfEventsPerBatch,]) rm(tmpEventsPerBatch) ``` #### 2.4.2) Missing values Number of missing values for **batch_id** variable: `r sum(is.na(trainMetaSet$batch_id))`. Number of missing values for **event_id** variable: `r sum(is.na(trainMetaSet$event_id))`. Number of missing values for **first_pulse_index** variable: `r sum(is.na(trainMetaSet$first_pulse_index))`. Number of missing values for **last_pulse_index** variable: `r sum(is.na(trainMetaSet$last_pulse_index))`. Number of missing values for **azimuth** variable: `r sum(is.na(trainMetaSet$azimuth))`. Number of missing values for **zenith** variable: `r sum(is.na(trainMetaSet$zenith))`. #### 2.4.3) Distribution of **azimuth** values ```{r echo=FALSE, cache=FALSE, fig.height=4} ggplot(trainMetaSet, aes(x=azimuth)) + labs(title = "Figure 6: Distribution of Azimuth Values", x = "Azimuth Angle (radians)", y = "Number of values") + scale_x_continuous(breaks = c(0., pi/2., pi, (3*pi)/2, 2*pi), labels = c("0", "π/2", "π", "3π/2", "2π")) + geom_histogram(color="red", fill="red", binwidth = 0.01) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) ``` #### 2.4.4) Distribution of **zenith** values ```{r echo=FALSE, cache=FALSE, fig.height=4} ggplot(trainMetaSet, aes(x=zenith)) + labs(title = "Figure 7: Distribution of Zenith Values", x = "Zenith Angle (radians)", y = "Number of values") + scale_x_continuous(breaks = c(0., pi/2., pi, (3*pi)/2, 2*pi), labels = c("0", "π/2", "π", "3π/2", "2π")) + geom_histogram(color="red", fill="red", binwidth = 0.01) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) ``` #### 2.4.5) Distribution of pulse durations ```{r echo=FALSE, cache=FALSE, fig.height=4} figure8 <- ggplot(trainMetaSet, aes(x=last_pulse_index - first_pulse_index)) + labs(title = "Figure 8: Distribution of Pulse Durations", x = "Pulse Duration", y = "Number of values") + scale_y_log10(breaks = c(10^2, 10^4, 10^6, 10^8), labels = c(expression(10^2), expression(10^4), expression(10^6), expression(10^8))) + geom_histogram(color="red", fill="red", bins = 1000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure8 ``` ### 2.5) The test_meta.parquet File. ```{r echo=FALSE, cache=FALSE} testMetaFilename <- file.path(dataDirectory, "test_meta.parquet") testMetaSet <- read_parquet(testMetaFilename) setDT(testMetaSet) setkeyv(testMetaSet, c("batch_id", "event_id")) ``` The test data set contains `r uniqueN(testMetaSet$batch_id)` batche(s) and `r uniqueN(testMetaSet$event_id)` event(s). Here’s the content of test_meta.parquet file. ```{r echo=FALSE, cache=FALSE} datatable(testMetaSet, rownames = FALSE, options = list(pageLength = 10)) ``` NB: Other quantities regarding the event, such as the interaction point in x, y, z (vertex position), the neutrino energy, or the interaction type and kinematics are not included in the dataset. ### 2.6) The [train/test]/batch_[n].parquet File. Each batch contains tens of thousands of events. Each event may contain thousands of pulses, each of which is the digitized output from a photomultiplier tube and occupies one row. ```{r echo=FALSE, cache=FALSE} explainBatch <- c("event_id (int)", "The event ID. Saved as the index column in parquet.", "time (int)", "The time of the pulse in nanoseconds in the current event time window. The absolute time of a pulse has no relevance, and only the relative time with respect to other pulses within an event is of relevance.", "sensor_id (int)", "The ID of which of the 5160 IceCube photomultiplier sensors recorded this pulse.", "charge (float32)", "An estimate of the amount of light in the pulse, in units of photoelectrons (p.e.). A physical photon does not exactly result in a measurement of 1 p.e. but rather can take values spread around 1 p.e. As an example, a pulse with charge 2.7 p.e. could quite likely be the result of two or three photons hitting the photomultiplier tube around the same time. This data has float16 precision but is stored as float32 due to limitations of the version of pyarrow the data was prepared with.", "auxiliary (bool)", "If True, the pulse was not fully digitized, is of lower quality, and was more likely to originate from noise. If False, then this pulse was contributed to the trigger decision and the pulse was fully digitized.") explainBatchTable <- data.frame(matrix(explainBatch, ncol = 2, byrow = T)) names(explainBatchTable) <- c("Field", "Description") knitr::kable(explainBatchTable) # Only few batch_* files will be used for the analyze. trainBatch1Filename <- file.path(dataDirectory, "train", "batch_1.parquet") trainBatch1Set <- read_parquet(trainBatch1Filename) setDT(trainBatch1Set) setkeyv(trainBatch1Set, c("event_id", "sensor_id")) trainBatch2Filename <- file.path(dataDirectory, "train", "batch_2.parquet") trainBatch2Set <- read_parquet(trainBatch2Filename) setDT(trainBatch2Set) setkeyv(trainBatch2Set, c("event_id", "sensor_id")) trainBatch659Filename <- file.path(dataDirectory, "train", "batch_659.parquet") trainBatch659Set <- read_parquet(trainBatch659Filename) setDT(trainBatch659Set) setkeyv(trainBatch659Set, c("event_id", "sensor_id")) trainBatch660Filename <- file.path(dataDirectory, "train", "batch_660.parquet") trainBatch660Set <- read_parquet(trainBatch660Filename) setDT(trainBatch660Set) setkeyv(trainBatch660Set, c("event_id", "sensor_id")) trainBatchSet <- rbindlist(list(trainBatch1Set, trainBatch2Set, trainBatch660Set, trainBatch659Set)) ``` Here are the first 10 records of the _batch_1.parquet_ file: ```{r echo=FALSE, cache=FALSE} head(trainBatch1Set, n = 10) ``` and the last 10 records of the file: ```{r echo=FALSE, cache=FALSE} tail(trainBatch1Set, n = 10) ``` For the rest of the analysis, I only use the content of the first two files and the last two files (_batch_1.parquet_, _batch_2.parquet_, _batch_659.parquet_, _batch_660.parquet_). Min/Max values for **event_id** variable are `r min(trainBatchSet$event_id)` and `r max(trainBatchSet$event_id)`. #### 2.6.1) Number of records per event In this section, we\'re looking for the number of records associated with the different events. ```{r echo=FALSE, cache=FALSE, fig.height=4,fig.width=8} tmpRecordsPerEvent <- trainBatchSet[, .(numberOfRecords = sum(.N)), by = event_id] head(tmpRecordsPerEvent) figure9 <- ggplot(tmpRecordsPerEvent, aes(x=numberOfRecords)) + labs(title = "Figure 9: Distribution of Number of Records per Event", x = "Number of records", y = "Count") + scale_x_log10(breaks = c(20, 50, 100, 1000, 10000, 50000, 150000), labels = c(20, 50, 100, 1000, 10000, 50000, 150000)) + scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure9 ``` Her\'s the list of events associated with the minimum number of records: ```{r echo=FALSE, cache=FALSE} tmpRecordsPerEvent[numberOfRecords == min(numberOfRecords),] ``` Here\'s the list of events associated with the maximum number of records: ```{r echo=FALSE, cache=FALSE} tmpRecordsPerEvent[numberOfRecords == max(numberOfRecords),] rm(tmpRecordsPerEvent) ``` #### 2.6.2) Number of sensors per event ```{r echo=FALSE, cache=FALSE} firstEvent <- trainBatchSet[1, event_id] recordsFirstEvent <- nrow(trainBatchSet[event_id == firstEvent,]) ``` The first event of the set is #`r firstEvent`, it's associated with `r recordsFirstEvent` records. It's not too much, we can display them. ```{r echo=FALSE, cache=FALSE} trainBatchSet[event_id==firstEvent,] ``` We notice the same sensor can be "illuminated" several times for a specific event. ```{r echo=FALSE, cache=FALSE} tmpSensorsPerEvent <- trainBatchSet[, .(numberOfSensors = .N, numberOfDistinctSensors = uniqueN(sensor_id)), by = event_id] minSensorsPerEvent <- tmpSensorsPerEvent[numberOfDistinctSensors == min(numberOfDistinctSensors),] ``` The lowest number of distinct sensors per event is `r minSensorsPerEvent[1, numberOfDistinctSensors]`, associated with event #`r minSensorsPerEvent[1, event_id]`. Here are the sensors associated with the lowest number of distinct sensors: ```{r echo=FALSE, cache=FALSE} head(setorderv(tmpSensorsPerEvent[numberOfDistinctSensors < 100,], cols = c("numberOfDistinctSensors"), order = c(1)), n=10) maxSensorsPerEvent <- tmpSensorsPerEvent[numberOfDistinctSensors == max(numberOfDistinctSensors),] ``` The highest number of distinct sensors per event is `r maxSensorsPerEvent[1, numberOfDistinctSensors]`, associated with event #`r maxSensorsPerEvent[1, event_id]`. Here are the sensors associated with the highest number of distinct sensors: ```{r echo=FALSE, cache=FALSE} head(setorderv(tmpSensorsPerEvent[numberOfDistinctSensors > 1000,], cols = c("numberOfDistinctSensors"), order = c(-1)), n=10) ``` So we can display the histograms of non distinct/distinct sensors per event. ```{r echo=FALSE, cache=FALSE, fig.height=4} figure10a <- ggplot(tmpSensorsPerEvent, aes(x=numberOfSensors)) + labs(title = "Figure 10a: Distribution of Number of Non Distinct Sensors per Event", x = "Number of Sensors", y = "Count") + scale_x_log10(breaks = c(10, 20, 50, 100, 500, 1000, 2000, 10000, 100000), labels = c(10, 20, 50, 100, 500, 1000, 2000, expression(10^4), expression(10^5))) + scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure10a figure10b <- ggplot(tmpSensorsPerEvent, aes(x=numberOfDistinctSensors)) + labs(title = "Figure 10b: Distribution of Number of Distinct Sensors per Event", x = "Number of Sensors", y = "Count") + scale_x_log10(breaks = c(10, 20, 50, 100, 500, 1000, 2000), labels = c(10, 20, 50, 100, 500, 1000, 2000)) + scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure10b rm(tmpSensorsPerEvent) ``` #### 2.6.3) Time per event In this section, we study the time associated with events. ```{r echo=FALSE, cache=FALSE} tmpTimePerEvent <- trainBatchSet[, .(timePerEvent = mean(time)), by = event_id] head(tmpTimePerEvent) ``` We can display the histogram of mean time per event. ```{r echo=FALSE, cache=FALSE, fig.height=4} figure11 <- ggplot(tmpTimePerEvent, aes(x=timePerEvent)) + labs(title = "Figure 11: Distribution of Time per Event", x = "Mean", y = "Count") + scale_x_log10(breaks = c(10^4, 1.2*10^4, 1.4*10^4, 2*10^4, 10^5), labels = c(10000, 12000, 14000, 20000, expression(10^5))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure11 ``` Here\'s the list of events associated with the minimum mean time value: ```{r echo=FALSE, cache=FALSE} tmpTimePerEvent[timePerEvent == min(timePerEvent),] ``` Here\'s the list of events associated with the maximum mean time value: ```{r echo=FALSE, cache=FALSE} tmpTimePerEvent[timePerEvent == max(timePerEvent),] rm(tmpTimePerEvent) ``` #### 2.6.4) Charge per event ```{r echo=FALSE, cache=FALSE, fig.height=4} tmpChargePerEvent <- trainBatchSet[, .(chargePerEvent = mean(charge)), by = event_id] head(tmpChargePerEvent) figure12 <- ggplot(tmpChargePerEvent, aes(x=chargePerEvent)) + labs(title = "Figure 12: Distribution of Charge per Event", x = "Mean", y = "Count") + scale_x_log10(breaks = c(0.5, 1, 5, 10, 20), labels = c(0.5, 1, 5, 10, 20)) + scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure12 ``` Here\'s the list of events associated with the minimum mean charge value: ```{r echo=FALSE, cache=FALSE} tmpChargePerEvent[chargePerEvent == min(chargePerEvent),] ``` Here\'s the list of events associated with the maximum mean charge value: ```{r echo=FALSE, cache=FALSE} tmpChargePerEvent[chargePerEvent == max(chargePerEvent),] rm(tmpChargePerEvent) ``` #### 2.6.5) Auxiliary ratio per event Here\'s auxiliary ratio for first events: ```{r echo=FALSE, cache=FALSE} tmpAuxiliaryRatioPerEvent <- trainBatchSet[, .(auxiliaryPerEvent = sum(auxiliary)/.N), by = event_id] head(tmpAuxiliaryRatioPerEvent) ``` ```{r echo=FALSE, cache=FALSE} minAuxiliaryEvent <- tmpAuxiliaryRatioPerEvent[auxiliaryPerEvent == min(auxiliaryPerEvent),] ``` The lowest auxiliary ratio value is `r minAuxiliaryEvent[1, auxiliaryPerEvent]`, associated with event #`r minAuxiliaryEvent[1, event_id]`. Here are the events associated with the lowest auxiliary ratio values: ```{r echo=FALSE, cache=FALSE} head(setorderv(tmpAuxiliaryRatioPerEvent[auxiliaryPerEvent < .2,], cols = c("auxiliaryPerEvent"), order = c(1)), n=10) maxAuxiliaryEvent <- tmpAuxiliaryRatioPerEvent[auxiliaryPerEvent == max(auxiliaryPerEvent),] ``` The highest auxiliary ratio value is `r maxAuxiliaryEvent[1, auxiliaryPerEvent]`, associated with event #`r maxAuxiliaryEvent[1, event_id]`. Here are the events associated with the highest auxiliary ratio values: ```{r echo=FALSE, cache=FALSE} head(setorderv(tmpAuxiliaryRatioPerEvent[auxiliaryPerEvent > .7,], cols = c("auxiliaryPerEvent"), order = c(-1)), n=10) ``` ```{r echo=FALSE, cache=FALSE, fig.height=4} figure13 <- ggplot(tmpAuxiliaryRatioPerEvent, aes(x=auxiliaryPerEvent)) + labs(title = "Figure 13: Distribution of Auxiliary Ratio per Event", x = "Auxiliary Ratio", y = "Count") + #scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + geom_histogram(color="red", fill="red", bins = 10000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure13 rm(tmpAuxiliaryRatioPerEvent) ``` #### 2.6.6) Used/Unused sensors On the IceCube website [REF1], it is written: ``` Once the refreezing process ended, which took a couple of weeks to stabilize, the failure rate of the instrumentation has been extremely low, fewer than 100 of the approximately 5,500 sensors are currently nonoperational. ``` So, even if our selected dataset does not include all data, it is interesting to be able to list all the sensors that are not associated to any event. ```{r echo=FALSE, cache=FALSE, fig.height=4} UnusedSensors <- merge(sensorGeometry, trainBatchSet[, .(batchRecords = .N), by = sensor_id], all.x=TRUE) UnusedSensors[is.na(batchRecords), batchRecords := 0] datatable(UnusedSensors[batchRecords == 0], rownames = FALSE, options = list(pageLength = 10)) rm(UnusedSensors) ``` This result seems to confirm the figure of about 100 sensors out of order in the network. #### 2.6.7) Auxiliary ratio per sensor Here\'s auxiliary ratio for first sensors: ```{r echo=FALSE, cache=FALSE, fig.height=4} auxiliaryRatioPerSensor <- trainBatchSet[, .(auxiliaryPerSensor = sum(auxiliary)/.N), by = sensor_id] setkey(auxiliaryRatioPerSensor, sensor_id) head(auxiliaryRatioPerSensor) figure14 <- ggplot(auxiliaryRatioPerSensor, aes(x=sensor_id, y=auxiliaryPerSensor)) + labs(title = "Figure 14: Auxiliary Ratio per Sensor", x = "Sensor", y = "Auxiliary Ratio") + #scale_y_log10(breaks = c(5, 10, 20, 50, 100, 10^3, 10^4, 10^5, 10^6), labels = c(5, 10, 20, 50, 100, expression(10^3), expression(10^4), expression(10^5), expression(10^6))) + #geom_histogram(color="red", fill="red") + geom_bar(stat = "identity", color="red", fill="red") + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure14 ``` I merge the auxiliary information with sensors location information. ```{r echo=FALSE, cache=FALSE} # Performing JOIN between auxiliaryRatioPerSensor and sensorGeometry tables. IceCubBadSensors <- merge(augmentedSensorGeometry, auxiliaryRatioPerSensor, all.x=TRUE) rm(auxiliaryRatioPerSensor) minAuxiliary <- IceCubBadSensors[auxiliaryPerSensor == min(auxiliaryPerSensor, na.rm = TRUE),] ``` The lowest auxiliary ratio value is `r minAuxiliary[1, auxiliaryPerSensor]`, associated with sensor #`r minAuxiliary[1, sensor_id]`. Here are the sensors associated with the lowest auxiliary ratio values: ```{r echo=FALSE, cache=FALSE} head(setorderv(IceCubBadSensors[auxiliaryPerSensor < .2,], cols = c("auxiliaryPerSensor"), order = c(1)), n=10) maxAuxiliary <- IceCubBadSensors[auxiliaryPerSensor == max(auxiliaryPerSensor, na.rm = TRUE),] ``` The highest auxiliary ratio value is `r maxAuxiliary[1, auxiliaryPerSensor]`, associated with sensor #`r maxAuxiliary[1, sensor_id]`. Here are the sensors associated with the highest auxiliary ratio values: ```{r echo=FALSE, cache=FALSE} head(setorderv(IceCubBadSensors[auxiliaryPerSensor > .7,], cols = c("auxiliaryPerSensor"), order = c(-1)), n=10) ``` Here are the DeepCore sensors associated with the highest auxiliary ratio values: ```{r echo=FALSE, cache=FALSE} head(setorderv(IceCubBadSensors[sensorSection == 'DeepCore' & !is.na(auxiliaryPerSensor),], cols = c("auxiliaryPerSensor"), order = c(-1)), n=10) ``` ```{r echo=FALSE, cache=FALSE, fig.height=8, fig.width=10, fig.align='center'} # ,fig.cap="\\label{fig:figs}FIGURE 15: The 3d representation of the sensors associated with auxiliary ratio"} IceCubeBadSensorsArchitecture <- plot_ly(IceCubBadSensors, x=~x, y=~y, z=~z, color=~sensorSection, colors=colorsSensors, size=~auxiliaryPerSensor, sizes=c(1, 100)*2, text = ~paste('Depth:', depth, '
String:', string, '
Scattering:', scatteringLength, '
Absorption:', absorptionLength, '
Sensor:', sensor_id, "
Auxiliary ratio:", sprintf("%.2f", auxiliaryPerSensor*100), "%")) IceCubeBadSensorsArchitecture <- IceCubeBadSensorsArchitecture %>% add_markers() IceCubeBadSensorsArchitecture <- IceCubeBadSensorsArchitecture %>% layout(title = "FIGURE 15: The 3d representation of the sensors associated with auxiliary ratio", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') IceCubeBadSensorsArchitecture ``` We observe that the highest values are concentrated in the area named "Dust Layer" which is located around depth -2000 above the DeepCore detector; this is probably not a coincidence. Only very few DeepCore sensors are associated with an important auxiliary ratio. The correlogram below (using "spearman" method) shows the influence of the ice quality on the "auxiliary" events recorded by the sensors. Note the very high correlation level of the **scatteringLength** and **absorptionLength** variables. ```{r echo=FALSE, cache=FALSE, fig.height=8, fig.width=10, fig.align='center'} IceCubBadSensors[is.na(auxiliaryPerSensor), auxiliaryPerSensor := 0.] corrSensors <- scale(IceCubBadSensors[, !c("sensor_id", "depth", "string", "sensorSection")]) corrplot(cor(corrSensors, method = "spearman"), method = "number", type = "lower") rm(IceCubBadSensors) ``` ### 2.7) Event Analysis In this section, we will study some events. #### 2.7.1) Event #24 Event #24 is the first event of the data set, it's associated with few sensors. Here are **azimuth** and **zenith** properties of this event: ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} currentEvent <- 24 trainMetaSet[event_id==currentEvent,] currentEventDataSet <- merge(sensorGeometry, trainBatch1Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[is.na(charge), charge := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[!is.na(time), opacity := 0.1] nbAuxRecords <- nrow(currentEventDataSet[auxiliary==TRUE,]) nbRecords <- nrow(currentEventDataSet[event_id==currentEvent,]) currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) # Processing azimuth and zenith information. azimuthEvent <- trainMetaSet[event_id==currentEvent,azimuth] zenithEvent <- trainMetaSet[event_id==currentEvent, zenith] xData = cos(azimuthEvent) * sin(zenithEvent) yData = sin(azimuthEvent) * sin(zenithEvent) zData = cos(zenithEvent) # Displaying graph with auxiliary records. currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 16a: The 3d representation of the event #24", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot ``` The red line is the trajectory of the neutrino based on azimuth and zenith values. Time color indicates arrival time, violet first, yellow last. For this event, we have `r sprintf("%2.2f", (nbAuxRecords/nbRecords)*100)`% auxiliary records; we display the graph after removing auxiliary=TRUE points. ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} # Displaying graph without auxiliary records. currentEventDataSet <- merge(sensorGeometry, trainBatch1Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[auxiliary==TRUE, c("time", "charge") := NA] currentEventDataSet[is.na(charge), charge := 0] #IceCubeEvent24[is.na(time), time := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[time != 0, opacity := 0.1] currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 16b: The 3d representation of the event #24\n (only Auxiliary=FALSE points)", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot rm(currentEventDataSet) ``` It is obvious that with so little information, it is very difficult to estimate a trajectory. So we will study events associated with much more information. #### 2.7.2) Event #2145471859 In this section, I will visualize the event #2145471859. Event #2145471859 is located in #660 batch file of the data set, it\'s associated with lot of sensors. Here are **azimuth** and **zenith** properties of this event: ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} currentEvent <- 2145471859 trainMetaSet[event_id==currentEvent,] currentEventDataSet <- merge(sensorGeometry, trainBatch660Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[is.na(charge), charge := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[!is.na(time), opacity := 0.1] nbAuxRecords <- nrow(currentEventDataSet[event_id==currentEvent & auxiliary==TRUE,]) nbRecords <- nrow(currentEventDataSet[event_id==currentEvent,]) currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) # Processing azimuth and zenith information. azimuthEvent <- trainMetaSet[event_id==currentEvent,azimuth] zenithEvent <- trainMetaSet[event_id==currentEvent, zenith] xData = cos(azimuthEvent) * sin(zenithEvent) yData = sin(azimuthEvent) * sin(zenithEvent) zData = cos(zenithEvent) # Displaying graph with auxiliary records. currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 17a: The 3d representation of the event #2145471859", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot ``` This event seems to be associated with a "double bang" signature, characteristic of a tauonic neutrino. For this event, we have `r sprintf("%2.2f", (nbAuxRecords/nbRecords)*100)`% auxiliary records; we can display the graph with only auxiliary=TRUE points. ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} # Displaying graph without auxiliary records. currentEventDataSet <- merge(sensorGeometry, trainBatch660Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[auxiliary==FALSE, c("time", "charge") := NA] currentEventDataSet[is.na(charge), charge := 0] #IceCubeEvent24[is.na(time), time := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[time != 0, opacity := 0.1] currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 17b: The 3d representation of the event #2145471859\n (only Auxiliary=TRUE points)", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot rm(currentEventDataSet) ``` We see that these points cover a large part of the network and do not allow to define a precise pattern. #### 2.7.3) Event #759693 Event #759693 is also associated with lot of sensors. Here are **azimuth** and **zenith** properties of this event: ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} currentEvent <- 759693 trainMetaSet[event_id==currentEvent,] currentEventDataSet <- merge(sensorGeometry, trainBatch1Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[is.na(charge), charge := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[!is.na(time), opacity := 0.1] nbAuxRecords <- nrow(currentEventDataSet[auxiliary==TRUE,]) nbRecords <- nrow(currentEventDataSet[event_id==currentEvent,]) currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) # Processing azimuth and zenith information. azimuthEvent <- trainMetaSet[event_id==currentEvent,azimuth] zenithEvent <- trainMetaSet[event_id==currentEvent, zenith] xData = cos(azimuthEvent) * sin(zenithEvent) yData = sin(azimuthEvent) * sin(zenithEvent) zData = cos(zenithEvent) # Displaying graph with auxiliary records. currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 18a: The 3d representation of the event #759693", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot rm(currentEventDataSet) ``` This event seems to be associated with a "track" signature, characteristic of a muon neutrino. For this event, we have `r sprintf("%2.2f", (nbAuxRecords/nbRecords)*100)`% auxiliary records; we display the graph with only auxiliary=TRUE points. ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} # Displaying graph without auxiliary records. currentEventDataSet <- merge(sensorGeometry, trainBatch1Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[auxiliary==FALSE, c("time", "charge") := NA] currentEventDataSet[is.na(charge), charge := 0] #IceCubeEvent24[is.na(time), time := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[time != 0, opacity := 0.1] currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 18b: The 3d representation of the event #759693\n (only Auxiliary=TRUE points)", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot ``` In this case, auxiliary points seem more interesting than previous case; several points are associated with a significant charge. #### 2.7.4) Event #5709504 Event #5709504 is also associated with lot of sensors. Here are **azimuth** and **zenith** properties of this event: ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} currentEvent <- 5709504 trainMetaSet[event_id==currentEvent,] currentEventDataSet <- merge(sensorGeometry, trainBatch2Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[is.na(charge), charge := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[!is.na(time), opacity := 0.1] nbAuxRecords <- nrow(currentEventDataSet[auxiliary==TRUE,]) nbRecords <- nrow(currentEventDataSet[event_id==currentEvent,]) currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) # Processing azimuth and zenith information. azimuthEvent <- trainMetaSet[event_id==currentEvent,azimuth] zenithEvent <- trainMetaSet[event_id==currentEvent, zenith] xData = cos(azimuthEvent) * sin(zenithEvent) yData = sin(azimuthEvent) * sin(zenithEvent) zData = cos(zenithEvent) # Displaying graph with auxiliary records. currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 19a: The 3d representation of the event #5709504", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot rm(currentEventDataSet) ``` This event seems to be associated with a "cascade" signature, typical of an electron neutrino. For this event, we have `r sprintf("%2.2f", (nbAuxRecords/nbRecords)*100)`% auxiliary records; we display the graph with only auxiliary=TRUE points. ```{r echo=FALSE, cache=FALSE, fig.height=10, fig.width=10, fig.align='center'} # Displaying graph without auxiliary records. currentEventDataSet <- merge(sensorGeometry, trainBatch2Set[event_id==currentEvent], all.x=TRUE) currentEventDataSet[auxiliary==FALSE, c("time", "charge") := NA] currentEventDataSet[is.na(charge), charge := 0] #IceCubeEvent24[is.na(time), time := 0] currentEventDataSet[, opacity := 0.9] currentEventDataSet[time != 0, opacity := 0.1] currentEventDataSet[, c("event_id", "sensorSection"):=NULL] eventKeycols = c("time", "sensor_id") setkeyv(currentEventDataSet, eventKeycols) currentEventPlot <- plot_ly(currentEventDataSet, x=~x, y=~y, z=~z, type="scatter3d", mode="markers", color=~time, opacity=~opacity, size=~charge, sizes=c(1, 200)*20, text = ~paste('Sensor:', sensor_id, "
Time:", time, "
Charge:", charge, "
Auxiliary:", auxiliary)) currentEventPlot <- currentEventPlot %>% layout(title = "FIGURE 19b: The 3d representation of the event #5709504\n (only Auxiliary=TRUE points)", scene = list(xaxis=list(title='X', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), yaxis=list(title='Y', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2), zaxis=list(title='Z', gridcolor = 'rgb(255, 255, 255)', ticklen = 5, gridwith = 2)), paper_bgcolor = '#d1d4d9', plot_bgcolor = '#d1d4d9') %>% add_trace(x = c(-xData*500, xData*500), y = c(-yData*500, yData*500), z = c(-zData*500, zData*500), type = "scatter3d", mode = "lines", name = "lines", text = ~paste('Azimuth:', sprintf("%.6f", azimuthEvent), "
Zenith:", sprintf("%.6f", zenithEvent)), showlegend = FALSE, inherit=FALSE, line=list(width=10, color='red')) currentEventPlot ``` In this case, the auxiliary points also seem to be associated with a "cascade" representation, around a central point. ## 3) How to Characterize the Events We know that there are three event signatures depending on the type of neutrino detected: - The “Cascade” signature is a typical signature of an electron neutrino. - The “Track” signature is characteristic of a muon neutrino. - The "Double Bang" is characteristic of a tauonic neutrino. It is therefore essential to be able to associate one of these three categories with an event. A fourth category "Unknown" can be defined for events whose signature is difficult to characterize, in particular because the number of associated detectors is too low. In this section, I will consider only auxiliary=FALSE data. ### 3.1) The "Cascade" Signature (Event #5709504) In this configuration, we have a single cluster composed of one or more large charges very close to the initial detection point with secondary charges that decrease in intensity as we move away from this point. ```{r echo=FALSE, cache=FALSE} currentEventDataSet <- trainBatch2Set[event_id == 5709504 & auxiliary == FALSE,] #head(currentEventDataSet) ``` Here are some additional features of this event. Number of records : `r nrow(currentEventDataSet)`. Records associated with the minimum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == min(charge)] ``` Records associated with the maximum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == max(charge)] ``` Quantile information for **charge** variable: ```{r echo=FALSE, cache=FALSE} quantile(currentEventDataSet$charge, probs = c(0, 1, 5, 25, 50, 75, 95, 99, 100)/100) quantile(currentEventDataSet$charge, 99/100) #head(currentEventDataSet[charge >= quantile(currentEventDataSet$charge, 99/100),], n=20) ``` ```{r echo=FALSE, cache=FALSE} figure20a <- ggplot(currentEventDataSet, aes(x=charge)) + labs(title = "Figure 20a: Distribution of Charges", x = "Charge", y = "Number of values") + scale_x_log10(breaks = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000), labels = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000)) + geom_histogram(color="red", fill="red", bins = 1000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure20a figure20b <- ggplot(currentEventDataSet, aes(x=time, y=charge)) + labs(title = "Figure 20b: Distribution of Charges in Function of Time", x = "Time", y = "Charge") + scale_y_log10(breaks = c(0, 1, 10, 100, 500, 1000, 2000), labels = c(0, 1, 10, 100, 500, 1000, 2000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure20b figure20c <- ggplot(currentEventDataSet[charge >= quantile(charge, 75/100),], aes(x=time, y=charge)) + labs(title = "Figure 20c: Distribution of Charges in Function of Time\n (Quantile greater than 75%)", x = "Time", y = "Charge") + scale_y_log10(breaks = c(200, 300, 400, 500, 1000, 2000), labels = c(200, 300, 400, 500, 1000, 2000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure20c ``` ### 3.2) The "Track" Signature (Event #759693) In this configuration, we do not really have separate clusters but a series of charges distributed along the path followed by the generated muon in the network. ```{r echo=FALSE, cache=FALSE} currentEventDataSet <- trainBatch1Set[event_id == 759693 & auxiliary == FALSE,] #head(currentEventDataSet) ``` Here are some additional features of this event. Number of records : `r nrow(currentEventDataSet)`. Records associated with the minimum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == min(charge)] ``` Records associated with the maximum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == max(charge)] ``` Quantile information for **charge** variable: ```{r echo=FALSE, cache=FALSE} quantile(currentEventDataSet$charge, probs = c(0, 1, 5, 25, 50, 75, 95, 99, 100)/100) quantile(currentEventDataSet$charge, 99/100) #head(currentEventDataSet[charge >= quantile(currentEventDataSet$charge, 99/100),], n=20) ``` ```{r echo=FALSE, cache=FALSE} figure21a <- ggplot(currentEventDataSet, aes(x=charge)) + labs(title = "Figure 21a: Distribution of Charges", x = "Charge", y = "Number of values") + scale_x_log10(breaks = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000), labels = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000)) + geom_histogram(color="red", fill="red", bins = 1000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure21a figure21b <- ggplot(currentEventDataSet, aes(x=time, y=charge)) + labs(title = "Figure 21b: Distribution of Charges in Function of Time", x = "Time", y = "Charge") + scale_y_log10(breaks = c(0, 1, 10, 100, 500, 1000), labels = c(0, 1, 10, 100, 500, 1000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure21b figure21c <- ggplot(currentEventDataSet[charge >= quantile(charge, 75/100),], aes(x=time, y=charge)) + labs(title = "Figure 21c: Distribution of Charges in Function of Time\n (Quantile greater than 75%)", x = "Time", y = "Charge") + scale_y_log10(breaks = c(200, 300, 400, 500, 1000, 2000), labels = c(200, 300, 400, 500, 1000, 2000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure21c ``` ### 3.3) The "Double Bang" Signature (Event #2145471859) The "Double Bang" signature is characteristic of a tauon neutrino. The latter interacts in the detector by producing first a hadronic cascade and then a tau, which decays almost immediately, thus triggering a second cascade. ```{r echo=FALSE, cache=FALSE} currentEventDataSet <- trainBatch660Set[event_id == 2145471859 & auxiliary == FALSE,] #head(currentEventDataSet) ``` Here are some additional features of this event. Number of records : `r nrow(currentEventDataSet)`. Records associated with the minimum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == min(charge)] ``` Records associated with the maximum value for **charge**: ```{r echo=FALSE, cache=FALSE} currentEventDataSet[charge == max(charge)] ``` Quantile information for **charge** variable: ```{r echo=FALSE, cache=FALSE} quantile(currentEventDataSet$charge, probs = c(0, 1, 5, 25, 50, 75, 95, 99, 100)/100) quantile(currentEventDataSet$charge, 99/100) #head(currentEventDataSet[charge >= quantile(currentEventDataSet$charge, 99/100),], n=20) ``` ```{r echo=FALSE, cache=FALSE} figure22a <- ggplot(currentEventDataSet, aes(x=charge)) + labs(title = "Figure 22a: Distribution of Charges", x = "Charge", y = "Number of values") + scale_x_log10(breaks = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000), labels = c(0.5, 1, 2, 3, 4, 5, 10, 100, 1000, 2000)) + geom_histogram(color="red", fill="red", bins = 1000) + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure22a figure22b <- ggplot(currentEventDataSet, aes(x=time, y=charge)) + labs(title = "Figure 22b: Distribution of Charges in Function of Time", x = "Time", y = "Charge") + scale_y_log10(breaks = c(0, 1, 10, 100, 500, 1000, 2000), labels = c(0, 1, 10, 100, 500, 1000, 2000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure22b figure22c <- ggplot(currentEventDataSet[charge >= quantile(charge, 75/100),], aes(x=time, y=charge)) + labs(title = "Figure 22c: Distribution of Charges in Function of Time\n (Quantile greater than 75%)", x = "Time", y = "Charge") + scale_y_log10(breaks = c(200, 300, 400, 500, 1000, 2000), labels = c(200, 300, 400, 500, 1000, 2000)) + geom_point(size=1, shape=3, color='red') + geom_smooth() + theme_bw() + theme(plot.title = element_text(hjust = 0.5)) figure22c ``` ### 4) References - [REF1] The website of the IceCube Neutrino Observatory: (https://icecube.wisc.edu/science/icecube/) - [REF2] Measurement of South Pole ice transparency with the IceCube LED calibration system (https://arxiv.org/pdf/1301.5361.pdf). - [REF3] The IceCube Masterclass website (french text): (https://masterclass.icecube.wisc.edu/fr/decouvrez/decouvrez-icecube-et-les-neutrinos) - [REF4] Neutrino Physics in Present and Future Kamioka Water-Čerenkov Detectors with Neutron Tagging, Pablo Fernández Menéndez (https://doi.org/10.1007/978-3-319-95086-0). - I drew my inspiration from this work: https://www.kaggle.com/code/utm529fg/eng-icecube-eda-understanding-train-data