Unfolding the muon neutrino energy spectrum with three years of IceCube data

Here presented is all information relevant to the work of unfolding the muon neutrino energy spectrum, developed for and with IceCube. 
Author: Mathis Börner
Contact person: Leonora Kardum
Reviewer: Richard Naab
Reviewer: Tom Gaisser

Introduction

The elusive neutrinos are impossible to detect directly. Instead, when neutrinos engage in an interaction, they create the lepton counterpart of the same flavor whose energy we can reconstruct by monitoring their propagation (and interactions) through ice. These secondary particles create Cherenkov radiation, if propagating at a speed higher than the speed of light in that medium. Resulting Cherenkov photons are captured in IceCube's DOMS and their known time, intensity, and position allow the reconstruction of the event.

Unfolding, or sometimes refered to as deconvolution, is the process of infering a physical quantity's distribution from other relating measures. Reconstructing interactions in IceCube's ice knowing a series of parameters affecting the results of detection enables infering the secondary lepton energy from which the primary neutrino energy can be deduced. Due to a complex nature of detection with tens or hundreds of relevant features, unfolding is a highly sensitive process.

Energy of neutrinos created in interactions of charged cosmic rays with our atmosphere prevail anything that can be created on Earth with today's technology. This gives an unprecedented opportunity in neutrino study. The energy spectrum ranging from GeV to PeV constitues from both the conventional and prompt atmospheric components and the astrophysical component in higher energies. The astrophysical component is expected to be flatter, and originating in extreme astrophysical sources. However, the component has not been characterized and IceCube remains the only experiment to ever detect a neutrino in the PeV range. This makes this particular energy range of special interest.

The analysis starts with a chain aimed at generating a data sample consisting solely of muon neutrinos through separation from the generated background. A clean dataset is of immense importance to the unfolding accuracy. The first section deals with preparation of the data for unfolding, including but not limited to feature preprocessing and handling background. The algorithm is then trained on Monte Carlo data to determine the response of the detector to certain events. The second section describes training the model and unfoulding the spectrum with all concepts neccessary. 

The CC interaction

Neutrinos interact exclusively through the electroweak interaction. Neutrino interacts with the quark in the nucleus by exchanging W or Z boson in a process called deep inelastic scattering. Exchanging the W boson is called a charged current interaction and results in the creation of a charged lepton. In Z exchanges, the hadronization of broken nucleus leads to cascades of radiation. However, leptons are not created and therefore the neutral current interaction does not create track-like events.

In this concept, we are interested only in track-like events caused by neutrino-induced propagating muons.

Detection in IceCube

The IceCube neutrino observatory is located at the geographic south pole made of three different components: the main detector, DeepCore, and IceTop. The main detector consists of 4680 digital optical modules (DOMs) housing a photomultiplier and suspended on 86 strings in ice at depths between 1.5 and 2.5 kilometers. The 8 DeepCore strings were built to lower the energy treshold to as low as 10 GeV. Since the lower energy is not of interest in this work, the DeepCore strings are excluded from the analysis. The total instrumentation size is one cubic kilometer. Photomultipliers detect photons and from the observed Cherenkov light, the cuasing particle's trajectory is reconstructed. 

Figure 1. Visualized reconstruction of an event detected in IceCube. Vertical lines denote the strings, and the timelapse of the event is color coded. Red bubbles are DOMs striked earlier than blue bubbles, from which the direction is inferred. Their size represents intensity. This example is an event from January 8th, 2014 
(RunID: 123751, EventID: 34239163)

Unfolding

If the relevant quantity cannot be determined directly, it ought to be inferred from related measures. Unfolding, or deconvolution, is the process of determining the function in inverse problems. The sought after quantity x, in this case the neutrino energy, is smeared with the detector response matrix A. The measured observables are given with y.
$$ g(y) = \int A(x,y) f(x) dx + b(y) $$ With unfolding, we aim to determine the quantity x from a set y of appropriately chosen observables, with a machine learning model trained to determine energy on a simulation set. 

Preprocessing and separation

Data

