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}{4 \pi }.$$ If you sample from $O(n)$, you can either discard the matrices with determinant $-1$, or multiply the first column by $-1$.
By way of illustration, this is the density profile for $n=100$ (computed by averaging over 1000 random matrices), when the peaks at $0$ and $\pm\pi$ have become very small.
