KAHIBARO
28.7. Number of Simulated Events
Table of Contents
Convergence
When you run a Geant4 simulation, each event is one random realization of what can happen. The result you care about, for example a dose in a voxel or a peak position in an energy spectrum, is an average over many events. As you increase the number of simulated events, this average converges to the true underlying value predicted by the physics in your model.
Convergence is the practical question: "Have I simulated enough events so that my result will not change in any important way if I keep running more?" A converged observable changes only within its expected statistical fluctuations as you continue to add events.
For a simple estimator like the mean of some quantity $x$ per event, computed as
$$
\bar{x}_N = \frac{1}{N}\sum_{i=1}^{N} x_i,
$$
convergence means that as $N$ grows, $\bar{x}_N$ approaches a stable value. In practice, you do not know the exact true value, so you look at how $\bar{x}_N$ behaves as you increase $N$.
A simple workflow to assess convergence is to run the simulation in stages, for example in blocks of $N_\text{block}$ events. After each block, you record the current estimate of your quantity of interest. If the estimate keeps shifting by amounts larger than your acceptable error, you have not yet converged. If changes become small and consistent with the expected statistical variation, you can consider the result converged for your application.
Convergence is always relative to a target precision. A medical physics dose calculation might need a few percent accuracy in a patient voxel, while a detector design study might tolerate 10 percent for a preliminary comparison. The same simulation could be considered converged for one goal and unconverged for another, depending on the required precision.
You will often see that some observables converge faster than others. For example, a total absorbed dose in a large volume converges faster than a finely binned energy spectrum. This happens because each bin in a histogram is effectively estimated from a smaller subset of events. If you care about the fine structure in a spectrum, you must consider convergence at the bin level, not only for the overall integral.
A practical approach is to plot your key observable as a function of the total number of events simulated. If the curve levels off and only fluctuates within your required tolerance, you have good evidence of convergence. If it continues to drift or shows clear trends as $N$ increases, you need more events or you may be seeing a systematic effect rather than a statistical one.
A simulation is only useful if the observables of interest have converged within the precision you require. Always study results as a function of the number of events before trusting them.
Statistical precision
The precision of a Monte Carlo result is governed by statistical uncertainties due to the finite number of events. For many quantities that are averages of independent events, the standard uncertainty scales approximately as
$$
\sigma \propto \frac{1}{\sqrt{N}},
$$
where $N$ is the number of simulated events.
For independent events, doubling the number of events improves the statistical uncertainty only by a factor of $\sqrt{2}$, and to reduce the uncertainty by a factor of 10 you need 100 times more events.
Consider a simple count, such as the number of photons detected in a given detector element. If the expected number of counts is $C$, and events are independent and rare enough, the distribution is approximately Poisson with variance $C$ and standard deviation
$$
\sigma_C \approx \sqrt{C}.
$$
The relative statistical uncertainty is then
$$
\frac{\sigma_C}{C} \approx \frac{1}{\sqrt{C}}.
$$
If each event has a small probability $p$ to contribute to this count, and you simulate $N$ events, then $C \approx pN$, so the relative uncertainty behaves like
$$
\frac{\sigma_C}{C} \approx \frac{1}{\sqrt{pN}} \propto \frac{1}{\sqrt{N}}.
$$
This relationship lets you estimate how many events you need. Suppose you aim for 1 percent relative statistical uncertainty in a simple count. You need
$$
\frac{1}{\sqrt{C}} \approx 0.01 \quad \Rightarrow \quad C \approx 10^4 \text{ counts}.
$$
If your selection is very rare, for example only 1 in $10^5$ events is counted, then to reach $10^4$ counts you need about $10^9$ events. This illustrates why rare processes are expensive to simulate with high precision.
For quantities like energy deposition per event, you often estimate the mean and its standard error. If $\sigma_x$ is the standard deviation of $x$ over events, the standard error of the mean is
$$
\sigma_{\bar{x}} = \frac{\sigma_x}{\sqrt{N}}.
$$
You can estimate $\sigma_x$ from your simulation and then decide how large $N$ must be for $\sigma_{\bar{x}}$ to be smaller than your target uncertainty. For example, if you estimate $\sigma_x \approx 2\ \text{MeV}$ and you want a standard error on the mean dose of $0.02\ \text{MeV}$, then
$$
\sigma_{\bar{x}} = \frac{2\ \text{MeV}}{\sqrt{N}} \le 0.02\ \text{MeV}
\quad \Rightarrow \quad
\sqrt{N} \ge 100 \quad \Rightarrow \quad N \ge 10^4.
$$
In practice you will often work with histograms. Each bin is like a separate count and has its own statistical uncertainty. For a bin with content $C_i$, the standard deviation is $\sigma_i \approx \sqrt{C_i}$. The relative uncertainty in that bin is $1/\sqrt{C_i}$. This means high energy tails or very narrow features in spectra, which receive few counts, have much larger uncertainties than regions around the main peak. If you need good statistical precision in these low count regions, you must strongly increase the total number of events.
It is important to separate statistical uncertainties from systematic ones. Statistical precision improves with more events, while systematic effects, for example an incorrect material definition or a missing physical process, do not improve by running longer. A highly precise but biased result is not valid. You should therefore first validate your model qualitatively, then refine the statistical precision.
A useful strategy is to perform short test runs with a modest number of events to estimate variances and relative uncertainties. From these test runs, you can extrapolate the number of events needed to reach your final precision goal, using the approximate $1/\sqrt{N}$ scaling. You then plan your full production runs accordingly, taking into account available computation time and resources.
Finally, remember that statistical precision is not uniform across all observables. When planning how many events to simulate, always base your decision on the most demanding observable that matters for your study, such as the most finely binned region of interest or the rarest process you need to characterize.
Views: 6