Note: ChocoPouce's answer is the same as mine but is more mathematical.
You have a (spherically symmetric) probability density distribution $\rho$ in space (which we get from the square of the amplitude). The "radial probability density" is roughly the chance that the electron is at a given radius, say $r = 0.1\mathrm{nm}$? In other words, how much of this distribution is in at the $0.1\mathrm{nm}$ shell? We won't exactly $0.1\mathrm{nm}$, so we take a thin shell: how much is between $0.1$ and $0.1+ε\mathrm{nm}$? Amount $=$ (average prob. density in shell)*(volume of shell). Since the shell is so thin ($\varepsilon$ is very small), the density will be almost constant and the shell's volume is given by $\mathrm{area}\times\mathrm{thickness} = 4\pi*(0.1nm^2)*\varepsilon nm$.
Thus the probability is $4\pi\rho(r)r^2\varepsilon$, where $\rho$ only depends on $r$ since it is spherically symmetric. But $\varepsilon$ is an arbitrary "small" thickness we defined, and it is best to divide by $\varepsilon$ to get the probability per unit radius, which is called the (radial) probability density $4\pi\rho(r)r^2$. The most likely radius (which is different from the average radius) maximizes this function.
To find this Fourier transform, you start with
$$\Psi(\vec{p})=\frac{1}{(2\pi \hbar)^{3/2}}\int{e^{-i\vec{p}\cdot\vec{r} / \hbar }\psi(\vec{r})\mathrm d\vec r},$$
as you've written, but you need to be a bit more careful after that.
To begin with, you can't assume that $\vec p$ is along the $z$ axis, because $\vec p$ has been given to you as the argument of $\Psi(\vec p)$, and you can't change it. What you can do, however, is set up the coordinate system for the $\vec r$ integration so that its axes simplify your life, i.e. with $z$ along the direction of $\vec p$.
(The reason you can do this, by the way, is because the ground state $\psi(\vec r)$ is spherically symmetric. If it's not symmetric, such as e.g. anything that's not an $s$ state, then you need to use a decomposition of the plane wave into partial (spherical) waves, like the one in Jackson, 2nd ed., p. 471, eq. (10.43).)
Once you've done that, then yes, you can write $\vec p = (0,0,p)$, and you can further use spherical polar coordinates $\vec r = r(\sin(\theta) \cos(\phi), \sin(\theta)\sin(\phi), \cos(\theta))$ for your position vector, write
\begin{align}
\Psi(\vec{p})
& =
\frac{1}{(2\pi \hbar)^{3/2}}
\int_0^\infty \int_0^\pi \int_0^{2\pi}
e^{-ipr\cos(\theta) / \hbar }\psi(\vec{r})
\:r^2\sin(\theta)\:\mathrm d\phi\mathrm d\theta\mathrm dr
\\& =
\frac{1}{(2\pi \hbar)^{3/2}}\sqrt{\frac{1}{\pi a_{0}^{3}}}
\int_0^\infty \int_0^\pi \int_0^{2\pi}
e^{-ipr\cos(\theta) / \hbar }e^{-r/a_0}
\:r^2\sin(\theta)\:\mathrm d\phi\mathrm d\theta\mathrm dr,
\end{align}
and integrate away.
Alternatively (and this is really for any future visitors who chance upon here looking for an actual answer for the title question), you can cheat and look up the answer, which is found in e.g.
The Momentum Distribution in Hydrogen-Like Atoms. B. Podolsky and L. Pauling. Phys. Rev., 34, 109 (1929).
(Not you, Marc, by the way. Go and do the integral.)
Best Answer
Personally, I’d do the approximation in the exponential,
$$e^\epsilon \approx 1 + \epsilon + \frac{\epsilon^2}{2!} + \cdots, $$
before doing the integral.
\begin{align} P &= \frac4{a_0^3} \int_0^R \mathrm dr\, r^2 \exp\frac{-2r}{a_0} \\ \frac{a_0^3}{4} P &\approx \int_0^R\mathrm dr\,r^2 - \frac2{a_0}\int_0^R \mathrm dr\,r^3 \\&= \frac{R^3}{3} - \frac2{a_0}\frac{R^4}{4} \\&= \frac{R^3}{3} \left( 1 - \frac32 \frac R{a_0} \right) \\ P &\approx \frac43 \left(\frac{R}{a_0}\right)^3 \end{align}
Throwing away the second term basically says “this volume is small enough we can assume the exponential is secretly a constant.”
This result suggests you are getting zeros because you’re only keeping terms that are quadratic in $R/a_0$, but apparently they all get cancelled out.