Sampling from energy spectra

For more complex simulations one often needs a beam that is not mono-energetic but follows a specific energy distribution. This feature is provided “out-of-the box” for CERN FLUKA by the
CERN-FLUKA-spectra-sampling.tgz (21.8 KB)
user routine package (version 1.4) that can be fully configured via the input file.
IMPORTANT NOTICE: With the release of FLUKA 4.0.0 the package had to be upgraded. In order to use it with FLUKA version >= 4.0.0 please re-download from the link above.

The provided code also supports the optional propagation of uncertainties of the input spectra to the simulation. In addition, the user can bias the provided distribution to “favor” certain energy regions while not changing the physics. More detailed explanations + an example can be found in the provided archive.

NOTE: Any configuration of the spatial distribution (e.g. volumetric sources) or divergence from the input file are taken into account.

The routines are mainly written in C++ and should also serve as a more complex example of how to combine C/C++ with FLUKA.


As I received a question regarding the units of the input spectra, that could be of interest for more people, I post the corresponding answer below:

The data in the file denotes a distribution which will reflect the probability of a particle having that energy. As a consequence it doesn’t matter which units you use: that could be 1/(cm2 Gev) or 1/cm2 or (1/cm2 * Gev * s) or arbitrary units. The only thing which is important is that relative to each other the values for each bin correctly represent the proportions of the probability of the energy with which you want to see the particles.

For example:

Bin #1 = 0.5
Bin #2 = 1

This means that the chance of having an energy within the limits of bin 1 is only half as big as for bin 2. Thus, 1/3 of the particles will have an energy within bin 1 and 2/3 an energy within bin 2. You could also have used:

Bin #1 = 100
Bin #2 = 200

to obtain the same result.

The shown integral is the integral of the input data over the energy --> Int (f(E) * dE). Naturally all units will be conserved. For example if you’re data were given in 1/(cm2 * GeV) then the integral would be in 1 / cm2 and denote the total fluence. If you’re data were given in 1 / (cm2 * GeV * s) then the integral would be 1 / (cm2 * s) which is the total fluence rate.

As mentioned above for the sampling, which is basically only a game of probability, it doesn’t change anything. But you should keep in mind that all results in FLUKA are given per primary particle. Therefore, if you might still want to normalize your end-results with the total number of particles or total number of particles per unit time, in order to compare them to a real-life experiment.