The Monte Carlo set used for training consists of the $CORSIKA$ set simulating atmospheric muons and the $nugen$ simulating neutrinos of all flavors. The CORSIKA is composed of three sets, simulating atmospheric with energy from 600 GeV to 100 TeV following a $E^{-2}$ spectrum (set 11499), energies from from 600 GeV to 100 TeV following a $E^{-2.6}$ spectrum (set 11058), and muons with highest energies from 100 TeV to 100 EeV following a $E^{-2}$ spectrum (set 11057). Muon and tau neutrinos are simulated with a $E^{-2}$ spectrum from 100 GeV to 100 PeV (in set 11374 and set 11477), and muon and electron neutrinos following a $E^{-1}$ spectrum from 100 GeV to 10 PeV (in set 11069 and set 11068)
The simulation set is reweighted to follow the distribution given in the results of (Aartsen et al, 2016).  

Measurement period starts at May 15th, 2012 and ends on May 18th, 2015 covering 3 years of operation (IC86-II, IC86-III, IC86-IV). No changes were made to the detector in the given timeframe, therefore datasets can be treated equally. Conditions for Runs of IceCube to be used in this analysis are given below.

GoodRunList
Measurement time > 12.5 min
No. active DOMs in main detector 3120

Burnsample is used to check the consistency of the approach on a smaller subsample. Here, 10 percent of first year (IC86-II) was used with a total of 9967 events.

Dataset Lifetime
IC86-II 331 d 15 h 24 min 
IC86-III 361 d 8 h 5 min 
IC86-IV 368 d 10 h 28 min
Burnsample 32 d 22 h 53 min

The processing of the data up to and including level 3 is handled by the collaboration and contains data reconstructed through various algorithms, of which SplineMPE is regarded as best for directional reconstruction and MuEXDifferentialEnergy and TruncatedEnergy for energy reconstruction.

Separation

Separation is aimed at creating a dataset consisting of only muon events that originated in a charged current interaction of a neutrino and can be split in four main steps:
• Making preliminary cuts
• Selecting relevant observables (which features best represent measured data?)
• Multivariate classification (how clearly can we
determine classification between signal and background?)
• Applying confidence cut (how pure our sample needs to be?) 

Preliminary cuts

• Drop all features that do not exist in both the training (Monte Carlo) and real (observed) datasets
• Drop features representing absolute timestamps
• Drop any features with a constant value, or containing a substantial amount of infinities and NaNs (not a number)
• Retain events with SplineMPE.fit_status == 0
• Apply a zenith cut to events with zenith over 86° (SplineMPE.dir.zenith >= 1.501)
• Drop all events with absolute angle reconstruction error larger than 5° (difference in simulated and reconstructed zenith)
• Drop features with Pearson correlation coefficient > 0.96 with other features
Events with simulated zenith angles lower than 86 degrees can be incorrectly reconstructed as events with sought after angles (over 86 degrees). The left panel on Figure 2. shows the events with a reconstructed zenith angle larger than 86°. The right panel shows the true (simulated) angle of the same set as in the left panel. It can be seen events with an angle lower than 86° can pass the zenith cut if their direction is ill-reconstructed. Therefore, we apply an additional cut to events with a higher than 5° reconstruction error in zenith.

Figure 2. After applying the cut on reconstructed zenith (left panel), it can be seen that the set still contains events with true (simulated) zenith under 86°. Therefore, an additional cut in the training dataset was introduced. Ill-reconstructed events are dropped from the training dataset by applying a cut of maximum 5° on the absolute reconstruction error.


Figure 3. Energy reconstruction error (angle between true simulated direction and the reconstructed direction) depending on the True simulated energy. Reconstruction error rises with energy, however most events belong to low energy region therefore contributing to higher absolute amount of ill-reconstructed events. CDF shows the ratio of events in certain energy ranges, normalized to y axis.


Figure 3a. Energy reconstruction error (angle between true simulated direction and the reconstructed direction) depending on the True simulated energy after applying the cuts. The 90% quantile shows lowering of the absolute reconstructon error for all energy regions in comparison to Level 3 data.

Data vs Monte Carlo mismatch

The aim of Monte Carlo data is to simulate the response of the detector as similar to real environment as possible. This is not always possible as simulations heavily depend on approximations. Given that, we want to use the subset of Monte Carlo that best represents the actual response. If a machine learning model can discriminate between real and simulated data, there exist intrinsical differences in features which allow for this discrimination to happen.

