How to set up phase space of the proton beam?

Closing as a duplicate of How to generate random numbers according to the 2D probability density function