As described here, I'm a developer of the TOPAS Monte Carlo simulation tool for radiotherapy research. To improve our beam source model, I wanted a method to randomly generate uniformly distributed directions that are constrained to a certain region of the unit sphere, centered upon the z-axis (the beam direction). Generated particles are assigned these random directions before they are transported through the simulated geometry, in order to account for beam divergence.
I've written a Jupyter Notebook describing how I tackled this issue. It can either be viewed statically below or interactively in binder. The interactive version provides a widget to adjust beam model parameters and discover how the sampling changes, and also visualizes 3D graphics that you can spin with your cursor.