\chapter{Data analysis}
\label{cha:data}
In lattice QCD, a large amount of data generated by Markow-chain-Monte-Carlo algorithms has to be analyzed. Statistical and systematic effects play a major role here and have to be considered appropriately. Fortunately, a broad range of mathematical concepts and methods exist that can be adopted for use in lattice QCD. A small fraction of these methods will be explained in the following chapter.
\section{Treating of statistical uncertainties}
\label{sec:stat_err}
We analyze data from a $N$-element Markow chain that samples the distribution that is given by the path integral. Due to the fact that only a finite subspace of the phase space is sampled, a statistical uncertainty is introduced, which can be estimated in several ways.\\
The elements of the Markow chain, which are gauge configurations $U_i$ in the context of lattice gauge theory, are used to measure observables $O[U_i]$. The expectation value of the observable is given by the mean value along the Markow chain:
\begin{align}
  \overline{O} = \frac{1}{N}\sum_{i = 1}^N O[U_i].
\end{align}
Of course, the distribution of the measurements has a width (standard deviation), which is given by the square root of the variance of the distribution. The variance is given by the mean value of the squared difference of the expectation value and the measurements, but the normalization is changed $\frac{1}{N-1}$ because the number of degrees of freedom is already reduced by one when computing the mean value. So the variance is given by
\begin{align}
\label{eq:variance}
  \sigma_O^2 = \frac{1}{N-1}\sum_{i = 1}^N \left(O[U_i]- \overline{O}\right)^2.
\end{align}
On a finite Markow chain, it is only possible to compute estimators for expectation value and variance. Of course the true value is only reached in the case of $N \to \infty$. But in most cases it can be assumed that $N$ is large enough to give a realistic value for those quantities.\\ 
In a Markow chain it is unavoidable that consecutive elements are dependent in some way. This effect is called autocorrelation and is a huge issue in simulation of lattice gauge theory. For example topological observables are plagued by large autocorrelation effects near the continuum limit~\cite{Eichhorn:2023uge}. This affects the estimation of expectation values and errors. The effect of autocorrelation can be quantified using the autocorrelation function (ACF) that gives (up to a normalization of $\sigma_O^2$) the correlation (statistical relationship) of two measurements on the Markow chain that are separated by $t$ steps:
\begin{align}
  \Gamma_O(t) = \lim_{K \to \infty}\frac{1}{K}\sum_{i = 1}^K\left(O[U_{i+t}] - \overline{O} \right)\left(O[U_{i}] - \overline{O} \right)
\end{align}
There is some circumstantial evidence that the autocorrelation is given by a sum of decaying exponential functions $A_i\exp{\frac{t}{\tau_i}}$ ($\tau_i > \tau_{i+1}$)~\cite{Schaefer_2011}. The most relevant contribution to the autocorrelation effects originates from the case $i = 0$ and we can assume that the other contributions are suppressed exponentially. It follows that the autocorrelation time (i.e., the time until the autocorrelation reaches $1/e$) can be estimated by the integral of the ACF. Since we only have access to the ACF for integers $t$, the integral is replaced by a sum. This quantity is called \textit{integrated autocorrelation time} and can be written as
\begin{align}
\label{eq:ACT}
  \tau_{int}(O) = \frac{1}{2} + \sum_{t=1}^\infty \frac{\Gamma_O(t)}{\Gamma_O(0)} = \frac{1}{2}\sum_{t = 1}^\infty\rho_O(t),
\end{align}
with $\rho_O(t) = \frac{\Gamma_O(t)}{\Gamma_O(0)}$ being the normalized autocorrelation function.
In general, one has $N<\infty$, thus it is not possible to compute the whole series, but from some point the ACF is dominated by noise. A simple  In this case, it is usual practice to stop the summation as soon as $\Gamma_O(t)< 0$. \\
The error of the normalized autocorrelation function can be computed via
\begin{align}
  \delta \rho_O(t) = \sqrt{\frac{1}{N}\sum_{m = 1}^W \left(\rho(m+t) + \rho(m-t) - 2\rho(m)\rho(t)\right)^2}
\end{align}
using the convention $\rho(-t) = \rho(t)$\cite{Wolff:2003sm,Joswig:2022qfe}.\\
With this proper estimation of autocorrelation effects it is finally possible to estimate the error on the mean value:
\begin{align}
  \delta O = \sqrt{\frac{2\sigma_O^2\tau_{int}(O)}{N}}.
\end{align}
For fully uncorrelated measurements equation \eqref{eq:ACT} will give an integrated autocorrelation time of $\frac{1}{2}$. Therefore, the naive error on the mean values is given by
\begin{align}
\label{eq:naive_error}
  \left(\delta O\right)_\mathrm{naive} = \frac{\sigma_O}{\sqrt{N}}
