Drawing random lines in a cylinder - How does the distribution of angles look like?
Let $X$ denote the point in the circle below, and $Y$ denote the point in the circle above. For any given $X$, with probability $1/2$ we would have that $Y$ lies in the semi-circle whose midpoint is directly above $X$.
This means with probability $1/2$, the angle lies between $0$ and $\arccos(2/\sqrt6)\simeq 35.26^{\circ}$. Hence, the assumption that the distribution of the angles would be uniform is flawed, and about half the time the angles should lie between $35.26^{\circ}$ and $45^{\circ}$ (which is far from half of the possible angles).
Here's an attempt at calculating things more explicitly. Let $X, Y$ be uniform in the circle $C=\left\{(x,y,0)\in\mathbb R^3\,|\,x^2+y^2=1\right\}$. Your random line has direction $\big(Y+(0,0,2)\big) - X$. With $W = Y-X$, $h=(0,0,2)$ and $Z$ the cosine of the angle between that line an the $z$-axis, we have
$$Z = \frac{\langle W+h, (0,0,1)\rangle}{\lVert W+h\rVert} = \frac2{\lVert W+h\rVert}.$$
This implies
$$Z^2 = \frac{4}{4 + {\lVert W\rVert}^2} = \frac1{1+{\left\lVert\frac{W}2\right\rVert}^2}.$$
If $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$, then $X\sim \big(\cos(\theta_1),\sin(\theta_1)\big)$ and $Y\sim \big(\cos(\theta_2),\sin(\theta_2)\big)$ so that $W\sim\big(\cos(\theta_1)-\cos(\theta_2),\sin(\theta_1)-\sin(\theta_2)\big)$. Hence,
\begin{align} {\left\lVert\frac{W}2\right\rVert}^2 &\sim \frac14\left( \underbrace{(\cos\theta_1)^2+(\cos\theta_2)^2}_1-2\cos\theta_1\cos\theta_2 + \underbrace{(\sin\theta_1)^2+(\sin\theta_2)^2}_1-2\sin\theta_1\sin\theta_2 \right) \\&= \frac14\left( 2 - 2\cos\theta_1\cos\theta_2 - 2\sin\theta_1\sin\theta_2 \right) \\&= \frac12\big( 1 - (\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2) \big) \\&= \frac12\big( 1 - \cos(\theta_1-\theta_2) \big) = \frac12 - \frac{\cos(\theta_1-\theta_2)}2 . \end{align}
Therefore
$$Z^2 \sim \frac{2}{3-\cos(\theta_1-\theta_2)}.$$
Observing that $Z\geq 0$, it follows that the angle $\alpha$ is distributed as
$$\alpha \sim \arccos\left(\sqrt{\frac{2}{3-\cos(\theta_1-\theta_2)}}\right),$$
where $\theta_1,\theta_2 \sim \text{Unif}([0,2\pi])$. From here on, we can calculate $\mathbb P(\alpha \leq x)$ for $0<x<\pi/4$ as some area in the square $[0,2\pi]^2$ contained in the $\theta_1$-$\theta_2$ plane.
Alternatively, we can observe that $\Theta = \theta_1-\theta_2$ has pdf with support in $[-2\pi,2\pi]$ given by
$$f_{\Theta}(z) = \left\{\begin{array}{cccc} \frac{2\pi+z}{4\pi^2}&&&-2\pi\leq z <0\\ \frac{2\pi-z}{4\pi^2}&&&0\leq z \leq 2\pi \end{array}\right. $$
Then, some algebraic manipulations show that
$$\mathbb P(\alpha \leq x) = \mathbb P\left(\cos(\Theta) \geq 3-\frac{2}{{\cos(x)}^2}\right) .$$
Using the pdf above, we can calculate this explicity in terms of integrals that simplify to
$$\mathbb P(\alpha \leq x) = \frac{1}{\pi}\,\arccos\left(3-\frac{2}{{\cos(x)}^2}\right) .$$
The density function $f_\alpha$ of $\alpha$ has support in $[0,\pi/4]$, and according to Wolfram we have (after simplifications considering its support) that it equals
$$f_{\alpha}(z) = \frac{2}{\pi}\,\frac{1}{\cos(z)\,\sqrt{\cos(2z)}},$$
which is way simpler than I expected, to be honest. Observe that it's indeed a probability density function, in that it is non-negative and its integral is $1$.
Your histogram does somewhat resemble the area under $f_{\alpha}$, as you can verify here, so this looks more like a case of 'wrong intuition' than of 'wrong code'.
I must point out that this sheds light to yet another rather curious identity involving $\pi$:
$$\int_{\sqrt{2}/2}^1\,\,\,\frac{1}{u\,\sqrt{2u^2-1}\,\sqrt{1-u^2}}\,\,\,du = \frac\pi2$$ $${}$$ $$\int_{1/2}^1\,\,\,\frac{1}{u\,\sqrt{2u-1}\,\sqrt{1-u}}\,\,\,du = \pi$$
It's possible to find the closed form for the cdf. The angle is $$A = \arccos \frac {((\cos v, \sin v, 1) - (\cos u, \sin u, -1)) \cdot (0, 0, 1)} {| (\cos v, \sin v, 1) - (\cos u, \sin u, -1) |} = \arccos \sqrt {\frac 2 {3 - \cos(v - u)}}.$$ Switching to $\xi = v + u, \eta = v - u$, the cdf is $$\operatorname P(A < x) = \int_0^{2\pi} \int_\eta^{4\pi - \eta} \frac 1 {(2 \pi)^2} [A < x] d\xi d\eta \overset * = \\ \int_0^\pi \frac 2 {(2 \pi)^2} (4\pi - 2\pi) [A < x] d\eta = \\ \frac 1 \pi \int_0^\pi [\eta \lt \arccos(3 - 2 \sec^2 x)] d\eta = \\ \frac 1 \pi \arccos(3 - 2 \sec^2 x),$$ where $[\text {condition}]$ is $1$ if the condition is true and $0$ otherwise. For $*$, $f$ is symmetric around $\pi$, therefore $\int_0^{2\pi} \eta f(\eta) d\eta = 2 \pi \int_0^\pi f(\eta) d\eta$.