# Generating a random number from a given distribution

I have a problem (a part of a Monte Carlo simulation) where I'm given the energy of an incoming particle, $\varepsilon$ and want to split this energy in two parts, randomly generating the fraction that goes into each part according to a given function

$P(\varepsilon_s) \propto %\frac{1}{\bar E \arctan\left(\frac{\varepsilon}{\bar E}\right)} \frac{1}{1+\left(\frac{\varepsilon_s}{\bar E}\right)^2}$

$\bar E$ is a known shape parameter. Now, the incoming particle can have energies ranging from 0 to a known maximum, $\varepsilon \in \left[0,\varepsilon_{\text{max}}\right]$, but I do not know the energy until the collision event, when I want to split it up, and therefore I can't figure out how to normalize the distribution function properly.

I have a pseudo-random number generator at my disposal which generates uniformly distributed numbers on $[0,1)$. How do I generate random numbers that properly describe how the two collisional products share the available energy?

I've tried reading up on methods for this in available books, but they're all a little to theoretical for me to be able to apply them to my scenario. It's been way too long since I took a proper statistics course, and now it seems I need it :)

Update: As requested, a more thorough description of $P$ (which will get into detail about the physics of this experiment...) follows here:

The physics of this problem is ionization - I have an incoming electron with kinetic energy $\varepsilon$, that collides with a molecule and ionizes it. After the collision event, I have two electrons with energies $\varepsilon_p$ (primary) and $\varepsilon_s$ (secondary), with $\varepsilon_p+\varepsilon_s=\varepsilon-\varepsilon_i$, the latter being the (known) ionization energy. Now, I don't really care which electron is which after the collision, so determining the two energies is simply a problem of dividing $\varepsilon-\varepsilon_i$ into two parts.

To do this, I am using the differential cross section of the collision, as suggested by Y. Tzeng and E. E. Kunhardt1, which is given by them as

$q(\varepsilon, \varepsilon_s) = \frac{C(\varepsilon)}{\bar E \arctan\left((\varepsilon-\varepsilon_i)/\bar E\right)}\cdot \frac{1}{1+\left({\varepsilon_s/\bar E}\right)^2}$

$C(\varepsilon)$ in this epxression is the total collision cross section, which is also known. I'm not 100% sure on this, but if I've understood the physics correctly then normalizing this function over the available energies should give me a probability density. As it turn out, this funciton is normalized so that

$\int_0^{\varepsilon-\varepsilon_i} q(\varepsilon, \varepsilon_s)d\varepsilon_s = C(\varepsilon)$

so normalizing w.r.t. the entire interval is just a matter of removing the factor $C(\varepsilon)$. However, if I insert $\varepsilon=\varepsilon_{\text{max}}$, that gives me a probability distribution which allows for energies higher than the available energy, and if I don't then I don't know how to generate random numbers from the distribution correctly.

1 I don't think I'm at liberty to relay the entire paper here, unfortunately...

-
Can you describe $P$ a bit more? –  oks Apr 12 '13 at 14:56
@Learner: $\varepsilon_{\text{max}}$ is known, but not the energy $\varepsilon$ of the incoming particle, which is an upper bound on the secondary. –  Tomas Lycken Apr 12 '13 at 15:11
@oks: I updated the question with the physics background. Please let me know if things are still unclear... –  Tomas Lycken Apr 12 '13 at 15:27

If my calculations are correct, the normalizing constant is $c \left( \varepsilon_{\max}, \bar{E} \right) = \bar{E} \arctan \left( \frac{\varepsilon_{\max}}{\bar{E}} \right)$, the CDF is $F \left( \varepsilon_s \right) = \frac{\arctan \left( \frac{\varepsilon_s}{\bar{E}} \right)}{\arctan \left( \frac{\varepsilon_{\max}}{\bar{E}} \right)}$ and thus the inverse CDF is $$\varepsilon_s = F^{- 1} \left( u \right) = \bar{E} \tan \left( \arctan \left( \frac{\varepsilon_{\max}}{\bar{E}} \right) u \right)$$ Thus draw a uniform number on $[0,1]$, set it to $u$ and compute the random draw $\varepsilon_s$ from it using the function $F^{-1}$.
Thanks a lot! Is $\varepsilon_{\text{max}}$ here the maximum possible energy available, or the energy available in this reaction? –  Tomas Lycken Apr 12 '13 at 15:39
If you can invert the cumulative distribution function for your function $P$, then you can use Inverse transform sampling. The first 2 paragraphs of the wiki article give you the idea.