\end{align}
and can serve as a lower bound for the error on the primary observable~\cite{Wolff:2003sm,schaefer2009investigating,Schaefer_2011,Joswig:2022qfe}.
\subsection{Blocking}
If a set of data has to be analyzed, blocking is a common strategy. It is done via replacing a subset of the data whose constitutions are similar in some parameters by their mean value~\cite{blocking}. This helps to estimate the statistical uncertainties properly. In the context of Markow chains one often blocks data that are related by autocorrelation effects. This is done by choosing $b \approx 2\tau$ and making blocks of length $b$ in the Markow chain, where $\tau$ is the autocorrelation time of the observable. The length of the chain then reduces from $N$ to $N/b$, so one has to take care that $b$ is a divisor of $N$. If this is not the cases, one can discard some on the measurements. Afterwards the blocks can be treated as independent measurements and the error can be computed naively (compare equation \eqref{eq:naive_error}). The blocking procedure is sketched in Figure \ref{fig:Blocking}.\\
\begin{figure}[H]
    \tikzstyle{block} = [rectangle, draw, fill=gray!20, text width=5.5em, text centered, rounded corners, minimum height=4em]
    \centering
    \begin{tikzpicture}
        \node[block] at (0,0) {\small Measure-ments \\ \tiny $x_{1},...,x_{N}$};
        \draw[->] (1.2,0.75) -- (2.75,2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Central sample};
        \draw[->] (1.2,0) -- (2.75,0) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Blocking};
        \draw[->] (1.2,-0.75) -- (2.75,-2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Blocking};
        % 	\draw[->] (5.2,0) -- (6.75,0);
        \draw[->] (5.2,2) -- (6.75,2);
        \draw[->] (5.2,0) -- (6.75,0);
        \draw[->] (5.2,-2) -- (6.75,-2);
        %\draw[decoration={brace,amplitude=10pt,raise=5pt},decorate] (7.4,2.55) -- node[right=14pt] {$\hat{O}_0\pm\sqrt{\frac{\text{Var}(O_{1,\dots,N/b})}{N/b-1}}$} (7.4,-2.55);%{$\hat{O}\pm \delta o$} (7.4,-2.55);
        \node[block] at (4,2) {\small Central sample \\ \tiny $x_{1},x_{2},...,x_{N}$};
        \node[block] at (4,0) {\small Block $1$ \\ \tiny $x_{1},...,x_{b}$};
        \node[block] at (4,-2) {\small Block $n$ \\ \tiny $x_{nb + 1},...,x_{(n+1)b}$};
        %\node[block] at (10,0) {Uncertainty};
        %
        \node at (7.2,2) {$\hat{O}_0$};
        \node at (7.2,0) {$\hat{O}_1$};
        \node at (7.2,-2) {$\hat{O}_n$};
        
        \draw[->] (9,2) -- (9.55,2);
        \node at (12,2) {$\langle O \rangle = \hat O_0$};
        % \node at (8,-1) {$\scalebox{5}{\}}$};
        \draw[decoration={brace,amplitude=7pt},decorate,very thick,xshift=10pt] (7.2, 0.2) -- (7.2, -2.2);
        \draw[->] (9,-1) -- (9.55,-1);
        \node at (12,-1) {$\delta O =\sqrt{\frac{1}{N(N-1)} \sum_{i = 1}^N (\hat O_i - \hat O_0)^2}$};
    \end{tikzpicture}
    \caption[Visualization of blocking]{Visualization of blocking procedure.}
    \label{fig:Blocking}
