Indeed, the distribution function of the eigenphases of a random matrix in SO$(n)$ has a peak at 0 and at $\pm\pi$. It only becomes uniform for large $n$. The joint distribution function for $n=2m$ even is given by
$$p(\theta_1, \cdots, \theta_m) = C \prod_{1 \leq k < j \leq m} (\cos\theta_k - \cos\theta_j)^2~.$$
As a simple example, for $n=4$ this gives the density profile
$$p(\theta)=\frac{\cos 2 \theta+2}{2 \pi }.$$
If you sample from $O(n)$, you can either discard the matrices with determinant $-1$, or multiply the first column by $-1$.