Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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...

share|cite|improve this question
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
up vote 2 down vote accepted

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}$.

share|cite|improve this answer
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
This seems to work, anyway. Thanks a lot! – Tomas Lycken Apr 12 '13 at 17:03

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.

share|cite|improve this answer
I've tried to implement that, but I can't get it to work. Part of the problem is that I don't know how to get the cumulative distribution function from my probability density, since I don't know how to normalize it when I don't know the upper limit of the energy. – Tomas Lycken Apr 12 '13 at 14:50
That is, I don't know the upper limit in advance, so I don't know how to do this analytically and insert the result in my program. – Tomas Lycken Apr 12 '13 at 14:51

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.