Stratified K-Folded Random Forest is trained to classify given events to real or simulated class with the goal of revealing inconsistencies rather than the purpose of classification. If the classification has an average accuracy higher than 0.5 (random guess), features with highest importance in this classification create discrepancies between simulated and real data. Over an iterative procedure, the forest is trained in classifications and features with high importance are discarded. Note: features with high importance are dropped in this case as classification is not desirable. This is not to be confused with the set of most important features later used in unfolding.

Figure 4. After training the Random Forest classifier, it is aplied to predict the class of real and simulated data iteratively. Classification scores given to particular events over 100 iterations are given in the histogram, where red represents the Burnsample, and green represents a Monte Carlo subset. Features that contribute to differentiating between Data and Monte Carlo will appear as outliers in the Feature importance plot (marked in red).

The Random Forest classifier is constructed of Decision Trees which make the decision based on some criterion, in this case gini impurity. For each Decision Tree, the decrease in the information impurity is calculated for introduction of each feature. The average ratio of impurity decrease over all Trees in a Random Forest is called Feature Importance. It is a relative measure of impurity decrease and attains values from 0 to 1. Sum of importance from a certain feature set will therefore add up to 1.
In this context, Feature importance F is the mean value of calculated importances for each iteration since this procedure is repeated iteratively.
To find features that are considered outliers, Median Absolute Deviation is considered. A feature is discarded from the set if it satisfies $$F_i > k_{Gauss} (\lambda) \cdot MAD$$ where $$ k_{Gauss} (\lambda) = 1.4826 ⋅ CDF_{Gauss}^{-1} (1 − \frac{\alpha}{2})$$ is the scale factor for normal distribuions scaling MAD to the central interval, with $ \alpha = 0.05 $ as the significance level and CDF being the Cumulative Distribution Function.

Figure 5. Histograms depicting distribution of same variable in Monte Carlo and Burnsample data. Monte Carlo data has been scaled to fit the size of the Burnsample. Examples show q_dir_pulses variable from reconstruction algorithms BestTrackDirectHitsE (upper) and BestTrackDirectHitsICE (lower). More examples are shown in Figure 7.

12 features were discarded with this approach that are given in table below. Plots portray the new classification scores and Feature Importances.

Figure 6. Predictions by the trained Random Forest show a better agreement between data and simulation after removing the feature with the MAD approach. Mean importance of remaining 311 observables are given on the right.


Discarded features
SplineMPE_MillipedeHighEnergyMIEFitParams.squared_residuals
SplineMPE_MillipedeHighEnergyMIEFitParams.logl
BestTrack.time
DataChallenge.Correlated
BestTrackDirectHitsE.q_dir_pulses
SplineMPEDirectHitsICE.q_dir_pulses
SplineMPE_MillipedeHighEnergyMIEFitParams.predicted_qtotal
SplineMPE_MillipedeHighEnergyMIEFitParams.chi_squared
SplineMPE_MillipedeHighEnergyMIEFitParams.ndof
SplineMPE_MillipedeHighEnergyMIEFitParams.qtotal
ProjectedQ.width
ProjectedQ.area

Table of features discarded with the MAD approach.

Figure 7. Comparison of simulated and real data on examplary discarded features.

Feature selection

The Minimum Redundancy Maximum Relevance is a dimensionality reduction algorithm encapsulating the necessary information for determining the target variable. In this work, the target variable that is sought to be determined as precisely as possible is the neutrino energy. In the simulation set, the energy of the simulated neutrino is given in the variable MCPrimary.energy.
Information theory is concerned about quantifying the uncertainty present in distributions, by measuring their entropy. MRMR assesses the mutual information both between the features themselves and between features and the sought-after target.

Low entropy is shown by distributions with a peak, i.e. distributions centered around a certain value, while flat distributions with many values of similar probabilities are high in entropy.

If we define the entropy of a distribution as
$$ H = - \sum_X p(x) log p (x) $$

and the mutual entropy of two distributions as
$$ H(X;Y) = - \sum_Y p(y) \sum_X p(x;y) log p(x;y) $$
then their mutual information is defined as
$$I (X;Y) = H(X) - H(X;Y) $$
Mutual information can be thought of as the amount of uncertainty in X that is removed by knowing Y.

