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…
Best Answer
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}$.