\end{figure}
Of course the choice of $\tau$ is not unambiguous. Typically, one chooses $\tau$ such that it is larger than the integrated autocorrelation time of a quantity with large autocorrelation effects~\cite{Borsanyi_2021}, for example the topological charge, which has shown to have dominating autocorrelation times~\cite{Eichhorn:2023uge,Eichhorn:2021ccz,Eichhorn:2022wxn,Rouenhoff:2022seh,Durr:2025qtq}.\\
Another way to estimate a proper value for $b$ is measuring the variance $\sigma^2$ of certain observables in dependence of $b$. If one finds a region where $\sigma^2 \propto \frac{1}{b}$ one can choose a $b$ from this region for an appropriate block size~\cite{GattringerLang}.\\
A disadvantage is that a pure blocking procedure is unpractical when computing the error of secondary observables, i.e., $f(O)$. In this case methods like Gaussian error propagation or automatic differentiation~\cite{Ramos:2018vgu,Ramos:2020scv} have to be applied, which cannot be generalized easily to arbitrary secondary observables~\cite{GattringerLang}. Therefore, one can use one of two other methods: Bootstrap and Jackknife. 
\subsection{Bootstrap resampling}
The method of bootstrapping uses the assumption that the sample of measurements, that we have access to, is representative for the whole set of possible measurements. So one draws (with replacement) set with $N$ elements (the same number as there are measurements in total) out of the sample and computes the corresponding mean value. This procedure is done $B$ times. These $B$ mean values can be either used to compute the error on the primary observable by calculating the standard deviation of the bootstrap samples or they can be used to compute secondary observables by evaluating $f(\hat O_i)$ as shown in Figure~\ref{fig:Bootstrap_resampling}. The error on the secondary observable can again be computed with the standard deviation.\\
To suppress autocorrelation effects one can also use a blocked bootstrapping. This means that a blocking procedure is applied to the data. The resulting mean values of the blocks are than sampled using the bootstrapping procedure~\cite{Michael_1994}.
\begin{figure}[H]
	\tikzstyle{block} = [rectangle, draw, fill=gray!20, text width=5.5em, text centered, rounded corners, minimum height=4em]
	\centering
	\begin{tikzpicture}
	    \node[block] at (0,0) {\small Measure-ments \\ \tiny $x_{1},...,x_{N}$};
		\draw[->] (1.2,0.75) -- (2.75,2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Central sample};
		\draw[->] (1.2,0) -- (2.75,0) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Resample};
		\draw[->] (1.2,-0.75) -- (2.75,-2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Resample};
		% 	\draw[->] (5.2,0) -- (6.75,0);
		\draw[->] (5.2,2) -- (6.75,2);
		\draw[->] (5.2,0) -- (6.75,0);
		\draw[->] (5.2,-2) -- (6.75,-2);
		%\draw[decoration={brace,amplitude=10pt,raise=5pt},decorate] (7.4,2.55) -- node[right=14pt] {$O_{0} \pm \sqrt{\text{Var}(O_{1,\dots,B})}$} (7.4,-2.55);
		\node[block] at (4,2) {\small Central Sample \\ \tiny $x_{1}, x_{2},..., x_{N}$};
		\node[block] at (4,0) {\small Sample $n$ \\ \tiny $x_{1,1},..., x_{1,N}$};
		\node[block] at (4,-2) {\small Sample $B$ \\ \tiny $x_{B,1},..., x_{B,N}$};
		%\node[block] at (10,0) {Uncertainty};
		%
		\node at (7.2,2) {$O_{0}$};
		\node at (7.2,0) {$O_{1}$};
		\node at (7.2,-2) {$O_{B}$};

        \draw[->] (9,2) -- (9.55,2);
        \node at (12,2) {$\langle f(O) \rangle = f(\hat O_0)$};
        % \node at (8,-1) {$\scalebox{5}{\}}$};
        \draw[decoration={brace,amplitude=7pt},decorate,very thick,xshift=10pt] (7.2, 0.2) -- (7.2, -2.2);
        \draw[->] (9,-1) -- (9.55,-1);
        \node at (12,-1) {$\delta f(O) =\sqrt{\frac{\sum_{i = 1}^B (f(\hat O_i) - f(\hat O_0))^2}{B-1} }$};
	\end{tikzpicture}
	\caption[Visualization of bootstrap resampling]{Visualization of bootstrap resampling.}
	\label{fig:Bootstrap_resampling}
\end{figure}
\subsection{Jackknife resampling}
\label{sec:Jackknife}
An alternative to the bootstrapping method is the Jackknife resampling, which is often referred to as ``linearized bootstrapping''. The Jackknife samples are constructed in the following way: To construct the first sample all measurements except for the first one are averaged, for the second sample all measurements except for the second are averaged and so on. Again on each Jackknife sample (and of course on the central sample), one can evaluate secondary observables $f(O)$. The value of the central sample can again be taken as the central value of the observable while the error can be computed from the standard deviation of the $J$ Jackknife samples according to Figure~\ref{fig:Jackknife_resampling}~\cite{10.1093/biomet/43.3-4.353,tukey}.\\
Similar to the Bootstrapping procedure one can include a blocking procedure before computing the Jackknife samples to handle autocorrelation effects. This also the method that is used in the papers the second part of this thesis is based on~\cite{Boccaletti:2024guq,FKAPI,Frech:2023bzy,Frech:2025fjp}.

\begin{figure}[H]
    \tikzstyle{block} = [rectangle, draw, fill=gray!20, text width=5.5em, text centered, rounded corners, minimum height=4em]
    \centering
    \begin{tikzpicture}
    	\node[block] at (0,0) {\small Measure-ments \\ \tiny $x_{1},...,x_{J}$};
    	\draw[->] (1.2,0.75) -- (2.75,2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Central sample};
    	\draw[->] (1.2,0) -- (2.75,0) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Resample};
    	\draw[->] (1.2,-0.75) -- (2.75,-2) node[fill=white,rounded corners=2pt,inner sep=1pt,midway,sloped] {\tiny Resample};
    	% 	\draw[->] (5.2,0) -- (6.75,0);
    	\draw[->] (5.2,2) -- (6.75,2);
    	\draw[->] (5.2,0) -- (6.75,0);
    	\draw[->] (5.2,-2) -- (6.75,-2);
    	%\draw[decoration={brace,amplitude=10pt,raise=5pt},decorate] (7.4,2.55) -- node[right=14pt] {$\hat{O}_0\pm \sqrt{J-1}\sqrt{\text{Var}(O_{1,\dots,J})}$} (7.4,-2.55);%{$\hat{O}\pm \delta o$} (7.4,-2.55);
    	\node[block] at (4,2) {\small Central sample \\ \tiny $x_{1},x_{2},...,x_{J}$};
    	\node[block] at (4,0) {\small Sample $1$ \\ \tiny $x_{2},...,x_{J}$};
        \node[block] at (4,-2) {\small Sample $n$ \\\tiny$...,x_{n-1},x_{n+1},...$};
    	%\node[block] at (10,0) {Uncertainty};
    	%
    	\node at (7.2,2) {$\hat{O}_0$};
    	\node at (7.2,0) {$\hat{O}_1$};
    	\node at (7.2,-2) {$\hat{O}_n$};
    
        \draw[->] (9,2) -- (9.55,2);
        \node at (12,2) {$\langle O \rangle = \hat O_0$};
        % \node at (8,-1) {$\scalebox{5}{\}}$};
        \draw[decoration={brace,amplitude=7pt},decorate,very thick,xshift=10pt] (7.2, 0.2) -- (7.2, -2.2);
        \draw[->] (9,-1) -- (9.55,-1);
        \node at (12,-1) {$\delta O =\sqrt{\frac{J-1}{J}\sum_{i = 1}^J (\hat O_i - \hat O_0)^2}$};
    \end{tikzpicture}
	\caption[Visualization of Jackknife resampling]{Visualization of Jackknife resampling.}
	\label{fig:Jackknife_resampling}
