2.7. Geant4 General Particle Source

2.7.1. Introduction

The G4GeneralParticleSource (GPS) is part of the Geant4 toolkit for Monte-Carlo, high-energy particle transport. Specifically, it allows the specifications of the spectral, spatial and angular distribution of the primary source particles. An overview of the GPS class structure is presented here. Section 2.7.2 covers the configuration of GPS for a user application, and Section 2.7.3 describes the macro command interface. Section 2.7.4 gives an example input file to guide the first time user.

G4GeneralParticleSource is used exactly the same way as G4ParticleGun in a Geant4 application. In existing application one can simply change yourPrimaryGeneratorAction by globally replacing G4ParticleGun with G4GeneralParticleSource. GPS may be configured via command line, or macro based, input. The experienced user may also hard-code distributions using the methods and classes of the GPS that are described in more detail in a technical note [1] .

The class diagram of GPS is shown in Figure 2.1. The G4GeneralParticleSource class can have multiple instanciations of the G4SingleParticleSource class, each with independent positional, angular and energy distributions as well as incident particle types.

The class diagram of G4GeneralParticleSource.

Figure 2.1. The class diagram of G4GeneralParticleSource.


2.7.2. Configuration

GPS allows the user to control the following characteristics of primary particles:

  • Spatial sampling: on simple 2D or 3D surfaces such as discs, spheres, and boxes.
  • Angular distribution: unidirectional, isotropic, cosine-law, beam or arbitrary (user defined).
  • Spectrum: linear, exponential, power-law, Gaussian, blackbody, or piece-wise fits to data.
  • Multiple sources: multiple independent sources can be used in the same run.

As noted above, G4GeneralParticleSource is used exactly the same way as G4ParticleGun in a Geant4 application, and may be substituted for the latter by "global search and replace" in existing application source code.

2.7.2.1. Position Distribution

The position distribution can be defined by using several basic shapes to contain the starting positions of the particles. The easiest source distribution to define is a point source. One could also define planar sources, where the particles emanate from circles, annuli, ellipses, squares or rectangles. There are also methods for defining 1D or 2D accelerator beam spots. The five planes are oriented in the x-y plane. To define a circle one gives the radius, for an annulus one gives the inner and outer radii, and for an ellipse, a square or a rectangle one gives the half-lengths in x and y.

More complicated still, one can define surface or volume sources where the input particles can be confined to either the surface of a three dimensional shape or to within its entire volume. The four 3D shapes used within G4GeneralParticleSource are sphere, ellipsoid, cylinder and parallelepiped.A sphere can be defined simply by specifying the radius. Ellipsoids are defined by giving their half-lengths in x, y and z. Cylinders are defined such that the axis is parallel to the z-axis, the user is therefore required to give the radius and the z half-length. Parallelepipeds are defined by giving x, y and z half-lengths, plus the angles α, θ, and φ (Figure 2.2).

The angles used in the definition of a Parallelepiped.

Figure 2.2. The angles used in the definition of a Parallelepiped.