Relevance of a feature to a target class C is given with

$$ D = \frac{1}{|S|}\sum_{X_{i}\in S} I (X_i, C) $$
and the redundancy with
$$ R = \frac{1}{|S|^{2}}\sum_{X_{i}, X_{j}\in S} I (X_i, X_{j}) $$
The aim is to minimize correlated features to reduce dimensionality, while not reducing features relevant in determining the target. The MRMR algorithm iteratively adds features to the set S. The set with the highest relevance to redundancy ratio is kept as the most relevant set.

The MRMR approach is used in both parts of this analysis, namely Separation and Unfolding. In the context of Separation, MRMR is applied with the target of events source (Data or Simulation) that are accordingly labelled depending on their origin. After running the MRMR alghoritm, prediction of the trained Random Forest is ran again, and the agreement and Feature Importances of the new set are given below. The AUC score for differentiating between real data and simulation is 0.544 for the final set, close to the desired value of 0.5 .

Figure 8. Classifying score of the trained Random Forest for the feature set after MRMR selection. The Feature Importances show no outliers.

Figure 9. AUC scores of classifier trained to differentiate between real and Monte Carlo data after applying the Pearson coefficient cut, MAD approach and finally the MRMR selection. Since the desired outcome is that the classifier is not able to differentiate between two types of data, we seek a lower AUC score.

Classification

Events originating from muon neutrinos whose direction reconstruction error is less than 5° and upgoing represent the signal class. All other events, including ill-reconstructed and events originating from atmospheric muons, are treated as background.

As opposed to previous analyses, training the signal class solely on well reconstructed muon neutrino events enables the model to recognize events with an appropriate response in the detector with the aim of achieving higher purity. This is justified through a k-fold training of the Classification Random Forest with different events used for signal and background and compared with the AUC score. The highest AUC score is attained by treating upgoing well reconstructed events as signal, and CORSIKA events (atmospheric muons) as background while disregarding downgoing muon neutrino events, electron and tau neutrino events. Figure 10. shows the mean AUC scores (over k folds) of Random Forests trained to classify with different signal and background definitions.

Figure 10. Mean AUC scores of Random Forest classifiers trained with different combinations of signal and background definitions. In the plot, green represents events used as signal and red represents events used as background. Events marked with gray are not considered in the respective classifier.

The spectral index of muon neutrinos in training is −2 compared to the physical expectation of −3.7 to ensure balance in quantity of training. In the opposite case, training on an unbalanced dataset with substantially lower number of high energy events creates a bias.

Using the trained Random Forest for classification, model is applied to the sample again to inspect the classification score given to each event to which an energy dependant cut is applied. Less background is expected in higher energies due to the astrophysical component exhibiting a flatter spectrum. Also, the mean reconstruction error descreases with energy so a more relaxed cut in high energy region can be applied. Less background is expected in higher zenith angles as well, therefore zenith is split in two regions and individual cuts are determined.
The set is split into 100 subsets in the energy of muon $E_\mu$, with 100 equidistant bins from 100 keV to TeV. For each subset, the classification cut is determined mandating a 99.7 purity ver a 100-cross-validation scheme. A cut for each energy subregion is the average of the 100 calculated cuts. Then, 100 obtained cuts over the whole energy range are put together to form an energy dependant classification cut. This is recalculated for different zenith regions aswell. These cuts are visualized as:

Figure 11. Energy dependant cuts for two zenith regions calculated as an average from 100 cuts mandating a 99.7% purity.

Only events with a classification score higher than the cut for their respetive energy and zenith angle are considered as signal in unfolding.

Particle Event rate Fraction
Atmospheric muons $$8.606 \cdot 10^{-6} Hz$$ $$0.2$$
 Muon neutrinos (of which): $$3.505 \cdot 10^{-3} Hz$$ $$99.74$$
Charged current originating
Upward 		$$3.496 \cdot 10^{-3} Hz$$ $$99.48$$
Chargen current originating
Downward 		$$8.864 \cdot 10^{-6} Hz$$ $$0.25 $$
Neutral current originating $$220.8 \cdot 10^{-9} Hz$$ $$6.28 \cdot 10^{-3} $$
Electron neutrinos $$200.7 \cdot 10^{-9} Hz $$ $$5.71 \cdot 10^{-3} $$
Tau neutrinos $$185.5 \cdot 10^{-9} Hz$$ $$5.28 \cdot 10^{-3} $$