\end{figure}
\section{Treating of systematic uncertainties}
\label{sec:syst_err}
In the context of lattice gauge theory systematic errors are the uncertainties that do not originate form the Monte-Carlo evaluation of the path integral. Typically they are of the same order of magnitude than the statistical error~\cite{Borsanyi_2021,Boccaletti:2024guq,FKAPI}. Systematic errors can have different origins. Typically one of most dominant error sources is the choice of different fitting methods and models. Fits are used to extrapolate the data to the continuum and the physical point (at which the simulation parameters take the physical values), but they also can be used to estimate the energies from correlation functions. In general, one does not know the correct model or method so a variety of different forms is used and the corresponding systematic error is estimated from the width of the resulting distribution. \\
Another dominant source of systematic errors is the setting of the scale and the physical point. Typically, a set of physical observables is computed in the lattice simulation and then fitted to the experimental (or physical) value of the corresponding observable. Of course every uncertainty in the physical value translates in an uncertainty in the final result. So what one typically does it to perform the analysis not using the exact physical value, but with the upper and lower $1\sigma$-deviations. From the spread of the two results, the corresponding systematic error can be estimated~\cite{Borsanyi_2021}.\\
Lattice computations are always done in a finite volume, so there are finite volume effects that have to be accounted for. Fortunately, the finite volume effects are often small enough that a conservative estimation of the error associated to the finite volume correction does not dominate the overall error of the analysis. Of course, there are even more sources of systematic uncertainties like the strong-isospin breaking effect, or effects from other interactions like QED.

\section{Fitting procedure}
\label{sec:fits}
A fit is a method to tune the parameters of a mathematical model in such a way that it aligns with the data from an experiment or a simulation. Usually fitting is done by minimizing some quantity for example the likelihood or the $\chi^2$.\\
Independent variables, for example lattice spacing $a$ (in case of continuum extrapolations) or Euclidean time $t$ (in case of mass/energy plateaus) will be labeled by $x_i$ within this section. Dependent variables (observables such as masses or decay constants) are labeled by $y_i$. The fit models are given by function $f_i$ (mostly polynomials or exponential functions in $x$) and the corresponding parameters that have to be tuned are called $c_i$.
\subsection{\texorpdfstring{$\chi^2$}{chi-squared}}
\label{sec:chi2}
The standard way to fit a model to a set of measurements is the minimization of the $\chi^2$. For a proper definition we need to generalize the concept of the variance mentioned in equation \eqref{eq:variance}: The covariance matrix of a set of statistical data $y_i$ is given by
\begin{align}
  \textrm{Cov}_{ij} = \frac{N}{N-1}\sum_{i = 1}^N \overline{\left(y_i - \overline{y_i}\right)\left(y_j - \overline{y_j}\right)}.
\end{align}
The $i$th diagonal element is given by variance of the quantity $y_i$. The off-diagonal element $(ij)$ of $\textrm{Cov}$ describes how the statistical fluctuation of $y_i$ affect the statistical fluctuation of $y_j$. In some contexts it might be also useful to remove the statistical fluctuation from this expression and only to look at the statistical correlations of certain observables. Therefore, one can define correlation matrix~\cite{Michael_1994,Bruno_2023}
\begin{align}
  \textrm{Cor}_{jk} = \frac{\textrm{Cov}_{jk}}{\sigma_i \sigma_j}
\end{align}
by normalizing the covariance matrix with the standard deviation of the corresponding line and column. If one explicitly wants to ignore all effects originating of correlation effects, one can set all off-diagonal elements to $0$ and receives
\begin{align}
  \Sigma_{ij} = \sigma_i^2 \delta_{ij}.
