I am not an expert but here is my understanding. It is best to first understand [MP12]. A lot of this is to give background and not exactly on-point, also for my reference.
Question: Why does SampleG return a spherical distribution?
Answer: SampleG returns a spherical distribution because the authors are using the convolution technique, introduced in [Pei10], Theorem 3.1, and concretely used in MP12, Section 5.4, point (2). Let's look at MP12 first. The theorem statement is as follows, where $\Sigma_*$ are positive definite matrices and $\eta_{\varepsilon}(\Lambda) = O(\sqrt{\ln(n / \varepsilon)})$ is the smoothing parameter.
Theorem ([Pei10], Thm. 3.1): Let $\color{red}{\Sigma = \Sigma_1 + \Sigma_2}$ and lattices $\Lambda_1, \Lambda_2$ such that $\sqrt{\Sigma_1} \geq \eta_{\varepsilon}(\Lambda_1)$ and $\sqrt{\left(\Sigma_1^{-1} + \Sigma_2^{-1}\right)^{-1}} \geq \eta_{\varepsilon}(\Lambda_2)$. Consider the distribution of the output of the experiment
\begin{align*}
\mathbf{x}_2 &\leftarrow \mathcal{D}_{\Lambda_2 + \mathbf{c}_2, \color{red}{\sqrt{\Sigma_2}}} \\
\mathbf{x}_1 &\leftarrow \mathbf{x}_2 + \mathcal{D}_{\Lambda_1 + \mathbf{c}_1 - \mathbf{x}_2, \color{red}{\sqrt{\Sigma_1}}} \\
&\mathrm{Output} \ \ \mathbf{x}_1
\end{align*}
Then, the (marginal) distribution of $\mathbf{x}_1$ is within statistical distance of $8\varepsilon$ of $\mathcal{D}_{\Lambda_1 + \mathbf{c}_1, \color{red}{\sqrt{\Sigma}}}$.
This theorem is significant in two ways, highlighted in $\color{red}{\text{red}}$. For one, it guarantees the distributions of summing two Gaussian variables. Looking the other direction, it also means you can split sampling from one distribution into sampling from two other distributions, if they are more convenient.
It is exactly what's done in MP12 and GM18: They are trying to sample preimages from $\mathcal{D}_{\Lambda^{\perp}_{\mathbf{u}}(\mathbf{G}), s}$ with covariance $\Sigma = s\mathbf{I}_n$, for which there is a straightforward way that outputs vectors with covariance $\Sigma_1 = \begin{bmatrix} \mathbf{R} \\ \mathbf{I} \end{bmatrix} \Sigma_{\mathbf{G}} \begin{bmatrix} \mathbf{R}^T & \mathbf{I} \end{bmatrix}$, so if we use the convolution theorem above, we can simply sample a perturbation $\mathbf{p}$ from a distribution with $\Sigma_2 = \Sigma - \Sigma_1$, then sample from the preimage of $\mathbf{u - Ap}$ instead of $\mathbf{u}$. More concretely, they sample directly from $\mathcal{D}_{\mathcal{L}(\Sigma_2), s\sqrt{\Sigma_2}}$, which there is an extremely simple sampler (Klein's, see [[MP18]], Figure 1) because the "bases are aligned"!
To show it even more clearly, here is an alternative presentation of SampleG:
\begin{align*}
\sigma &:= s / r \\
\mathbf{p} &\leftarrow \operatorname{Perturb}(\sigma) \\
\mathbf{c} &:= \mathbf{B}_{b^k}^{-1}(\mathbf{u} - \mathbf{p}) \\
\mathbf{z} &\leftarrow \operatorname{SampleD}(\sigma, \mathbf{c}) \\
\mathbf{t} &:= \mathbf{B}_q\mathbf{z} + \mathbf{u} \\
&\text{Output} \ \mathbf{t}
\end{align*}
Where Perturb and SampleD are algorithms such that $\mathbf{p} \sim \mathcal{D}_{\mathcal{L}(\Sigma_2), \sigma\sqrt{\Sigma_2}}$ and $D\mathbf{z} \sim \mathcal{D}_{\mathcal{L}(D), -\mathbf{c}, \sigma}$. To simplify it even more, note that
$$
\mathbf{B}_q\mathbf{z} + \mathbf{u} = \mathbf{B}_{b^k}(D\mathbf{z} + \mathbf{B}_{b^k}^{-1}\mathbf{u}) = \mathbf{B}_{b^k}(D\mathbf{z} + \mathbf{c}) + \mathbf{p}
$$
so the algorithm can be written as
\begin{align*}
\sigma &:= s / r \\
\mathbf{p} &\leftarrow \mathcal{D}_{\mathcal{L}(\Sigma_2), \sigma\sqrt{\Sigma_2}} \\
\mathbf{z}' &\leftarrow \mathcal{D}_{\mathcal{L}(\mathbf{B}_q) + \mathbf{u} - \mathbf{p}, \color{blue}{\sigma\mathbf{B}_{b^k}}} \\
\mathbf{t} &:= \mathbf{p} + \mathbf{z}' \\
&\text{Output} \ \mathbf{t}
\end{align*}
where, of course, $\color{blue}{\sigma\mathbf{B}_{b^k}} = \sigma\sqrt{\Sigma_1}$
Question: What is the significance of $\beta$ in Perturb?
Answer: It is an "unrolled" version of the randomised nearest plane algorithm, described in [GPV08] section 4.2, combined with the compact upper triangular representation of $\sqrt{\Sigma_2}$ in ([GM18]) section 3.1. You can derive it yourself by unwrapping the dot products in the [GPV08] description. This is exactly why having a simple (sparse) $\sqrt{\Sigma}$ is important - it allows for efficient sampling!
Question: Why is the lattice $\mathcal{L}(\Sigma_2)$ in Perturb?
Answer: I think I answered above. It can be chosen arbitrarily, since you can just adjust the other matrix in the convolution theorem above. Also, sampling for $\mathcal{D}_{\mathcal{L}(\Sigma), \sqrt{\Sigma}}$ is extremely quick.
I highly recommend first understanding [MP12]'s framework especially with perturbation sampling (section 5.2). The difference between [MP12] and [GM18] essentially lies in finding an efficient factorisation of $\sqrt{\Sigma_2}$, but the sampling is essentially the same.
Best,