Expected event rates from different primaries after application of the described Random Forest classification.

Figure 12. Distribution of classification scores for simulated data (Monte Carlo) in blue, and Burnsample data in green.

Unfolding

The general unfolding problem given in Introduction, aimed at determining the density function f(x) from known g(y) can be simplified by discreticizing these functions, effectively translating the problem to a Linear model of the form
$$ \vec{g} = A_{m,n} \vec{f} $$
where the background b(y) can be dropped due to working with a highly pure set (99.7 ± 0.3)%.

The neutrino energy flux f is not measured directly and is smeared through the stochastic processes involved in its detection. These processes are represented in the matrix A while the measured quantities then become g.

The migration function A, or detector response, is now a matrix of size m,n where m is the dimension of g  and n is the number of bins in the target space f. It indicates the conditional probability to measure that a given true quantity x is measured with y. A has to map the complete detection process, from the propagation and interactions of primary neutrinos and emerged leptons to the whole detector response described in measured features.

Due to the high condition number of matrix A, the problem cannot be solved by simply looking for the inverse matrix, as this usually results in values ith large errors (the higher the condition number of matrix, the higher the uncertainty its inverse can be determined). These type of problems are commonly referred to as ill-conditioned.

Migration matrix can be determined empirically from the simulation (Monte Carlo) set and subsequently divided by the energy density function of the training dataset to make it independent of the training energy distribution.

The searched spectrum is estimated with the Maximum likelihood estimation where all bins are assumed Poissonian. The likelihood for each bin for some expected value $\lambda$ is given with  

The expected value $\lambda$ corresponds to the linear model given above. The likelihood is also multiplied with a regularization term that is explained in the subsequrnt section. For simplicity, the logarithm of the likelihood function is searched for, with constant terms disregarded. This leads to an alternative form of the seaked likelihood:  

Regularization

When using the likelihood function, high variance can arise due to the ill-conditioning of the migration matrix. These fluctuations are supressed through regularization.

The target distribution f is expected to be smooth. Additionally, we know that the muon neutrino energy flux is a combination of power spectra with different indexes, atleast three for the atmospheric, prompt and the astrophysical component. The logarithm of the power spectrum is a linear function with some expected first derivative pointing to its slope, and with a second derivative close to zero. However, since the detector has a limited sensitivity in the low energy region, it leads to a widening of the linear function in the leftest most part of the spectrum which subsequently leads to some small second derivative. To encompass these expectations, a multivariate normal distribution of form

is added for regularization with C being the Thikonov regularization. Since a small second derivative is expected, the normal distribution is centered around 0. This way, the likelihood becomes a prior with expectation zero for the second derivative of the logarithm of the solution which favors power laws. For the covariance of this expression, a diagonal matrix with regularization strength $\tau$ is introduced. This regularization strength can be adjusted and must be chosen so that unphysical solutions are supressed without deteriorating the results.


Certain ranges, and therefore bins, of the target distribution f can be empty. This would lead to a logarithm of zero in the given expression, which is not defined. For this reason, a small offset d is added to the event spectrum. Both the regularization strength and the offset are determined through an iterative method. The muon neutrino energy flux is a product of the event spectrum and the effective area A_eff. Simplified, the likelihood logarithm is now in form

Feature selection

For the purpose of unfolding, a feature selection has to be run again. The MRMR approach is used again, explained in the previous section, but with a target feature now being the primary neutrino energy fom the Monte Carlo set as this is the sought after value.