\end{align}
In this case, the correlation matrix will given by $\textrm{Cor}_{ij} = \delta_{ij}$. To fit a model $f$ depending on parameters $c = (c_i)_i$ and independent variables $x = (x_i)_i$ to a set of dependent variables $y = (y_i)_i$, one minimizes the $\chi^2$ of the model and the data set. Generally spoken, the $\chi^2$ is the relative squared difference of the model and the dependent variables. Mathematically spoken, the correlated $\chi^2$ is given by~\cite{Michael_1994,Bruno_2023}
\begin{align}
  \chi^2_{\textrm{corr}}(c) =\sum_i (y_i - f(x_i,c))\textrm{Cor}_{ij}^{-1}(y_j - f(x_j,c)).
\end{align}
Furthermore, the uncorrelated $\chi^2$ can be defined by letting out the off-diagonal elements in the covariance matrix~\cite{Bruno_2023}:
\begin{align}
  \chi^2_{\textrm{uncorr}} = (y_i - f(x_i))\Sigma_{ij}^{-1}(y_j - f(x_j))
\end{align}
Using standard optimizers one can tune the parameters $c_i$ such that the corresponding $\chi^2$ is minimized. Assuming $y_i$ is chosen from the distribution $\mathcal{N}(f(x_i,c),\sigma_{y_i})$ the $\chi^2$ is distributed as the sum of squares of standard normal distributed random variables. The resulting distribution is therefore called the $\chi^2$-distribution, which is given by
\begin{align}
\label{eq:chi2dist}
  p_n^{\chi^2}(x) = \frac{2^{-\frac{n}{2}}}{\Gamma\left( \frac{n}{2}\right)}x^{\frac{n}{2}-1}\exp\left(-\frac{x}{2} \right).
\end{align}
Here, $\Gamma$ denotes the Gamma function that generalizes the factorial to the complex plane. Within this context $n$ is often referred to the reduced number of degrees of freedom, which is given by the difference of the number independent variables $n_{\text{fit}}$ and the number of parameters $n_{\text{param}}$. Also the expectation value of $ p_n^{\chi^2}$ is given by $n$. Therefore, often the reduced chi-squared $\frac{\chi^2}{n}$ is used, whose expectation value is one. This can be used for a quick check if the fitting procedure is properly or if the model or the set of variables have to be adjusted.\\
When dealing with correlated fits another issue occurs: The covariance matrix is made from a (comparatively small) set of Jackknife samples and therefore might be ill-conditioned, which means that condition number, i.e., the ratio of smallest and largest eigenvalue of the covariance matrix is large in some sense. So there might be small modes, that are nothing but noise, which dominate the $\chi^2$ due to the inversion of the covariance matrix. Usually there are two simple ways to deal with this problem: First, one can reduce $n_{\text{fit}}$. Since the covariance matrix is a $n_{\text{fit}}\times n_{\text{fit}}$ matrix, there are less modes and the covariance matrix can be estimated properly even with a small set of samples. Of course, the removal of some of the data points has to be justified in some way. For example in a plateau fit (which can be used to measure mass states) the points that go in to the fit should be equally distributed along the fitting window.\\
Furthermore, one can also perform a regularization on the covariance matrix. Basically, this is nothing else than moving the smallest eigenvalues up to a certain bound. In detail, one has to compute the singular value decomposition (SVD) of the covariance matrix
\begin{align}
  \textrm{Cov} = U \alpha V
\end{align}
where $U$ and $V$ are unitary matrices and $\alpha$ is a real diagonal matrix, whose diagonal elements are called $(\alpha_i)_{i = 1,...,n_{\text{fit}}}$ with $\alpha_i > \alpha_{i+1}$. A minimal value for the singular values $\alpha_{\text{min}}$ is set, for example $\alpha_{\text{min}} = 10^{-8}\alpha_1$. Whenever a singular value $\alpha_i$ is smaller than $\alpha_{\text{min}}$ or smaller than twice its statistical error it is replaced by $\alpha_{i-1}$, otherwise it is simply replaced by $\vert \alpha_i \vert$. Note, that it is important to proceed in increasing order in $i$ to get fully correct results. Finally, the tweaked singular values are inverted and the singular value decomposition is inverted. Since the smallest modes are removed, noise will not dominate the $\chi^2$ of the correlated fits. A big advantage compared to the first method is that none of the data points has to be removed, which will give a more reliable result~\cite{PhysRevD.106.054502,PhysRevD.100.094508}.
\subsection{Evidence for the quality of fits}
Of course, fits can have different qualities. A set of data points with a strong curvature cannot be described properly with a straight line. On the other hand, fitting data that are within errors on a straight line with a high-order polynomial can introduce unwanted high deviations in the extrapolation regime (over-fitting). Thus, there are indeed fits with different qualities. However, some methods are needed to quantify the quality of a fit. As mentioned before the expectation value of the $\chi^2$-distribution is expected to be equal to the number of degrees of freedom, but the quantity $\chi^2/n$ is not a good parameter to describe the quality. For example this quantity tends to overvalue fits with too many parameters. But fortunately there are better alternatives.
\subsubsection{Kolmogorow-Smirnow test}
The Kolmogorow-Smirnow (KS) test is a parameter free check to test if a cumulative probability density function $F_0$ is consistent with a discrete test sample $\{X_i\}_{i = 1..n}$. To perform a KS test, an empirical distribution function
\begin{align}
  \tilde F_n(x) = \frac{1}{n} \sum_{i = 1}^n \theta(X_i - x)
