I am not expecting a "realistic" answer to my question, since it is based on an impossible scenario. What I'm waiting for is a theoretic explanation/interpretation so that I can sleep at night :)
Let's take n reals from the interval (0; 1), evenly distributed. What is the probability of the sum of squares being less than 1? With a geometric approach it is relatively easy to see that the solution is $$P\left(\sum_{i=1}^nx_i^2<1\right)={\pi^{\frac n 2}\over 2^n\cdot{\frac n 2}!}=\frac {volume \;of\;hypersphere}{volume\;of\;hypercube}$$ This solution works fine for all nonnegative integer n (for odd n we compute the factorial using the $\Gamma$ function).
But what if n could be something else, namely a real from (0; 1)? In that case the resulting probability is greater than 1, (properly) indicating that something went wrong. Is there a way to make sense of this result?
The result can make sense, but I think the probability part is a red herring. You have a function which is relevant to your probability function for integer $n$: \begin{equation}f\left(n\right)=\frac{\pi^{\frac{n}{2}}}{2^{n}\cdot\frac{n}{2}!}\end{equation} where $\frac{n}{2}!$ is defined for odd $n$ based on $\frac{1}{2}!=\frac{\sqrt{\pi}}{2}$. If you choose to extend this with the $\Gamma$ function, then $f$ goes above $1$.
However, the important point to note is that this isn't a completely arbitrary choice/occurrence. If you want to extend $f$ to a nice (i.e. analytic) function on $\mathbb{R}$ (or some big interval thereof), then since $f(0)=f(1)=1$, and it can't be constant on $[0,1]$, it's forced to go above $1$ somewhere.
As an aside, if you want an aesthetic reason to use the Gamma function in particular, the Bohr–Mollerup theorem says that the Gamma function (on the positive reals) is the only function $g$ with $g(1)=1$ (corresponding to $0!=1$), obeys the functional identity $g(x+1)=xg(x)$ (which applies even for half-integer factorials), and is "log-convex" (so that $\log g(x)$ is convex).