Distribution of $\sqrt{X^2+Y^2}$ when $X,Y$ are independent $U(0,1)$ variables 
As a routine exercise, I am trying to find the distribution of $\sqrt{X^2+Y^2}$ where  $X$ and $Y$ are independent $ U(0,1)$ random variables.

The joint density of $(X,Y)$ is  $$f_{X,Y}(x,y)=\mathbf 1_{0<x,y<1}$$
Transforming to polar coordinates $(X,Y)\to(Z,\Theta)$ such that $$X=Z\cos\Theta\qquad\text{ and }\qquad Y=Z\sin\Theta$$ 
So, $z=\sqrt{x^2+y^2}$ and $0< x,y<1\implies0< z<\sqrt 2$.
When $0< z<1$, we have $0< \cos\theta<1,\,0<\sin\theta<1$ so that $0<\theta<\frac{\pi}{2}$.
When $1< z<\sqrt 2$, we have $z\cos\theta<\implies\theta>\cos^{-1}\left(\frac{1}{z}\right)$, as $\cos\theta$ is decreasing on $\theta\in\left[0,\frac{\pi}{2}\right]$; and $z\sin\theta<1\implies\theta<\sin^{-1}\left(\frac{1}{z}\right)$, as $\sin\theta$ is increasing on $\theta\in\left[0,\frac{\pi}{2}\right]$.
So, for $1< z<\sqrt 2$, we have $\cos^{-1}\left(\frac{1}{z}\right)<\theta<\sin^{-1}\left(\frac{1}{z}\right)$.
The absolute value of jacobian of transformation is $$|J|=z$$
Thus the joint density of $(Z,\Theta)$ is given by 
$$f_{Z,\Theta}(z,\theta)=z\mathbf 1_{\{z\in(0,1),\,\theta\in\left(0,\pi/2\right)\}\bigcup\{z\in(1,\sqrt2),\,\theta\in\left(\cos^{-1}\left(1/z\right),\sin^{-1}\left(1/z\right)\right)\}}$$
Integrating out $\theta$, we obtain the pdf of $Z$ as
$$f_Z(z)=\frac{\pi z}{2}\mathbf 1_{0<z<1}+\left(\frac{\pi z}{2}-2z\cos^{-1}\left(\frac{1}{z}\right)\right)\mathbf 1_{1<z<\sqrt 2}$$

Is my reasoning above correct? In any case, I would like to avoid this method and instead try to find the cdf of $Z$ directly. But I couldn't find the desired areas while evaluating $\mathrm{Pr}(Y\le \sqrt{z^2-X^2})$ geometrically.

EDIT.
I tried finding the distribution function of $Z$ as
\begin{align}
F_Z(z)&=\Pr(Z\le z)
\\&=\Pr(X^2+Y^2\le z^2)
\\&=\iint_{x^2+y^2\le z^2}\mathbf1_{0<x,y<1}\,\mathrm{d}x\,\mathrm{d}y
\end{align}
Mathematica says this should reduce to
$$F_Z(z)=\begin{cases}0 &,\text{ if }z<0\\ \frac{\pi z^2}{4} &,\text{ if } 0< z<1\\ \sqrt{z^2-1}+\frac{z^2}{2}\left(\sin^{-1}\left(\frac{1}{z}\right)-\sin^{-1}\left(\frac{\sqrt{z^2-1}}{z}\right)\right) &,\text{ if }1< z<\sqrt 2\\ 1 &,\text{ if }z>\sqrt 2 \end{cases}$$
which looks like the correct expression. Differentiating $F_Z$ for the case $1< z<\sqrt 2$ though brings up an expression which doesn't readily simplify to the pdf I already obtained.
Finally, I think I have the correct pictures for the CDF:
For $0<z<1$ :
 
And for $1<z<\sqrt 2$ :

Shaded portions are supposed to indicate the area of the region $$\left\{(x,y):0<x,y< 1\,,\,x^2+y^2\le z^2\right\}$$
The picture immediately yields 
\begin{align}
F_Z(z)&=\Pr\left(-\sqrt{z^2-X^2}\le Y\le\sqrt{z^2-X^2}\right)
\\&=\begin{cases}\frac{\pi z^2}{4} &,\text{ if } 0<z<1\\\\ \sqrt{z^2-1}+\int_{\sqrt{z^2-1}}^1 \sqrt{z^2-x^2}\,\mathrm{d}x &,\text{ if }1< z<\sqrt 2 \end{cases}
\end{align}
, as I had previously found.
 A: 
$f_z(z)$ :

So, for $1\le z<\sqrt 2$, we have
  $\cos^{-1}\left(\frac{1}{z}\right)\le\theta\le\sin^{-1}\left(\frac{1}{z}\right)$

You can simplify your expressions when you use symmetry and evaluate the expressions for $\theta_{min} < \theta < \frac{\pi}{4}$. Thus, for half of the space and then double the result.
Then you get:
$$P(Z \leq r) = 2 \int_0^r z \left(\int_{\theta_{min}}^{\frac{\pi}{4}}d\theta\right) dz = \int_0^r z \left(\frac{\pi}{2}-2\theta_{min}\right)  dz$$ 
and your $f_z(z)$ is 
$$f_z(z) =  z \left(\frac{\pi}{2}-2\theta_{min}\right) = \begin{cases} z\left(\frac{\pi}{2}\right) & \text{ if } 0 \leq z \leq 1 \\ z \left(\frac{\pi}{2} - 2 \cos^{-1}\left(\frac{1}{z}\right)\right) & \text{ if } 1 < z \leq \sqrt{2} \end{cases}$$