The final set of features used in unfolding are:
SplineMPEDirectHitsICE.n_dir_doms - The total number of DOMs, which have at least one direct pulse. A direct pulse has a time residual, that is inside the given direct hits time window.
VariousVariables.Cone_Angle - Angle between the reconstruction of the LineFit and the SplineMPE algorithm.
SplineMPECramerRaoParams.variance_theta - variance of the estimated Zenith in the Cramer Rao reconstruction algorithm in squared radians
Borderness.Q_ratio_in_border - Fraction of charge measured in the outermost layer of the detector.
SplineMPETruncatedEnergy_SPICEMie_DOMS_Neutrino.energy - reconstructed neutrino energy in GeV
SplineMPEDirectHitsICB.n_late_doms - the total number of DOMs, which have at least one late pulse. A late pulse has a time residual, that is after the direct hits time window.
Dustyness.n_doms_in_dust - Number of DOMs that have detected photons and are located within the dust layer.
LineFitGeoSplit1Params.n_hits - Number of reconstructed hits used to calculate the fit.
SplineMPETruncatedEnergy_SPICEMie_BINS_MuEres.value - energy resolution of the primary
SplineMPEDirectHitsICC.dir_track_hit_distribution_smoothness - The smoothness value, based on the direct hit DOMs: How uniformly distributed are the direct hit DOM projections onto the track.
SPEFit2GeoSplit1BayesianFitParams.logl - Negative log likelihood for this reconstruction alghoritm.
SplineMPECharacteristicsIC.avg_dom_dist_q_tot_dom - The average DOM distance from the track weighted by the total charge of each DOM.

MCMC sampling

The Monte Carlo Markov Chain is a Monte Carlo sampling method in which the current sample drawn depends on the previous point sampled (hence the Markov chain property). However, any sample depends only on the one sampled before it, and is fully independent of all samples before and the full evolution of the chain. Starting from a random guess, the chain moves to another point by adding a noise from a predetermined distribution centered around zero. If the posteriori probability at the new point is higher than at the previous, the jump is accepted and the point is included in the sample.

The a-posteriori distribution is expected to be determined with the Bayes probability


where the a-priori distribution of f is set to be uniform



MCMC returns points of the a-posteriori distribution for every bin (energy range), and it's median is used as the estimation for the event number in that respective bin. In this work, the energy range is split into 14 bins between 125 GeV and 2 PeV with two additional bins for underflow and overflow.  

Estimation of Systematic uncertainties

Parameters used in the development of the simulation are subject to uncertainty which leads to a neccessity to handle the uncertainty arising from these systematic settings. Considering systematics introduces additional parameters $\xi$ to the migration matrix $$ A \rightarrow A(\xi)$$ In this analysis, following systematics are cosidered: $$ \epsilon_{DOM} - DOM \ efficiency $$ $$ \alpha_{ice}, s_{ice} - relative \ absorption \ coefficient \ and \ relative \ scattering \ coefficient $$ where the ice parameters are considered jointly. Ice properties describe the absorption and scattering degree and are strongly correlated. The DOM efficiency is a parameter of the simulation that determines how many photons are detected depending on the energy deposited in the detector.
For each systematic parameter $\xi_j$, new simulation set is used. Then, for each observable a function $\omega_i (\xi_j)$ can be interpolated that presents the relative change of the bin content depending on the change of the systematic. To consider the full set of parameters, the function is calculated as a product of the individual reweighting functions as in $$ \omega_i (\xi) = \prod_j \omega_i (\xi_j)$$ where j is the set of considered systematic parameters.
There are three sets with three values of the DOM efficiency $$\epsilon_{DOM} ∶ {0.89, 0.99, 1.09} $$ where the middle value is called the baseline. For a higher efficiency, more events would pass the detection treshold. For two non-baseline sets, each observable bin value is divided by the corresponding value in the baseline set producing a function of three points that give observable ratios to the baseline set $$ \omega_{u}^{\epsilon} = \left \{ \frac{g_u (\epsilon_{DOM} = 0.89)}{g_u (\epsilon_{DOM} = 0.99)} , 1 , \frac{g_u (\epsilon_{DOM} = 0.109)}{g_u (\epsilon_{DOM} = 0.99)} \right \} $$ for some observable bin $g_u$. To this, a linear function can be fit of the form $$ \omega_{u} (\epsilon_{DOM}) = m \cdot \epsilon_{DOM} + l $$ which describes the reweighting function in relation to DOM efficiency.
In a similar manner, a function is fit to describe the dependancy on ice parameters. However, since two parameters are considered jointly, the fit function will be a plane. Available sets have the parameters with values of $$ (\alpha_{ice}, s_{ice} ) ∶ {(1.00, 1.00), (1.10, 1.00), (1.00, 1.10), (0.93, 0.93)} $$ where the first set is the baseline. As for the DOM efficiency, the reweighting function for some observable bin $g_u$ will be $$ \omega_{u}^{\alpha, s} = \left \{ 1, \frac{g_u (\alpha_{ice} = 1.10, s_{ice} = 1.00)}{g_u (\alpha_{ice} = 1.10, s_{ice} = 1.00)} , \frac{g_u (\alpha_{ice} = 1.00, s_{ice} = 1.10)}{g_u (\alpha_{ice} = 1.00, s_{ice} = 1.00)}, \frac{g_u (\alpha_{ice} = 0.93, s_{ice} = 0.93)}{g_u (\alpha_{ice} = 1.00, s_{ice} = 1.00)} \right \} $$ and the fit function will now be of the form $$ \omega_{u} (\alpha_{ice}, s_{ice}) = a \cdot \alpha_{ice} + b \cdot s_{ice} + c $$ Finally, the product of the reweighting function is used to multiply the migration matrix $$ A(\xi) = Diag(\omega (\xi)) A $$ where A is the old migration matrix. Examples of these interpolated functions for random observables can be seen in the following plots.

