The Veto Algorithm

The `radioactive decay' type of problems is very common, in particular in parton showers, but it is also used, e.g. in the multiple interactions description in PYTHIA. In this kind of problems there is one variable $t$, which may be thought of as giving a kind of time axis along which different events are ordered. The probability that `something will happen' (a nucleus decay, a parton branch) at time $t$ is described by a function $f(t)$, which is non-negative in the range of $t$ values to be studied. However, this naïve probability is modified by the additional requirement that something can only happen at time $t$ if it did not happen at earlier times $t' < t$. (The original nucleus cannot decay once again if it already did decay; possibly the decay products may decay in their turn, but that is another question.)

The probability that nothing has happened by time $t$ is expressed by the function ${\cal N}(t)$ and the differential probability that something happens at time $t$ by ${\cal P}(t)$. The basic equation then is

$${\cal P}(t) = - \frac{\d {\cal N}}{\d t} = f(t) \, {\cal N}(t) ~.$$ (4)

For simplicity, we shall assume that the process starts at time $t = 0$, with ${\cal N}(0) = 1$.

The above equation can be solved easily if one notes that $\d {\cal N} / {\cal N} = \d\ln {\cal N}$:

$${\cal N}(t) = {\cal N}(0) \exp \left\{ - \int_0^t f(t') \, ... ...\right\} = \exp \left\{ - \int_0^t f(t') \, \d t' \right\} ~,$$ (5)

and thus

$${\cal P}(t) = f(t) \exp \left\{ - \int_0^t f(t') \, \d t' \right\} ~.$$ (6)

With $f(t) = c$ this is nothing but the textbook formulae for radioactive decay. In particular, at small times the correct decay probability, ${\cal P}(t)$, agrees well with the input one, $f(t)$, since the exponential factor is close to unity there. At larger $t$, the exponential gives a dampening which ensures that the integral of ${\cal P}(t)$ never can exceed unity, even if the integral of $f(t)$ does. The exponential can be seen as the probability that nothing happens between the original time 0 and the final time $t$. In the parton-shower language, this corresponds to the so-called Sudakov form factor.

If $f(t)$ has a primitive function with a known inverse, it is easy to select $t$ values correctly:

$$\int_0^t {\cal P}(t') \, \d t' = {\cal N}(0) - {\cal N}(t) = 1 - \exp \left\{ - \int_0^t f(t') \, \d t' \right\} = 1 - R ~,$$ (7)

which has the solution

$$F(0) - F(t) = \ln R ~~~ \Longrightarrow ~~~ t = F^{-1}(F(0) - \ln R) ~.$$ (8)

If $f(t)$ is not sufficiently nice, one may again try to find a better function $g(t)$, with $f(t) \leq g(t)$ for all $t \geq 0$. However to use method 3 with this $g(t)$ would not work, since the method would not correctly take into account the effects of the exponential term in ${\cal P}(t)$. Instead one may use the so-called veto algorithm:

24.

start with $i = 0$ and $t_0 = 0$;

25.

add 1 to $i$ and select $t_i = G^{-1}(G(t_{i-1}) - \ln R)$, i.e. according to $g(t)$, but with the constraint that $t_i > t_{i-1}$,

26.

compare a (new) $R$ with the ratio $f(t_i)/g(t_i)$; if $f(t_i)/g(t_i) \le R$, then return to point 2 for a new try;

27.

otherwise $t_{i}$ is retained as final answer.

It may not be apparent why this works. Consider, however, the various ways in which one can select a specific time $t$. The probability that the first try works, $t = t_1$, i.e. that no intermediate $t$ values need be rejected, is given by

$${\cal P}_0 (t) = \exp \left\{ - \int_0^t g(t') \, \d t' \ri... ...t)} = f(t) \exp \left\{ - \int_0^t g(t') \, \d t' \right\} ~,$$ (9)

where the exponential times $g(t)$ comes from eq. () applied to $g$, and the ratio $f(t)/g(t)$ is the probability that $t$ is accepted. Now consider the case where one intermediate time $t_1$ is rejected and $t = t_2$ is only accepted in the second step. This gives

$${\cal P}_1 (t) = \int_0^t \d t_1 \exp \left\{ - \int_0^{t_... ..._{t_1}^t g(t') \, \d t' \right\} g(t) \, \frac{f(t)}{g(t)} ~,$$ (10)

where the first exponential times $g(t_1)$ gives the probability that $t_1$ is first selected, the square brackets the probability that $t_1$ is subsequently rejected, the following piece the probability that $t = t_2$ is selected when starting from $t_1$, and the final factor that $t$ is retained. The whole is to be integrated over all possible intermediate times $t_1$. The exponentials together give an integral over the range from 0 to $t$, just as in ${\cal P}_0$, and the factor for the final step being accepted is also the same, so therefore one finds that

$${\cal P}_1 (t) = {\cal P}_0 (t) \int_0^t \d t_1 \left[ g(t_1) - f(t_1) \right] ~.$$ (11)

This generalizes. In ${\cal P}_2$ one has to consider two intermediate times, $0 \leq t_1 \leq t_2 \leq t_3 = t$, and so

$${\cal P}_2 (t) = {\cal P}_0 (t) \int_0^t \d t_1 \left[ g(t_1) - f(t_1) \right] \int_{t_1}^t \d t_2 \left[ g(t_2) - f(t_2) \right]$$

$$= {\cal P}_0 (t) \frac{1}{2} \left( \int_0^t \left[ g(t') - f(t') \right] \d t' \right)^2 ~.$$ (12)

The last equality is most easily seen if one also considers the alternative region $0 \leq t_2 \leq t_1 \leq t$, where the rôles of $t_1$ and $t_2$ have just been interchanged, and the integral therefore has the same value as in the region considered. Adding the two regions, however, the integrals over $t_1$ and $t_2$ decouple, and become equal. In general, for ${\cal P}_i$, the $i$ intermediate times can be ordered in $i!$ different ways. Therefore the total probability to accept $t$, in any step, is

$${\cal P} (t) = \displaystyle{ \sum_{i=0}^{\infty} {\cal P}_i (t) = {\cal P}_0 (t... ...ty} \frac{1}{i!} \left( \int_0^t \left[ g(t') - f(t') \right] \d t' \right)^i }$$

$$= \displaystyle{ f(t) \exp \left\{ - \int_0^t g(t') \, \d t' \right\} \exp \left\{ \int_0^t \left[ g(t') - f(t') \right] \d t' \right\} }$$

$$= \displaystyle{ f(t) \exp \left\{ - \int_0^t f(t') \, \d t' \right\} ~, }$$ (13)

which is the desired answer.

If the process is to be stopped at some scale $t_{\mathrm{max}}$, i.e. if one would like to remain with a fraction ${\cal N}(t_{\mathrm{max}})$ of events where nothing happens at all, this is easy to include in the veto algorithm: just iterate upwards in $t$ at usual, but stop the process if no allowed branching is found before $t_{\mathrm{max}}$.

Usually $f(t)$ is a function also of additional variables $x$. The methods of the preceding section are easy to generalize if one can find a suitable function $g(t,x)$ with $f(t,x) \leq g(t,x)$. The $g(t)$ used in the veto algorithm is the integral of $g(t,x)$ over $x$. Each time a $t_i$ has been selected also an $x_i$ is picked, according to $g(t_i,x) \, dx$, and the $(t,x)$ point is accepted with probability $f(t_i,x_i)/g(t_i,x_i)$.