$F_z(z)$ :
You can use the indefinite integral:
$$\int z \cos^{-1}\left(\frac{1}{z}\right) = \frac{1}{2} z \left( z \cos^{-1}\left(\frac{1}{z}\right) - \sqrt{1-\frac{1}{z^2}} \right) + C $$
note $\frac{d}{du} \cos^{-1}(u) = - (1-u^2)^{-0.5}$
This leads straightforward to something similar as Xi'ans expression for $Pr(Z \leq z)$ namely 
if $1 \leq z \leq \sqrt{2}$ then:
$$F_z(z) = {z^2} \left(\frac{\pi}{4}-\cos^{-1}\left(\frac{1}{z}\right) + z^{-1}\sqrt{1-\frac{1}{z^2}} \right)$$

The relation with your expression is seen when we split up the $cos^{-1}$ into two $cos^{-1}$ expressions, and then converted to different $sin^{-1}$ expressions.
for $z>1$ we have
$$\cos^{-1}\left(\frac{1}{z}\right) = \sin^{-1}\left(\sqrt{1-\frac{1}{z^2}}\right) = \sin^{-1}\left(\frac{\sqrt{z^2-1}}{z}\right) $$
and
$$\cos^{-1}\left(\frac{1}{z}\right) = \frac{\pi}{2} -\sin^{-1}\left(\frac{1}{z}\right) $$
so
$$\begin{array}\\
\cos^{-1}\left(\frac{1}{z}\right) & = 0.5 \cos^{-1}\left(\frac{1}{z}\right)  + 0.5 \cos^{-1}\left(\frac{1}{z}\right) \\ 
 & = \frac{\pi}{4} - 0.5 \sin^{-1}\left(\frac{1}{z}\right)  + 0.5 \sin^{-1}\left(\frac{\sqrt{z^2-1}}{z}\right) \end{array} $$
which results in your expression when you plug this into the before mentioned $F_z(z)$ for $1<z<\sqrt{2}$
A: That the pdf is correct can be checked by a simple simulation
samps=sqrt(runif(1e5)^2+runif(1e5)^2)
hist(samps,prob=TRUE,nclass=143,col="wheat")
df=function(x){pi*x/2-2*x*(x>1)*acos(1/(x+(1-x)*(x<1)))}
curve(df,add=TRUE,col="sienna",lwd=3)


Finding the cdf without the polar change of variables goes through
\begin{align*}
\mathrm{Pr}(\sqrt{X^2+Y^2}\le z) &= \mathrm{Pr}(X^2+Y^2\le z^2)\\
&= \mathrm{Pr}(Y^2\le z^2-X^2)\\
&=\mathrm{Pr}(Y\le \sqrt{z^2-X^2}\,,X\le z)\\
&=\mathbb{E}^X[\sqrt{z^2-X^2}\mathbb{I}_{[0,\min(1,z)]}(X)]\\
&=\int_0^{\min(1,z)} \sqrt{z^2-x^2}\,\text{d}x\\
&=z^2\int_0^{\min(1,z^{-1})} \sqrt{1-y^2}\,\text{d}y\qquad [x=yz\,,\ \text{d}x=z\text{d}y]\\
&=z^2\int_0^{\min(\pi/2,\cos^{-1} z^{-1})} \sin^2{\theta} \,\text{d}\theta\qquad [y=\cos(\theta)\,,\ \text{d}y=\sin(\theta)\text{d}\theta]\\
&=\frac{z^2}{2}\left[ \min(\pi/2,\cos^{-1} z^{-1}) - \sin\{\min(\pi/2,\cos^{-1} z^{-1})\}\cos\{\min(\pi/2,\cos^{-1} z^{-1}\}\right]\\
&=\frac{z^2}{2}\begin{cases}
\pi/2 &\text{ if }z<1\\
\cos^{-1} z^{-1}-\sin\{\cos^{-1} z^{-1})\}z^{-1}&\text{ if }z\ge 1\\
\end{cases}\\
&=\frac{z^2}{2}\begin{cases}
\pi/2 &\text{ if }z<1\\
\cos^{-1} z^{-1}-\sqrt{1-z^{-2}}z^{-1}&\text{ if }z\ge 1\\
\end{cases}
\end{align*}
which ends up with the same complexity! (Plus potential mistakes of mine along the way!)
A: For $0 \leq z \leq 1$, $P\left(\sqrt{X^2+Y^2} \leq z\right)$ is just the area of the quarter-circle of radius $z$ which is $\frac 14 \pi z^2$.  That is, 
$$\text{For }0 \leq z \leq 1, ~\text{area of quarter-circle} = \frac{\pi z^2}{4} = P\left(\sqrt{X^2+Y^2} \leq z\right).$$
For $1 < z \leq \sqrt{2}$, the region over which we need to integrate to find $P\left(\sqrt{X^2+Y^2} \leq z\right)$can be divided into two right triangles $\big($one of them has vertices $(0,0), (0,1)$ and $(\sqrt{z^2-1}, 1)$ while the other has vertices $(0,0), (1,0)$ and $(1, \sqrt{z^2-1})$ $\big)$ together with a sector of a circle of radius $z$ and included angle $\frac{\pi}{2}-2\arccos\left(\frac{1}{z}\right)$. The area of this region (and hence the value of $\left( P(\sqrt{X^2+Y^2} \leq z\right)$) is easily found. We have that for $1 < z \leq \sqrt{2}$,
\begin{align}\text{area of region} &= \text{area of two triangles plus area of sector}\\
&=\sqrt{z^2-1} + \frac 12 z^2\left( \frac{\pi}{2}-2\arccos \left(\frac{1}{z}\right)\right)\\
&= \frac{\pi z^2}{4} + \sqrt{z^2-1} - z^2\arccos \frac{1}{z}\\
&= \left( P(\sqrt{X^2+Y^2} \leq z\right)\end{align}
which is the result in Martijn Wetering's answer.