\end{align}
is defined. The null-hypothesis states that $F_0$ is the correct function to describe the distribution of the $X_i$. The supremum of the differences of $F_0$ and the empirical distribution function
\begin{align}
  D_n = \sup_{x \in \mathbb{R}}\vert \tilde F_{n}(x) - F_0(x) \vert
\end{align}
is called the KS-statistic.
For a given significancy value $\alpha$, the hypothesis is accepted if the KS-statistic is smaller than
\begin{align}
  d_\alpha = \sqrt{\frac{-\frac{1}{2}\log \frac{\alpha}{2}}{n}},
\end{align} 
which is a useful approximation of the cumulated Kolmogorow distribution.\\
$d_\alpha$ also serves as the probability density function of the empirical distribution function.
Also, one-sided KS-tests can be defined in order to test if the empirical distribution function is above
\begin{align}
  D_n^u = \sup_{x \in \mathbb{R}}\vert \tilde F_{n}(x) - F_0(x) \vert\times\theta( \tilde F_{n}(x) - F_0(x))
\end{align}
or below
\begin{align}
  D_n^l = \sup_{x \in \mathbb{R}}\vert \tilde F_{n}(x) - F_0(x) \vert\times \theta(F_0(x) -\tilde F_{n}(x))
\end{align}
the distributions function from the null-hypothesis.
\subsubsection{P-value}
The P-value is a measure used in statistics to quantify the strength of evidence against a null hypothesis. It indicates how likely it is that the observed data (or even more extreme data) occurs given the null hypothesis. A low P-value means that the observed data is unlikely under the null hypothesis, which calls the null hypothesis into question. The P-value is a number between zero and one. Typically one chooses a lower bound from which the fits are trusted~\cite{Fisher1992,Wasserstein02042016}.\\
The most common example is the P-value according to the $\chi^2$-distribution of a fit. It tells how likely it is that another fit with the same degrees of freedom is worse than the fit from the null-hypothesis. In practice this means that one has to compute the value of the cumulative distribution function of the $p_n^{\chi^2}$ probability density at the $\chi^2$ of the null-hypothesis. So the P-value (for the $\chi^2$-distribution) can be computed by dividing the incomplete $\Gamma$-function with argument $\frac{n}{2}$ and $\frac{\chi^2}{2}$ divided by $\Gamma\left( \frac{n}{2}\right)$:
\begin{align}
\label{eq:pval}
  P(n,\chi^2) = \frac{\Gamma\left(\frac{n}{2},\frac{\chi^2}{2}\right)}{\Gamma\left(\frac{n}{2}\right)} = \frac{\int_{\chi^2/2}^\infty \mathrm{d}t\, t^{\frac{n}{2}-1}e^{-t} }{\int_{0}^\infty \mathrm{d}t\, t^{\frac{n}{2}-1}e^{-t} }
\end{align}
\subsubsection{Q-value}
\label{sec:qval}
Bruno and Sommer suggested not to use the P-value to estimate the quality of a fit but to replace Equation~\eqref{eq:pval} by another one~\cite{Bruno_2023}. Usually the name Q-value is also used for what I call the P-value. Within in this thesis, I will use the Bruno/Sommer definition
\begin{align}
  Q(\nu,\chi^2) = 
  \int \mathrm{d}^{n_\text{data}}y\,p(y) \theta\left( \chi^2(y) - \chi^2\right),\\
  p(y) = \frac{\exp\left(-\frac{1}{2}\chi^2(y)\right)}{\sqrt{(2\pi)^{n_\text{data}}\det C}},
\end{align}
where $\chi^2(y)$ is the (not optimized) $\chi^2$ for the data $y$ compared to the optimized fit and $C$ is the covariance matrix of the given data. Using substitutions this can be reformulated to
\begin{align}
Q =  \int \mathrm{d}^{n_\text{data}}z\, (2\pi)^{-\frac{n_\text{data}}{2}} \exp\left(-\frac{1}{2}\Vert z \Vert^2\right)\times\theta\left(z^T\cdot \nu \cdot z - \chi^2\right),
\end{align}
where $\nu$ is a symmetric positive definite matrix, that is defined by
\begin{align}
  \nu = C^{1/2}W(1-\mathcal{P})WC^{1/2},