Figure 16. Linear fit to the relative change of an observable in datasets with varied DOM efficiency. The function must pass through 1.00 where the baseline compares to itself.

Figure 17. Twodimensional fit to the relative change of an observable with varied ice scattering and absorption coefficients. The red dot represents point (1.00, 1.00) where the baseline compares to itself.

The effects of systematics on unfolding results have been tested in other works such as: Reproducibility of unfolding results among different unfolding results: Data mining on the rocks : A measurement of the atmospheric muon neutrino flux using IceCube in the 59-string configuration and a novel data mining based approach to unfolding , Application of machine learning algorithms in imaging Cherenkov and neutrino astronomy. Effect of systematic parameters change on unfolding IceCube data: Measurement of the extraterrestrial and atmospheric muon neutrino energy spectrum with IceCube-79 , Entfaltung des Energiespektrums der von IceCube in der 86 String-Konfiguration gemessenen Myon-Neutrinos . Influence of the systematic effects on the muon cross section: Systematic uncertainties of high energy muon propagation using the lepton propagator PROPOSAL

In this analysis, the unfolding operator is optimized by fitting the three systematic parameters of interest. For the blind dataset, the best estimates correspond to $$ \epsilon_{DOM} = 101.7 \pm 0.3$$ $$ \alpha_{ice}, s_{ice} = 97.3 \pm 0.3, 97.5 \pm 0.3 $$

Expected results

Overall, the unfolding is done in the following steps presented previously:

  • Feature selection for energy unfolding
  • Rebinning the observable space
  • Determining regularization parameters
  • MCMC sampling
  • Estimation of systematic uncertanties

The result of MCMC procedure is the neutrino event spectrum.
Due to existing zenith dependancy of the flux models, the unfolding is done over different zenith regions. Since a subdivision of zenith will have less statistics, the results are rebinned to broader binning.
The unfolding is performed over the whole zenith range (86° to 180°) as well as in following zenith regions:

86° to 107° 107° to 130° 130° to 180°

For each considered region, the regularization parameters have to be redetermined as they contain different events.

Test unfolding

To test the procedure, the unfolding has been done on random subsets of Monte Carlo data. Tested Monte Carlo subsets contain 1333638 entries and are unfolded with 250 Markov Chain walkers each making 200000 steps.

Figure 18. Examples of random Monte Carlo subsets unfolded with the given procedure. The median of the 5000000 point distribution in the Markov Chain (from 200000 steps per 250 walkers) is used as a prediction and its corresponding MAD is visualized as uncertainty. The true distribution of the Monte Carlo subset is shown in orange. The underflow and overflow bins are portrayed as shaded area, while red bars show the uncertainty.

Following the same approach, the described Burnsample is unfolded. Features used are given in the 'Feature selection' chapter and the same parameters for the Markov Chain as in Monte Carlo unfolding were used.

Figure 19. Unfolding of the measured data subset recorded for 32 d 22 h 53 min. The red bars are the uncertainty around the Markov Chain prediction. Larger uncertainties in the unfolded Burnsample compared to Monte Carlo unfolding are caused by the limited statistics of the subset, but are expected to shrink with more data.

.