Integrating over a hypercube, not a hypersphere Denote $\square_m=\{\pmb{x}=(x_1,\dots,x_m)\in\mathbb{R}^m: 0\leq x_i\leq1,\,\,\forall i\}$ be an $m$-dimensional cube.
It is all too familiar that $\int_{\square_1}\frac{dx}{1+x^2}=\frac{\pi}4$.

QUESTION. If $\Vert\cdot\Vert$stands for the Euclidean norm, then is this true?
  $$\int_{\square_{2n-1}}\frac{d\pmb{x}}{(1+\Vert\pmb{x}\Vert^2)^n}
=\frac{\pi^n}{4^nn!}.$$

 A: Yes, this is true, it follows from formulas in Higher-Dimensional Box Integrals, by Jonathan M. Borwein, O-Yeat Chan, and R. E. Crandall. 
\begin{align}
&\text{define}\;\;C_{m}(s)=\int_{[0,1]^m}(1+|\vec{r}|^2)^{s/2}\,d\vec{r},\;\;\text{we need}\;\;C_{2n-1}(-2n).\\
&\text{the box integral is}\;\;B_m(s)=\int_{[0,1]^m}|\vec{r}|^s\,d\vec{r},\\
&\text{related to our integral by}\;\;2n C_{2n-1}(-2n)=\lim_{\epsilon\rightarrow 0}\epsilon B_{2n}(-2n+\epsilon)\equiv {\rm Res}_{2n}.
\end{align}
The box integral $B_n(s)$ has a pole at $n=-s$, with "residue" ${\rm Res}_n$. Equation 3.1 in the cited paper shows that the residue is a piece (an $n$-orthant) of the surface area of the unit $n$- sphere, so it is readily evaluated,
\begin{align}
&{\rm Res}_{n}=\frac{1}{2^{n-1}}\frac{\pi^{n/2}}{\Gamma(n/2)},\\
&\text{hence}\;\;C_{2n-1}(-2n)=\frac{1}{2n}\frac{1}{2^{2n-1}}\frac{\pi^{n}}{\Gamma(n)}=\frac{1}{4^n}\frac{\pi^n}{n!}\;\;\text{as in the OP}.
\end{align}
More generally, we can consider
$$C_{p-1}(-p)=\int_{[0,1]^{p-1}}(1+|\vec{r}|^2)^{-p/2}\,d\vec{r}=\frac{1}{p}\,{\rm Res}_p=\frac{1}{p} \frac{1}{2^{p-1}}\frac{\pi^{p/2}}{\Gamma(p/2)}. $$
This formula holds for even and odd $p\geq 2$, as surmised by Gro-Tsen in a comment.

There are similarly remarkable hypercube integrals where this came from, for example
$$\int_{[0,\pi/4]^k}\frac{d\theta_1 d\theta_2\cdots d\theta_k}{(1+1/\cos^2\theta_1+1/\cos^2\theta_2\cdots+1/\cos^2\theta_k)^{1/2}}=\frac{k!^2\pi^k}{(2k+1)!}.$$
A: Here is an elementary proof. Using the formula 
\begin{equation}
 \frac1{a^n}=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u\,a}\,du \tag{1}
\end{equation}
with $a=1+\|x\|^2$, 
denoting your integral by $J_n$, letting $I:=[0,1]$, $\phi(z):=\frac1{\sqrt{2\pi}}\,e^{-z^2/2}$, and $\Phi(z):=\int_{-\infty}^z\phi(t)\,dt$, and making substitutions $x_1=z_1/\sqrt{2u}$ and then $u=z^2/2$, we have 
\begin{align}
 J_n&=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u}\,du \int_{I^{2n-1}}e^{-u\,\|x\|^2}\,dx \\ 
 &=\frac1{\Gamma(n)}\,\int_0^\infty u^{n-1} e^{-u}\,du \Big(\int_I e^{-u\,x_1^2}\,dx_1\Big)^{2n-1} \\ 
 &=\frac{\pi^{n-1/2}}{\Gamma(n)}\,\int_0^\infty u^{-1/2} e^{-u}\,du \Big(\Phi(\sqrt{2u})-\frac12\Big)^{2n-1} \\ 
 &=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,\phi(z)\,dz \\ 
 &=\frac{2\pi^{n}}{\Gamma(n)}\,\int_0^\infty \Big(\Phi(z)-\frac12\Big)^{2n-1}\,d\Phi(z) \\ 
 &=\frac{\pi^{n}}{4^n \Gamma(n+1)}, 
\end{align}
as desired. 
This derivation obviously holds whenever $2n-1$ is a positive integer. 

Further comments:
More generally, in view of Bernstein's theorem on completely monotone functions, one can quite similarly express the box integral 
\begin{equation}
 \int_{I^n} g(\|x\|^2)\,dx
\end{equation}
as an ordinary integral over $[0,\infty)$ -- for any completely monotone function $g$. Indeed, Bernstein's theorem states that any completely monotone function is a positive mixture of decreasing exponential functions; identity (1) is then such an explicit Bernstein representation of the completely monotone function $a\mapsto \frac1{a^n}$. 
Now, furthermore, one does not have to insist that the mixture be positive, since we are dealing here with identities rather than inequalities. So, any function $g$ that can be represented as the mixture 
\begin{equation}
 g(a)=\int_0^\infty e^{-u\,a}\,\mu(du)
\end{equation}
of decreasing exponential functions $a\mapsto e^{-u\,a}$ 
for $a>0$, where $\mu$ is any finite, possibly signed ("mixing") measure, will do just as well. 