\end{align}
with $W^2$ being either the covariance matrix (in case of correlated fits) or the matrix with the inverse squared errors on the diagonal (in case of uncorrelated fits). $\mathcal{P}$ is the projection on the $n_\text{fit}$ dimensional subspace that is spanned by the image of $W$. It was mentioned before that non-positive eigenvalues of the (measured) correlation matrix can appear, due to noise in the Jackknife estimation. So if $\lambda_i$ ($i = 1...n_\nu \leq n = n_\text{data} - n_\text{fit}$) are the positive eigenvalues of $\nu$ the Q-value can be given as~\cite{Bruno_2023}
\begin{align}
  Q(\chi^2 ,\nu) = \int\left(\prod_{i = 1}^{n_\nu}\frac{\mathrm{d}z_i\, e^{-\frac{1}{2}z_i^2}}{\sqrt{2\pi}}\right)\times\theta\left(\sum_{j = 1}^{n_\nu} \lambda_j z_j^2 - \chi^2 \right). 
\end{align}
In general it is not possible to compute the Q-value analytically but a Monte-Carlo evaluation of the integral with $\mathcal{O}(10^5)$ evaluation is precise enough to give a good estimator for the Q-value.
\subsubsection{Akaike-Information criterion (AIC)}
Usually, a lot of different fit models is used to describe the data from the simulation. Of course one could use the model with the best $P/Q$-value for the final result. But it is still possible that the best model is wrong. So a better approach is to sample a lot of different function and compute their mean value using weights that respect the quality of the fit in some way. Akaike came up with the following approach~\cite{Akaike:1974vps}: A fit through $n_\text{data}$ data points with $n_\text{fit}$ fit parameters and an optimized $\chi^2$, has a weight of
\begin{align}
  w = \mathcal{N}\exp\left(-\frac{\chi^2}{2}-n_\text{fit} + n_\text{data} \right).
\end{align}
The normalization factor $\mathcal{N}$ is chosen such that the sum of all weights is equal to one. In~\cite{Borsanyi_2021} it was suggested to modify this weight by replacing $n_\text{fit}$ by $\frac{n_\text{fit}}{2}$ so that the overall weight is given by
\begin{align}
  w = \mathcal{N}\exp\left(-\frac{\chi^2 + n_\text{fits} - 2n_\text{data}}{2}\right).
\end{align}
With this approach, fit models with many parameters are suppressed. This prevents from including models that suffer from over-fitting.

\section{Statistical and systematic covariances}
\label{sec:cov}
It is not trivial to compute systematic and statistical errors (and covariances) from the data. As mentioned in the previous sections, $N_s$ different systematics are assigned to the corresponding weights and the statistical effects are considered using $N_J$ blocked Jackknife samples and the central sample. To compute the total error and in addition separate it into its statistical and systematic contribution, the following approach was suggested in~\cite{Borsanyi_2021}: For each systematic $s = 1...N_s$, the statistical contribution is approximated by a Gaussian function with the mean value $\mu_s$ taken from the central sample and the width $\sigma_s$ calculated from the Jackknife samples. All of these Gaussian functions are added using their respective weights. Since this might be impractical for a large amount of systematics, the Gaussian functions are put in a histogram. This new total histogram again has a mean $\mu_\text{tot}$ and a width $\sigma_\text{tot}$. The variance $\sigma_\text{tot}^2$ is the sum of the variance of the systematics $\sigma_\text{sys}^2$ and the statistical variance $\sigma_\text{stat}^2$. To figure out how the total variance is divided into these two values, the procedure is repeated by using $\sqrt{\lambda}\sigma_s$ as the width for the different systematics. The resulting width of the histogram is called $\sigma_{\text{tot},\lambda}^2$. Then, the following system of linear equations can be solved:
\begin{align}
  \begin{aligned}
    \sigma_\text{sys}^2 + \sigma_\text{stat}^2 = \sigma_\text{tot}^2\\
    \sigma_\text{sys}^2 + \lambda\sigma_\text{stat}^2 = \sigma_{\text{tot},\lambda}^2
  \end{aligned},
\end{align}
The total variances and $\lambda$ are known so the system can be solved for the statistical and systematic variances~\cite{Borsanyi_2021}.\\
It is also interesting to compute statistical and systematic covariances of two observables. Of course one could fill a two dimensional histogram instead of a one dimensional one, but from experience one can know that this gets expensive in terms of computation effort because the overall number of bins can easily exceed a million if for the one dimensional case a few thousands are used. But luckily there is an easy formula to compute the covariance of two observables
\begin{align}
    \mathrm{cov}(Y_1,Y_2) = \frac{1}{2}\left(\sigma^2_{Y_1 + Y_2} - \sigma^2_{Y_1} - \sigma^2_{Y_2} \right).