To allow easy definition of the sources, the planes and shapes are assumed to be orientated in a particular direction to the coordinate axes, as described above. For more general applications, the user may supply two vectors (x' and a vector in the plane x'-y') to rotate the co-ordinate axes of the shape with respect to the overall co-ordinate system (Figure 2.3). The rotation matrix is automatically calculated within G4GeneralParticleSource. The starting points of particles are always distributed homogeneously over the 2D or 3D surfaces, although biasing can change this.

An illustration of the use of rotation matrices. A cylinder is defined with its axis parallel to the z-axis (black lines), but the definition of 2 vectors can rotate it into the frame given by x', y', z' (red lines).

Figure 2.3. An illustration of the use of rotation matrices. A cylinder is defined with its axis parallel to the z-axis (black lines), but the definition of 2 vectors can rotate it into the frame given by x', y', z' (red lines).


2.7.2.2. Angular Distribution

The angular distribution is used to control the directions in which the particles emanate from/incident upon the source point. In general there are three main choices, isotropic, cosine-law or user-defined. In addition there are options for specifying parallel beam as well as diversed accelerator beams. The isotropic distribution represents what would be seen from a uniform 4π flux. The cosine-law represents the distribution seen at a plane from a uniform 2π flux.

It is possible to bias (Section 2.7.2.4) both θ and φ for any of the predefined distributions, including settin lower and upper limits to θ and φ. User-defined distributions cannot be additionally biased (any bias should obviously be incorporated into the user definition).

Incident with zenith angle θ=0 means the particle is travelling along the -z axis. It is important to bear this in mind when specifying user-defined co-ordinates for angular distributions. The user must be careful to rotate the co-ordinate axes of the angular distribution if they have rotated the position distribution (Figure 2.3).

The user defined distribution requires the user to enter a histogram in either θ or φ or both. The user-defined distribution may be specified either with respect to the coordinate axes or with respect to the surface-normal of a shape or volume. For the surface-normal distribution, θ should only be defined between 0 and π/2, not the usual 0 to π range.

The top-level /gps/direction command uses direction cosines to specify the primary particle direction, as follows:

Px = - sin θ cos φ
Py = - sin θ sin φ
Pz = - cos θ

Equation 2.1. Direction cosine calculations


2.7.2.3. Energy Distribution

The energy of the input particles can be set to follow several built-in functions or a user-defined one, as shown in Table 2.1. The user can bias any of the pre-defined energy distributions in order to speed up the simulation (user-defined distributions are already biased, by construction).

SpectrumAbbreviationFunctional FormUser Parameters
mono-energeticMonoI ∝ δ(E-E0)Energy E0
linearLinI ∝ I0 + m × EIntercept I0, slope m
exponentialExpI ∝ exp(-E/E0)Energy scale-height E0
power-lawPowI ∝ EαSpectral index α
GaussianGaussI = (2πσ) exp[-(E-E0)² / σ²]Mean energy E0, standard deviation σ
bremsstrahlungBremI = ∫ 2E² [ h²c² (exp(-E/kT) - 1)]-1Temperature T
black bodyBbodyI ∝ (kT) E exp(-E/kT)Temperature T (see note below)
cosmic diffuse gamma rayCdgI ∝ [ (E/Eb)α1 + (E/Eb)α2 ]-1Energy range Emin to Emax; break energy Eb and indices α1 and α2 are fixed (see note below)

Table 2.1. 


There is also the option for the user to define a histogram in energy ("User") or energy per nucleon ("Epn") or to give an arbitrary point-wise spectrum ("Arb")that can be fit with various simple functions. The data for histograms or point spectra must be provided in ascending bin (abscissa) order. The point-wise spectrum may be differential (as with a binned histogram) or integral (a cumulative distribution function). If integral, the data must satsify s(e1) ≥ s(e2) for e1 < e2 when entered; this is not validated by the GPS code. The maximum energy of an integral spectrum is defined by the last-but-one data point, because GPS converts to a differential spectrum internally.

Unlike the other spectral distributions it has proved difficult to integrate indefinitely the black-body spectrum and this has lead to an alternative approach. Instead it has been decided to use the black-body formula to create a 10,000 bin histogram and then to produce random energies from this.

Similarly, the broken power-law for cosmic diffuse gamma rays makes generating an indefinite integral CDF problematic. Instead, the minimum and maximum energies specified by the user are used to construct a definite-integral CDF from which random energies are selected.

2.7.2.4. Biasing

The user can bias distributions by entering a histogram. It is the random numbers from which the quantities are picked that are biased and so one only needs a histogram from 0 to 1. Great care must be taken when using this option, as the way a quantity is calculated will affect how the biasing works, as discussed below. Bias histograms are entered in the same way as other user-defined histograms.

When creating biasing histograms it is important to bear in mind the way quantities are generated from those numbers. For example let us compare the biasing of a θ distribution with that of a φ distribution. Let us divide the θ and φ ranges up into 10 bins, and then decide we want to restrict the generated values to the first and last bins. This gives a new φ range of 0 to 0.628 and 5.655 to 6.283. Since φ is calculated using φ = 2π × RNDM, this simple biasing will work correctly.

If we now look at θ, we expect to select values in the two ranges 0 to 0.314 (for 0 ≤ RNDM ≤ 0.1) and 2.827 to 3.142 (for 0 ≤ RNDM ≤ 0.9). However, the polar angle θ is calculated from the formula θ = cos-1(1 - 2×RNDM). From this, we see that 0.1 gives a θ of 0.644 and a RNDM of 0.9 gives a θ of 2.498. This means that the above will not bias the distribution as the user had wished. The user must therefore take into account the method used to generate random quantities when trying to apply a biasing scheme to them. Some quantities such as x, y, z and φ will be relatively easy to bias, but others may require more thought.

2.7.2.5. User-Defined Histograms

The user can define histograms for several reasons: angular distributions in either θ or φ; energy distributions; energy per nucleon distributions; or biasing of x, y, z, θ, φ, or energy. Even though the reasons may be different the approach is the same.

To choose a histogram the command /gps/hist/type is used (Section 2.7.3). If one wanted to enter an angular distribution one would type "theta" or "phi" as the argument. The histogram is loaded, one bin at a time, by using the /gps/hist/point command, followed by its two arguments the upper boundary of the bin and the weight (or area) of the bin. Histograms are therefore differential functions.

Currently histograms are limited to 1024 bins. The first value of each user input data pair is treated as the upper edge of the histogram bin and the second value is the bin content. The exception is the very first data pair the user input whose first value is the treated as the lower edge of the first bin of the histogram, and the second value is not used. This rule applies to all distribution histograms, as well as histograms for biasing.

The user has to be aware of the limitations of histograms. For example, in general θ is defined between 0 and π and φ is defined between 0 and 2π, so histograms defined outside of these limits may not give the user what they want (see also Section 2.7.2.4).

2.7.3. Macro Commands

G4GeneralParticleSource can be configured by typing commands from the /gps command directory tree, or including the /gps commands in a g4macro file.

2.7.3.1. G4ParticleGun equivalent commands

CommandArgumentsDescription and restrictions
/gps/particlenameDefines the particle type [default geantino], using Geant4 naming convention.
/gps/ionZ A Q EAfter /gps/particle ion, sets the properties (atomic number Z, atomic mass A, ionic charge Q, excitation energy E in keV) of the ion.
/gps/List List available incident particles
/gps/energyE unitSets the energy [default 1 MeV] for mono-energetic sources. The units can be eV, keV, MeV, GeV, TeV or PeV.
/gps/positionX Y Z unitSets the centre co-ordinates (X,Y,Z) of the source [default (0,0,0) cm]. The units can be micron, mm, cm, m or km.
/gps/directionPx Py PzSet the momentum direction [default (1,0,0)] of generated particles using direction cosines (Equation 2.1).
/gps/timet0 unitSets the primary particle (event) time [default 0 ns]. The units can be ps, ns, us, ms, or s.
/gps/polarizationPx Py PzSets the polarization vector of the source, which does not need to be a unit vector.
/gps/numberNSets the number of particles [default 1] to simulate on each event.
/gps/verboselevelControl the amount of information printed out by the GPS code. Larger values produce more detailed output.

Table 2.2. 


2.7.3.2. Multiple source specification

CommandArgumentsDescription and restrictions
/gps/source/list List the particle sources defined.
/gps/source/addintensityAdd a new particle source with the specified intensity
/gps/source/deleteindexRemove the specified particle source.
/gps/source/clear remove all defined particle sources.
/gps/source/show display the current particle source
/gps/source/setindexSelect the specified particle source as the current one.
/gps/source/multiplevertexflagSpecify true for simulaneous generation of mutiple vertices, one from each specified source. False [default] generates a single vertex, choosing one source randomly.
/gps/source/intensityintensityReset the current source to the specified intensity

Table 2.3. 


2.7.3.3. Source position and structure

CommandArgumentsDescription and restrictions
/gps/pos/typedistSets the source positional distribution type: Point [default], Plane, Beam, Surface, Volume.
/gps/pos/shapeshapeSets the source shape type, after /gps/pos/type has been used. For a Plane this can be Circle, Annulus, Ellipse, Square, Rectangle. For both Surface or Volume sources this can be Sphere, Ellipsoid, Cylinder, Para (parallelpiped).
/gps/pos/centreX Y Z unitSets the centre co-ordinates (X,Y,Z) of the source [default (0,0,0) cm]. The units can only be micron, mm, cm, m or km.
/gps/pos/rot1R1x R1y R1zDefines the first (x' direction) vector R1 [default (1,0,0)], which does not need to be a unit vector, and is used together with /gps/pos/rot2 to create the rotation matrix of the shape defined with /gps/shape.
/gps/pos/rot2R2x R2y R2zDefines the second vector R2 in the xy plane [default (0,1,0)], which does not need to be a unit vector, and is used tohgether with /gps/pos/rot1 to create the rotation matrix of the shape defined with /gps/shape.
/gps/pos/halfxlen unitSets the half-length in x [default 0 cm] of the source. The units can only be micron, mm, cm, m or km.
/gps/pos/halfylen unitSets the half-length in y [default 0 cm] of the source. The units can only be micron, mm, cm, m or km.
/gps/pos/halfzlen unitSets the half-length in z [default 0 cm] of the source. The units can only be micron, mm, cm, m or km.
/gps/pos/radiuslen unitSets the radius [default 0 cm] of the source or the outer radius for annuli. The units can only be micron, mm, cm, m or km.
/gps/pos/inner_radiuslen unitSets the inner radius [default 0 cm] for annuli. The units can only be micron, mm, cm, m or km.
/gps/pos/sigmarsigma unitSets the transverse (radial) standard deviation [default 0 cm] of beam position profile. The units can only be micron, mm, cm, m or km.
/gps/pos/sigmaxsigma unitSets the standard deviation [default 0 cm] of beam position profile in x-direction. The units can only be micron, mm, cm, m or km.
/gps/pos/sigmaysigma unitSets the standard deviation [default 0 cm] of beam position profile in y-direction. The units can only be micron, mm, cm, m or km.
/gps/pos/paralpalpha unitUsed with a Parallelepiped. The angle [default 0 rad] α formed by the y-axis and the plane joining the centre of the faces parallel to the zx plane at y and +y. The units can only be deg or rad.
/gps/pos/parthetheta unitUsed with a Parallelepiped. Polar angle [default 0 rad] θ of the line connecting the centre of the face at z to the centre of the face at +z. The units can only be deg or rad.
/gps/pos/parphiphi unitUsed with a Parallelepiped. The azimuth angle [default 0 rad] φ of the line connecting the centre of the face at z with the centre of the face at +z. The units can only be deg or rad.
/gps/pos/confinenameAllows the user to confine the source to the physical volume name [default NULL].

Table 2.4. 


2.7.3.4. Source direction and angular distribution

CommandArgumentsDescription and restrictions
/gps/ang/typeAngDisSets the angular distribution type (iso [default], cos, planar, beam1d, beam2d, focused, user) to either isotropic, cosine-law or user-defined.
/gps/ang/rot1AR1x AR1y AR1zDefines the first (x' direction) rotation vector AR1 [default (1,0,0)] for the angular distribution and is not necessarily a unit vector. Used with /gps/ang/rot2 to compute the angular distribution rotation matrix.
/gps/ang/rot2AR2x AR2y AR2zDefines the second rotation vector AR2 in the xy plane [default (0,1,0)] for the angular distribution, which does not necessarily have to be a unit vector. Used with /gps/ang/rot2 to compute the angular distribution rotation matrix.
/gps/ang/minthetaMinTheta unitSets a minimum value [default 0 rad] for the θ distribution. Units can be deg or rad.
/gps/ang/maxthetaMaxTheta unitSets a maximum value [default π rad] for the θ distribution. Units can be deg or rad.
/gps/ang/minphiMinPhi unitSets a minimum value [default 0 rad] for the φ distribution. Units can be deg or rad.
/gps/ang/maxphiMaxPhi unitSets a maximum value [default 2π rad] for the φ distribution. Units can be deg or rad.
/gps/ang/sigmarsigma unitSets the standard deviation [default 0 rad] of beam directional profile in radial. The units can only be deg or rad.
/gps/ang/sigmaxsigma unitSets the standard deviation [default 0 rad] of beam directional profile in x-direction. The units can only be deg or rad.
/gps/ang/sigmaysigma unitSets the standard deviation [default 0 rad] of beam directional profile in y-direction. The units can only be deg or rad.
/gps/ang/focuspointX Y Z unitSet the focusing point (X,Y,Z) for the beam [default (0,0,0) cm]. The units can only be micron, mm, cm, m or km.
/gps/ang/user_coorboolCalculate the angular distribution with respect to the user definded co-ordinate system (true), or with respect to the global co-ordinate system (false, default).
/gps/ang/surfnormboolAllows user to choose whether angular distributions are with respect to the co-ordinate system (false, default) or surface normals (true) for user-defined distributions.

Table 2.5. 


2.7.3.5. Energy spectra

CommandArgumentsDescription and restrictions
/gps/ene/typeEnergyDisSets the energy distribution type to one of (Table 2.1):
Mono (mono-energetic, default)
Lin (linear)
Pow (power-law)
Exp (exponential)
Gauss (Gaussian)
Brem (bremsstrahlung)
Bbody (black-body)
Cdg (cosmic diffuse gamma-ray)
User (user-defined histogram)
Arb (point-wise spectrum)
Epn (energy-per-nucleon histogram)
/gps/ene/minEmin unitSets the minimum [default 0 keV] for the energy distribution. The units can be eV, keV, MeV, GeV, TeV or PeV.
/gps/ene/maxEmax unitSets the maximum [default 0 keV] for the energy distribution. The units can be eV, keV, MeV, GeV, TeV or PeV.
/gps/ene/monE unitSets the energy [default 1 MeV] for mono-energetic sources. The units can be eV, keV, MeV, GeV, TeV or PeV.
/gps/ene/sigmasigma unitSets the standard deviation [default 0 keV] in energy for Gaussian or Mono energy distributions. The units can be eV, keV, MeV, GeV, TeV or PeV.
/gps/ene/alphaalphaSets the exponent α [default 0] for a power-law distribution.
/gps/ene/tempTSets the temperature in kelvins [default 0] for black body and bremsstrahlung spectra.
/gps/ene/ezeroE0Sets scale E0 [default 0] for exponential distributions.
/gps/ene/gradientgradientSets the gradient (slope) [default 0] for linear distributions.
/gps/ene/interceptinterceptSets the Y-intercept [default 0] for the linear distributions.
/gps/ene/calculate Prepares integral PDFs for the interally-binned cosmic diffuse gamma ray (Cdg) and black body (Bbody) distributions.
/gps/ene/emspecboolAllows user to specify distributions are in momentum (false) or energy (true, default). Only valid for User and Arb distributions.
/gps/ene/diffspecboolAllows user to specify whether a point-wise spectrum is integral (false) or differential (true, default). The integral spectrum is only usable for Arb distributions.

Table 2.6. 


2.7.3.6. User-defined histograms and interpolated functions

CommandArgumentsDescription and restrictions
/gps/hist/typetypeSet the histogram type: predefined biasx [default], biasy, biasz, biast (angle θ), biasp (angle φ), biaspt (position θ), biaspp (position φ), biase; user-defined histograms theta, phi, energy, arb (point-wise), epn (energy per nucleon).
/gps/hist/resettypeRe-set the specified histogram: biasx [default], , biasy, biasz, biast, biasp, biaspt, biaspp, biase, theta, phi, energy, arb, epn.
/gps/hist/pointEhi WeightSpecify one entry (with contents Weight) in a histogram (where Ehi is the bin upper edge) or point-wise distribution (where Ehi is the abscissa). The abscissa Ehi must be in Geant4 default units (MeV for energy, rad for angle).
/gps/hist/intertypeSets the interpolation type (Lin linear, Log logarithmic, Exp exponential, Spline cubic spline) for point-wise spectra. This command must be issued immediately after the last data point.

Table 2.7. 


2.7.4. Example Macro File

	# Macro test2.g4mac
	/control/verbose 0
	/tracking/verbose 0
	/event/verbose 0
	/gps/verbose 2
	/gps/particle gamma
	/gps/pos/type Plane
	/gps/pos/shape Square
	/gps/pos/centre 1 2 1 cm
	/gps/pos/halfx 2 cm
	/gps/pos/halfy 2 cm
	/gps/ang/type cos
	/gps/ene/type Lin
	/gps/ene/min 2 MeV
	/gps/ene/max 10 MeV
	/gps/ene/gradient 1
	/gps/ene/intercept 1
	/run/beamOn 10000

The above macro defines a planar source, square in shape, 4 cm by 4 cm and centred at (1,2,1) cm. By default the normal of this plane is the z-axis. The angular distribution is to follow the cosine-law. The energy spectrum is linear, with gradient and intercept equal to 1, and extends from 2 to 10 MeV. 10,000 primaries are to be generated.

Figure 4. Energy, position and angular distributions of the primary particles as generated by the macro file shown above.

Figure 2.4. Figure 4. Energy, position and angular distributions of the primary particles as generated by the macro file shown above.


The standard Geant4 output should show that the primary particles start from between 1, 0, 1 and 3, 4, 1 (in cm) and have energies between 2 and 10 MeV, as shown in Figure 2.4, in which we plotted the actual energy, position and angular distributions of the primary particles generated by the above macro file.



[1] General purpose Source Particle Module for Geant4/SPARSET: Technical Note, UoS-GSPM-Tech, Issue 1.1, C Ferguson, February 2000.