\end{align}
Thus, if the variances of $Y_1$ and $Y_2$ are already known the covariance can be calculated with the same computational costs as the variances~\cite{Davier:2023cyp}. The remaining problem is the computation of the variance of $Y_1 + Y_2$. If both of the quantities come from the same fit and therefore are associated to the same weights, it is very easy. For each systematic $s$, one chooses $(Y_1 + Y_2)_s = {Y_1}_s + {Y_2}_s$ for all statistical samples and the weights $w_s$. But if both quantities have completely unrelated systematics, a different approach is needed. For example, one can draw a systematic $s_1$ for $Y_1$ according to the weights $\{w^1_s\}$ of the systematics of $Y_1$. The same is done for $Y_2$ and then ${Y_1}_{s_1} + {Y_2}_{s_2}$ is included in the histogram with the weight one. This procedure has to be repeated a certain amount of times $D$. $D$ can be estimated by setting $Y_2 = Y_1$. Then it can be checked for which value of $D$ the desired precision is reached, i.e., if the result is compatible with the  variance from the standard procedure.  If the two observables share some (flat-weighted) systematics that are aligned to an external input like different choices of a scale setting variable one can respect this correlation in the following way: Let us call these systematics $s'$ and the remaining $s_1''$ resp. $s_2''$ . Due to the flat-weighting, all systematics $s'$ have the same weight. Then it is looped over $s'$ for both observables parallel and from the other observables $s_1''$ and $s_2''$ it is chosen according to the weight. Th Gaussian functions in the histogram are weighted with one.  
\section{Generalized Eigenvalue problems}
\label{sec:GEVP}
In lattice QCD, a correlation function of an operator $O$ takes the form
\begin{align}
  C_O(t) = \sum_i \vert a^O_i\vert^2 \exp(-E_i^O t),
\end{align}
where $E_0^O < E_1^O < E_2^O...$ are the different energies contributing the operator~\cite{Varnhorst2019lecture,GattringerLang}. In many cases all states above the ground state\footnote{The ground state is the state with the lowest energy, i.e., $i = 0$.} are exponentially suppressed after a certain time $t_\text{min}$ so that for $t > t_\text{min}$ it can be assumed that the correlation function only consists of the ground state's exponential function. But this is not always the case and the following problem can appear: The coefficients of the lowest energy states can be very small which will lead to the problem that $t_\text{min}$ is very large. Therefore, it might already be in a region where the signal-to-noise ratio is too bad to get the desired precision. This can be the case if the vector-correlator, which has small contributions from some two-pion systems that can have lower energies than the mass of the vector correlator~\cite{Lahert:2023ore,Frech:2023bzy} is considered. \\
Luckily, there is a general strategy to solve this problem. To resolve $N_E$ different energy states, one needs $N_E$ operators with this energy states but contributions of different strengths. Basically this means that we look for a set of operators $\{O_j\}_{j = 1...N_E}$ that have the same quantum numbers. Now, not only the correlation functions $C_{j}(t) =\langle O_j^\dagger(t)O_j(0)\rangle$ have to be measured but also the cross-correlation functions $C_{jk}(t) =\langle O_j^\dagger(t)O_k(0)\rangle$~\cite{DeTar_2015}. Thus, the correlation function $C = \{C^{jk}\}_{j,k = 1...N_E}$ with
\begin{align}
  C^{jk}(t) = \sum_i (a^{j}_i)^*a^k_i \exp(-E_i t)
\end{align}
can be used to extract the modes $\exp(-E_it)$ using a Generalized Eigenvalue Problem (GEVP). Now a parameter $t_0$ is fixed in order to have
\begin{align}
  C^{jk}(t) &= \sum_i \underbrace{\left[(a^{j}_i)^*a^k_i\exp(-E_it_0)\right]}_{B_i = \{b^{jk}_i\}_{jk}}\times \exp(-E_i (t-t_0)).
\end{align}
Each matrix $B_l$ is constructed to be the outer product of a vector $a_l$ with itself. Therefore, it is Hermitian with rank one, and thus it has only one eigenvector $v_l$ with a non-zero eigenvalue. This eigenvector is given by $(v_l)^m = a_l^m$. It cannot be assumed that the vectors $\{v_l\}$ build an orthogonal system, but we can generally expect that $v_l^\dagger v_k$ is small compared to $v_l^\dagger v_l$ if $k \neq l$. Thus, the multiplication of $C(t)$ with $v_l$ can be approximated in the following way~\cite{Aubin_2011,Varnhorst:2020dba}:
\begin{align}
\label{eq:GEVP_aproxiamtion}
  C(t)\cdot v_l \approx B_l \cdot v_l \times \exp(-E_l(t-t_0)),
\end{align}
or, as it is written usually:
\begin{align}
  C(t)v_l &\approx \lambda_l(t,t_0)C(t_0)v_l,\\
  \lambda(t,t_0) &= \exp(-E_l(t-t_0)).
\end{align}
A GEVP can be solved by computing the eigenvalues and eigenvectors of $C(t_0)^{-1}\cdot C(t)$. To extract the functions $\lambda_l(t,t_0)$ for a fixed $t_0$ one could solve the GEVP for all $t$. However, the apparent order of the eigenvalues $\lambda_l$ may differ from the true order due to statistical fluctuations and systematic errors (from the approximation in Equation~\eqref{eq:GEVP_aproxiamtion}). Therefore, it is better to extract the eigenvectors $v_l$ for one specific choice of $t_0$ and $t$, and then compute the eigenvalues via
\begin{align}
  \lambda_l(t,t_0) = v_l^\dagger \cdot C(t) \cdot v_l.
\end{align}

