1. User guide¶

1.1. Introduction & installation¶

1.1.1. Why this new project ?¶

The GATE project is more than 15 years old. During this time, it evolves a lot, it now allows to perform a wide range of medical physics simulations such as various imaging systems (PET, SPECT, Compton Cameras, X-ray, etc) and dosimetry studies (external and internal radiotherapy, hadrontherapy, etc). This project led to hundreds of scientific publications, contributing to help researchers and industry.

GATE fully relies on Geant4 for the Monte Carlo engine and provides 1) easy access to Geant4 functionalities, 2) additional features (e.g. variance reduction techniques) and 3) collaborative development to shared source code, avoiding reinventing the wheel. The user interface is done via so-called macro files (.mac) that contain Geant4 style macro commands that are convenient compared to direct Geant4 C++ coding. Note that other projects such as Gamos or Topas also rely on similar principles.

Since the beginning of GATE, a lot of changes have happened in both fields of computer science and medical physics, with, among others, the rise of machine learning and Python language, in particular for data analysis. Also, the Geant4 project is still very active and is guaranteed to be maintained at least for the ten next years (as of 2020).

Despite its usefulness and its unique features (collaborative, open source, dedicated to medical physics), we think that the GATE software in itself, from a computer science programming point of view, is showing its age. More precisely, the source code has been developed during 15 years by literally hundreds of different developers. The current GitHub repository indicates around 70 unique contributors, but it has been set up only around 2012 and a lot of early contributors are not mentioned in this list. This diversity is the source of a lot of innovation and experiments (and fun!), but also leads to maintenance issues. Some parts of the code are “abandoned”, some others are somehow duplicated. Also, the C++ language evolves tremendously during the last 15 years, with very efficient and convenient concepts such as smart pointers, lambda functions, ‘auto’ keyword … that make it more robust and easier to write and maintain.

Keeping in mind the core pillars of the initial principles (community-based, open-source, medical physics oriented), we decide to start a project to propose a new way to perform Monte Carlo simulations in medical physics. Please, remember this is an experimental (crazy ?) attempt, and we are well aware of the very long and large effort it requires to complete it. At time of writing, it is not known if it can be achieved, so we encourage users to continue using the current GATE version for their work. Audacious users may nevertheless try this new system and make feedback. Mad ones can even contribute …

Never stop exploring !

(2020)

1.1.2. Goals and features¶

The purpose of this software is to facilitate the creation of Geant4-based Monte Carlo simulations for medical physics using Python as the primary scripting language. The user interface has been redesigned to allow for direct creation of simulations in Python, rather than using macro files.

Some key features of this software include:

Use of Python as the primary scripting language for creating simulations
Multithreading support for efficient simulation execution
Native integration with ITK for image management
Compatibility with Linux, Mac, and potentially Windows operating systems
Convenient installation via a single pip install opengate command
…

1.1.3. Installation (for users, not for developers)¶

You only have to install the Python module with, the –pre option is mandatory to have the latest release:

pip install --pre opengate

Then, you can create a simulation using the opengate module (see below). For developers, please look the developer guide for the developer installation.

{tip} We highly recommend creating a specific python environment to 1) be sure all dependencies are handled properly and 2) don’t mix with your other Python modules. For example, you can use venv. Once the environment is created, you need to activate it:

python -m venv opengate_env
source opengate_env/bin/activate
pip install --pre opengate

or with conda environment:

conda create --name opengate_env python=3.10
conda activate opengate_env
pip install --pre opengate

Maybe you need to upgrade the pip module with:

pip install --upgrade pip

If you already installed opengate, just upgrade it with:

pip install --upgrade --pre opengate

Once installed, we recommend to check the installation by printing GATE information and running the tests:

opengate_info
opengate_tests

WARNING 1 The first time a simulation is executed, the Geant4 data must be downloaded and installed. This step is automated but can take some times according to your bandwidth. Note that this is only done once. Running opengate_info will print some details and the path of the data.

WARNING 2 With some linux systems (not all), you may encounter an error similar to “cannot allocate memory in static TLS block”. In that case, you must add a specific path to the linker as follows:

export LD_PRELOAD=<path to libG4processes>:<path to libG4geometry>:${LD_PRELOAD}

The libraries (libG4processes and libG4geometry) are usually found in the Geant4 folder, something like ~/build-geant4.11.0.2/BuildProducts/lib64.

1.1.4. Additional command lines tools¶

There is some additional commands lines tools that can also be used, see the addons section.

1.1.5. Teaching resources and examples¶

Warning they are only updated infrequently, you may have to adapt them to changes in the gate version.

exercices (initially developed for DQPRM, French medical physics diploma)
exercices (initially developed for Opengate teaching)

1.2. Simulation¶

1.2.1. Units values¶

The Geant4 physics units can be retrieved with the following:

import opengate as gate

cm = gate.g4_units.cm
eV = gate.g4_units.eV
MeV = gate.g4_units.MeV
x = 32 * cm
energy = 150 * MeV
print(f'The energy is {energy/eV} eV')

The units behave like in Geant4 system of units.

1.2.2. Main Simulation object¶

Any simulation starts by defining the (unique) Simulation object. The generic options can be set with the user_info data structure (a kind of dictionary), as follows. You can print this user_info data structure to see all available options with the default value print(sim.user_info).

sim = gate.Simulation()
ui = sim.user_info
print(ui)
ui.verbose_level = gate.LOG_DEBUG
ui.running_verbose_level = gate.LOG_EVENT
ui.g4_verbose = False
ui.g4_verbose_level = 1
ui.visu = False
ui.visu_verbose = False
ui.random_engine = 'MersenneTwister'
ui.random_seed = 'auto'
ui.number_of_threads = 1

A simulation must contains 4 main elements that define a complete simulation:

Geometry: all geometrical elements that compose the scene, such as phantoms, detectors, etc.
Sources: all sources of particles that will be created ex-nihilo. Each source may have different properties (location, direction, type of particles with their associated energy ,etc).
Physics: describe the properties of the physics models that will be simulated. It describes models, databases, cuts etc.
Actors : define what will be stored and output during the simulation. Typically, dose deposition or detected particles. This is the generic term for ‘scorer’. Note that some Actors can not only store and output data, but also interact with the simulation itself (hence the name ‘actor’).

Also, you can use the command line opengate_user_info to print all default and possible parameters.

Each four elements will be described in the following sections. Once they have be defined, the simulation must be initialized and can be started. The initialization corresponds to the Geant4 step needed to create the scene, gather cross-sections, etc.

output = sim.start()

1.2.2.1. Random Number Generator¶

The RNG can be set with ui.random_engine = "MersenneTwister". The default one is “MixMaxRng” and not “MersenneTwister” because it is recommended by Geant4 for MT.

The seed of the RNG can be set with self.random_seed = 123456789, with any number. If you run two times a simulation with the same seed, the results will be exactly the same. There are some exception to that behavior, for example when using PyTorch-based GAN. By default, it is set to “auto”, which means that the seed is randomly chosen.

1.2.2.2. Run and timing¶

The simulation can be split into several runs, each of them with a given time duration. Geometry can only be modified between two runs, not within one. By default, the simulation has only one run with a duration of 1 second. In the following example, we defined 3 runs, the first has a duration of half a second and start at 0, the 2nd run goes from 0.5 to 1 second. The 3rd run starts later at 1.5 second and lasts 1 second.

sim.run_timing_intervals = [
    [ 0, 0.5 * sec],
    [ 0.5 * sec, 1.0 * sec],
    # Watch out : there is (on purpose) a 'hole' in the timeline
    [ 1.5 * sec, 2.5 * sec],
    ]

1.2.2.3. Verbosity (for debug)¶

The verbosity, i.e. the messages printed on the screen, are controlled via various parameters.

ui.verbose_level: can be DEBUG or INFO. Will display more or less messages during initialization
ui.running_verbose_level: can be RUN or EVENT. Will display message during simulation run
ui.g4_verbose: (bool) enable or disable the Geant4 verbose system
ui.g4_verbose_level: level of the Geant4 verbose system
ui.visu_verbose: enable or disable Geant4 verbose during visualisation

1.2.2.4. Visualisation¶

Visualisation is enabled with ui.visu = True. Then, you have the choice to choose between qt, vrml or gdml interface.

1.2.2.4.1. QT¶

It will start a Qt interface with ui.visu_type = "qt". By default, the Geant4 visualisation commands are the ones provided in the file opengate\mac\default_visu_commands.mac. It can be changed with self.visu_commands = gate.read_mac_file_to_commands('my_visu_commands.mac').

The visualisation with qt is still work in progress. First, it does not work on some linux systems (we don’t know why yet). When a CT image is inserted in the simulation, every voxel should be drawn which is highly inefficient and cannot really be used.

1.2.2.4.2. VRML¶

You can choose vrml visualization with ui.visu_type = "vrml". Opengate uses pyvista for the GUI, so you need to install it before with pip install pyvista. Alternatively, if you want to use an external VRML viewer, you can save a VRML file with ui.visu_type = "vrml_file_only". In such case, the GUI is not open, and you do not need pyvista. In both cases, you need to set ui.visu_filename = "geant4VisuFile.wrl" to save the VRML file.

If you want to personalized the pyvista GUI, you can set ui.visu_type = "vrml_file_only" and execute you own code in your python script. You can find an example in test004_simple_visu_vrml.py

1.2.2.4.3. GDML¶

With GDML visualization, you can only view the geometry, not the paths of the particles. It is enabled with ui.visu_type = "gdml". GDML visualization needs to be enabled in Geant4 with GEANT4_USE_GDML=ON during the compilation but you need to have xerces-c available on your computer (install it with yum, brew, or apt-get, …). Opengate uses pyg4ometry for the GUI, so you need to install it with pip install pyg4ometry. pyg4ometry uses opencascade librairy, so install opencascade with your package manager. If you want to use an external GDML viewer, you can save the visualization to a GDML file with ui.visu_type = "gdml_file_only". In such case, the GUI is not open, and you do not need pyg4ometry. In both cases, you need to set ui.visu_filename = "geant4VisuFile.gdml" to save the GDML file.

1.2.2.5. Multithreading¶

Multithreading is enabled with ui.number_of_threads = 4 (larger than 1). When MT is enabled, there will one run for each thread, running in parallel.

Warning, the speedup is not optimal in all situations. First, it takes time to start a new thread, so it the simulation is short, MT does not bring any speedup. Second, if the simulation contains several runs (for moving volumes for example), all runs will be synchronized, i.e. the master thread will wait for all threads to terminate the run before starting another one. This synchronisation takes times and impact the speedup.

However, for other cases, MT is very efficient and brings almost linear speedups, at least for a “low” numbers of threads (we tested it with 8 threads on dose computation, leading to almost x8 time gain).

1.2.3. Starting and SimulationEngine¶

Once all simulation elements have been described (see next sections for more details), the Geant4 engine can be initialized and started. This is done by one single command:

output = sim.start()

Geant4 engine is designed to be the only one instance of the engine, and thus prevent to run two simulations in the same process. In most of the cases, this is not an issue, but sometimes, for example in notebook, we want to run several simulations during the same process session. This can be achieved by setting the following option that will start the Geant4 engine in a separate process and copy back the resulting output in the main process. This is the task of the SimulationEngine object.

se = gate.SimulationEngine(sim, start_new_process=True)
output = se.start()
# or shorter :
output = sim.start(start_new_process=True)

1.2.4. After the simulation¶

Once the simulation is terminated (after the start()), user can retrieve some actor outputs via the output.get_actor function. Note that some output data can be un-available when the simulation is run in a separate process. For the moment, G4 objects (ROOT output) and ITK images cannot be copied back to the main process, e.g. ITK images and ROOT files should be written on disk to be accessed back.

This behavior may change in the future.

1.2.4.1. Multiprocessing (advanced use)¶

The Geant4 simulation engine has a limitation where it can only run one single simulation and cannot be reused in the same process. This can be a problem in certain contexts, such as when using Python Notebooks. To overcome this limitation, multiprocessing can be used to run the simulation in a separate process (not a thread) with its own memory space. The following option can be used to achieve this:

output = sim.start(start_new_process=True)

When this option is used, the Geant4 engine will be created and run in a separate process, which will be terminated after the simulation is finished. The output of the simulation will be copied back to the main process that called the start() function. This allows for the use of Gate in Python Notebooks, as long as this option is not forgotten.

For advanced usage, you can explicitly create the engine for the simulation with:

se = gate.SimulationEngine(sim)
se.start_new_process = True
se.user_fct_after_init = my_function
output = se.start(True)

Here user can also define a function (my_function in the above example) that will be called after the Geant4 engine is initialized, and before it starts the simulation. This function will be called in the newly created process, so all data it accesses must be serializable (Python’s pickable) to be copied to the new process.

1.3. Geometry and volumes¶

Gate fundamentally relies on the geometry principles of Geant4, but provides the user with an easy-to-use interface to set up the geometry of a simulation. In this part of the Gate user guide, we explain how a simulation geometry is set up in Gate.

Under the hood, geometry is handled and parametrized by Geant4. GATE just sets it up for you. Therefore, it might be worthwhile looking at the Geant4 user guide as well.

1.3.1. Overview: Volumes¶

Volumes are the components that make up the simulation geometry. Following Geant4 logic, a volume contains information about its shape, its placement in space, its material, and possibly settings about physics modeling within that volume. In Gate, all these properties are stored and handled in a single volume object, e.g. a BoxVolume, SphereVolume, ImageVolume.

Volumes are managed by the VolumeManager and can be created in two ways:

… with the add_volume command, providing the type of volume as string argument. It is mandatory to provide a unique name as well. A volume is created according to the specified volume type and the volume object is returned. Example:

import opengate as gate
sim = gate.Simulation()
myboxvol = sim.add_volume('Box', name='mybox')

Most users will opt for this way of creating volumes.

… by calling the volume class. In this case, the volume is created, but not yet added to the simulation. It has to be added to the simulation explicitly. Example:

import opengate as gate
sim = gate.Simulation()
myspherevol = gate.geometry.volumes.SphereVolume(name='mysphere')
sim.add_volume(myspherevol)

This second way of creating volumes is useful in cases where the volume is needed but should not be part of the simulation. For example, if it serves as basis for a boolean operation, e.g. to be intersected with another volume.

Note that the add_volume command in the second example does not require the name because the volume already exists and already has a name. For the same reason, the add_volume command does not return anything, i.e. it returns None.

Note: Every simulation has a default volume called world (lowercase) which is automatically created.

The parameters of a volume can be set as follows:

import opengate as gate
sim = gate.Simulation()
vol = sim.add_volume('Box', 'mybox')
vol.material = 'G4_AIR'
vol.mother = 'world'  # by default
cm = gate.g4_units.cm
mm = gate.g4_units.mm
vol.size = [10 * cm, 5 * cm, 15 * mm]

To get an overview of all the properties of a volume, simply print it:

vol = sim.add_volume('Box', 'mybox')
print(vol)

In an interactive python console, e.g. ipython, you can also type help(vol) to get an explanation of the parameters.

To dump a list of all available volume types:

print('Volume types :')
print(sim.volume_manager.dump_volume_types())

1.3.2. Volume hierarchy¶

All volumes have a parameter mother which contains the name of the volume to which they are attached. You can also pass a volume object to the mother parameter and Gate will extract its name from it. By default, a volume’s mother is the world volume (which has the name world). Gate creates a hierarchy of volumes based on each volume’s mother parameter, according to Geant4’s logic of hierarchically nested volumes. The volume hierarchy can be inspected with the command dump_volume_tree of the volume manager. Example:

import opengate as gate
sim = gate.Simulation
b1 = sim.add_volume('Box', name='b1')
b1_a = sim.add_volume('Box', name='b1_a')
b1_b = sim.add_volume('Box', name='b1_b')
b1_a.mother = b1
b1_b.mother = b1
sim.volume_manager.dump_volume_tree()

1.3.3. Common parameters¶

Some of the parameters are common to all volumes, while others are specific to a certain type of volume. Use print(vol) to display the volume’s parameters and their default values.

Common parameters are:

mother: the name of the mother volume (world by default) in the hierarchy of volumes. Volumes are always positioned with respect to the reference frame of the mother volume and therefore moves with the mother volume.
material: the name of the material that composes the volume, e.g. G4_WATER. See section Materials
translation: list of 3 numerical values, e.g. [0, 2*cm, 3*mm]. It defines the translation of the volume with respect to the reference frame of the mother volume. Note: the origin of the reference frame is always at the center of the shape in Geant4.
rotation: a 3x3 rotation matrix. Rotation of the volume with respect to the mother volume. We advocate the use of scipy.spatial.transform.Rotation to manage the rotation matrix.
color: a list of 4 values (Red, Green, Blue, Opacity) between 0 and 1, e.g. [1, 0, 0, 0.5]. Only used when visualization is on.

Take a look at test007 as example for simple volumes.

1.3.4. Utility properties¶

Volume objects come with several properties which allow you to extract information about the volume. The following description assumes that you have created a volume already, i.e.

import opengate as gate
sim = gate.Simulation()
mysphere = sim.add_volume('SphereVolume', name='mysphere')

You can use the following properties to obtain information about the volume mysphere:

mysphere.volume_depth_in_tree: this yields the depth in the hierarchy tree of volumes where 0 is the world, 1 is a volume attached to the world, 2 the first-level subvolume of another volume, and so forth.
mysphere.world_volume: returns the world volume to which this volume is linked through the volume hierarchy. Useful in a simulation with parallel worlds.
mysphere.volume_type: returns the volume type, e.g. “BoxVolume”, “BooleanVolume”, “ImageVolume”. Technically speaking, it yields the name of the volume’s class.
mysphere.bounding_limits: returns the corner coordinates (3 element list: (x,y,z)) of the bounding box of the volume
mysphere.bounding_box_size: returns the size of the bounding box along x, y, z

Note that the above properties are read-only - you cannot set their values.

1.3.5. Materials¶

From the simulation point of view, a material is a set of parameters describing its chemical composition and physical properties such as its density.

Geant4 defines a set of default materials which are also available in GATE. A prominent example is “G4_WATER”. The full list of Geant4 materials is available here.

On top of that, Gate provides different mechanisms to define additional materials. One option is via a text file which can be loaded with

sim.volume_manager.add_material_database("GateMaterials.db")

All material names defined in the “GateMaterials.db” can then be used for any volume. Please check the file in tests/data/GateMaterials.db for the required format of database file.

1.3.6. Image volumes¶

An image volumes is essentially a box filled with a voxelized volumetric (3D) image. The box containing the image behaves pretty much like a BoxVolume and its size is automatically adjusted to match the size of the input image. The image should be provided in a format readable by the itk package and the path to the image file is set via the parameter image. In general, we advocate the use of the mhd/raw file format, but other itk-compatible file formats can be used as well. The image must be 3D, with any pixel type (float, int, char, etc).

From the simulation point of view, a voxel is like a small box through which particles need to be transported. Therefore, in order for Gate/Geant4 to make use of the image, the image values need to be mapped to materials to be associated with the corresponding voxel. To this end, you need to provide a lookup table via the parameter voxel_materials, which is a list of 3-item-lists, each defining a value range and the material name to be used. Take the following example:

import opengate as gate
sim = gate.Simulation()
patient = sim.add_volume("Image", name="patient")
patient.image = "data/myimage.mhd"
patient.mother = "world"
patient.material = "G4_AIR"  # material used by default
patient.voxel_materials = [
  [-2000, -900, "G4_AIR"],
  [-900, -100, "Lung"],
  [-100, 0, "G4_ADIPOSE_TISSUE_ICRP"],
  [0, 300, "G4_TISSUE_SOFT_ICRP"],
  [300, 800, "G4_B-100_BONE"],
  [800, 6000, "G4_BONE_COMPACT_ICRU"],
]
patient.dump_label_image = "labels.mhd"

In the example above, the material “Lung” will be assigned to every voxel with a value between -900 and -100. Voxels whose value does not fall into any of the intervals are considered to contain the volume’s default material, i.e. patient.material = "G4_AIR" in the example above. If a path is provided as dump_label_image parameter of the image volume, an image will be written to the provided path containing material labels. Label 0 stands for voxels to which the default material was assigned, and labels greater than 1 represent all other materials, in ascending order of the lower interval bounds provided in voxel_materials. In the example above, voxels with label 3 correspond to “G4_ADIPOSE_TISSUE_ICRP”, voxels with label 4 correspond to “G4_TISSUE_SOFT_ICRP”, and so forth. See test test009 as an example simulation using an Image volume.

The frame of reference of an Image is linked to the bounding box and treated like other Geant4 volumes, i.e. by default, the center of the image box is positioned at the origin of the mother volume’s frame of reference. Important: Currently, the origin provided by the input image (e.g. in the DICOM or mhd file) is ignored. If you want to place the Image volume according to the origin and rotation provided by the input image, you need to extract that information and set it via the translation and rotation parameters of the image volume. A future version of Gate 10 might provide an option to do this automatically. If you are motivated, you can implement that feature and contribute it to the opengate package.

There is a helper function HounsfieldUnit_to_material to create an interval-material list that can be used as input to the voxel_materials parameter, specifically for CT images expressed in Hounsfield Units:

import opengate as gate
sim = gate.Simulation()
gcm3 = gate.g4_units.g_cm3
f1 = "PATH_TO_OPENGATE/tests/data/Schneider2000MaterialsTable.txt"
f2 = "PATH_TO_OPENGATE/tests/data/Schneider2000DensitiesTable.txt"
tol = 0.05 * gcm3
voxel_materials, materials = gate.geometry.materials.HounsfieldUnit_to_material(sim, tol, f1, f2)

The function HounsfieldUnit_to_material returns two objects:

A list of intervals and material names which can be used as parameter voxel_materials
A list of materials for other use

The input parameters of the function HounsfieldUnit_to_material are

An existing simulation (here sim)
The density tolerance (in g/cm3)
The path to a file containing a list of reference materials
The path to a file containing a list of reference densities

Examples of such files can be found in the opengate/tests/data folder. See test test009 as example.

1.3.7. Tesselated (STL) volumes¶

It is possible to create a tesselated volume shape based on an Standard Triangle Language (STL) data file. Such a file contains a mesh of triangles for one object. It is a typical output format of Computer Aided Design (CAD) software. To create such a volume add a volume of type “Tesselated”. Please keep in mind, that no material information is provided, it has to be specified by the user. A Tesselated volume inherits the the same basic options as other solids described above such as translation or rotation. A basic example how to import an STL file into a geometry “MyTesselatedVolume” and assign the material G4_WATER to it can be found below. In order to verify the correct generation of the solid, one could look at the volume.

import opengate as gate
sim = gate.Simulation()
tes = sim.add_volume("Tesselated", name="MyTesselatedVolume")
tes.material = "G4_WATER"
tes.mother = "world"  # by default
tes.file_name = "myTesselatedVolume.stl"
#to read the volume of the generated solid
print("volume: ",sim.volume_manager.get_volume(
        "MyTesselatedVolume"
    ).solid_info.cubic_volume)
#an alternative way read the volume of the generated solid
print("same volume: ",tes.solid_info.cubic_volume)

See test test067_stl_volume for example.

1.3.8. Repeated volumes¶

The first method, described in this section, is controlled via the translation and rotation parameters. To instruct Geant4 to repeat a volume in multiple locations, it is sufficient to provide a list of translation vectors to the volume parameter translation. Gate will make sure that a G4PhysicalVolume is created for each entry. Consequently, the length of the list of translations determines the number of copies. If only a single rotation matrix is provided as volume parameter rotation, this will be used for all copies. If each copies requires a separate individual rotation, e.g. when repeating volume around a circle, then the volume parameter rotation should receive a list of rotation matrices. Obviously, the number of rotations and translation should match.

Each volume copy corresponds to a G4PhysicalVolume in Geant4 with its own unique name. Gate automatically generates this name. It can be obtained from a given copy index (counting starts at 0) via the method get_repetition_name_from_index(). Or vice versa, the copy index can be obtained from the copy name via get_repetition_index_from_name().

Gate comes with utility functions to generate translation and rotation parameters for common types of volume repetitions - see below.

import opengate as gate
from scipy.spatial.transform import Rotation

cm = gate.g4_units.cm
crystal = sim.add_volume("Box", "crystal")
crystal.size = [1 * cm, 1 * cm, 1 * cm]
crystal.material = "LYSO"
m = Rotation.identity().as_matrix()
crystal.translation = [[1 * cm, 0 * cm, 0],
                       [0.2 * cm, 2 * cm, 0],
                       [-0.2 * cm, 4 * cm, 0]]
print(f"The crystal is repeated in {crystal.number_of_repetitions} locations. ")
print(f"Specified by the following translation vectors: ")
for i, t in enumrate(crystal.translation):
    print(f"Repetition {crystal.get_repetition_name_from_index(i)}: {t}. ")

In this example, the volume named crystal, with the shape of a box, a size of 1x1x1 cm³, and made of LYSO, is repeated in 3 positions. In this example, only the translation is modified, the rotation is set to the same default identity matrix.

There are utility functions that help you to generate lists of translations and rotations. For example:

import opengate as gate
mm = gate.g4_units.mm
crystal = sim.add_volume("Box", "crystal")
translations_grid = gate.geometry.utility.get_grid_repetition(size=[1, 4, 5], spacing=[0, 32.85 * mm, 32.85 * mm])
crystal.translation = translations_grid
# or
detector = sim.add_volume("Box", "detector")
translations_circle, rotations_circle = gate.geometry.utility.get_circular_repetition(number_of_repetitions=18, first_translation=[391.5 * mm, 0, 0], axis=[0, 0, 1])
detector.translation = translations_circle
detector.rotation = rotations_circle

To get help about the utility functions, do:

import opengate as gate
help(gate.geometry.utility.get_grid_repetition)
help(gate.geometry.utility.get_circular_repetition)

You can also have a look at the philipsvereos.py and siemensbiograph.py examples in the opengate/contrib/pet/ folder.

You are obviously free to generate your own list of translations and rotations to suit your needs and they do not need to be regularly spaced and/or follow any spatial pattern such as a grid or ring. Just remember that Geant4 does not allow volumes to overlap and make sure that repetitions to not geometrically interfere (overlap) with each other.

Volume repetitions controlled via the translation and rotation parameter are a convenient and generic way to construct a “not too large” number of repeated objects. In case of “many” repetitions, the Geant4 tracking engine can become slow. In that case, it is better to use parameterised volumes described in the next section. It is not easy to quantify “not too many” repetitions. Based on our experience, a few hundred is still acceptable, but you might want to check in your case. Note that, if the volume contains sub-volumes (via their mother parameter, everything will be repeated, albeit in an optimized and efficient way.

1.3.9. Repeat Parametrised Volumes¶

In some situations, the repeater concept explained in the previous section is not sufficient and can be inefficient when the number of repetitions is large. A specific example is a collimator for SPECT imaging containing a large number of holes. RepeatParametrisedVolume is an alternative repeated volume type which suits this use case. See this example:

import opengate as gate
mm = gate.g4_units.mm
crystal = sim.add_volume("Box", "crystal")
param_vol = sim.add_volume("RepeatParametrised", f"my_param")
param_vol.repeated_volume_name = "crystal"
param_vol.translation = None
param_vol.rotation = None
size = [183, 235, 1]
tr = [2.94449 * mm, 1.7 * mm, 0]
param_vol.linear_repeat = size
param_vol.translation = tr
param_vol.start = [-(x - 1) * y / 2.0 for x, y in zip(size, tr)]
param_vol.offset_nb = 1
param_vol.offset = [0, 0, 0]

Note that the RepeatParametrisedVolume is still partly work in progress. The user guide on this will soon be updated and extended.

param = sim.add_volume("RepeatParametrised", f"my_param")
param.repeated_volume_name = "crystal"
param.translation = None
param.rotation = None
size = [183, 235, 1]
tr = [2.94449 * mm, 1.7 * mm, 0]
param.linear_repeat = size
param.translation = tr
param.start = [-(x - 1) * y / 2.0 for x, y in zip(size, tr)]
param.offset_nb = 1
param.offset = [0, 0, 0]

1.3.10. Boolean volumes¶

Geant4 provides a mechanism to combine volumetric shapes (Solids in Geant4) into new ones via boolean operations, i.e. union, intersection, and subtraction. In GATE, the details of this mechanism are taken care of under the hood and the user can directly combine compatible volumes. For example:

import opengate as gate
from scipy.spatial.transform import Rotation

sim = gate.Simulation()
cm = gate.g4_units.cm
b = gate.geometry.volumes.BoxVolume(name="box")
b.size = [10 * cm, 10 * cm, 10 * cm]
s = gate.geometry.volumes.SphereVolume(name="sph")
s.rmax = 5 * cm
t = gate.geometry.volumes.TubsVolume(name="t")
t.rmin = 0
t.rmax = 2 * cm
t.dz = 15 * cm

combined_b_s = gate.geometry.volumes.unite_volumes(b, s, translation=[0, 1 * cm, 5 * cm])
final_vol = gate.geometry.volumes.subtract(combined_b_s, t, rotation=Rotation.from_euler("x", 3, degrees=True).as_matrix())

final_vol.translation = [5 * cm, 5 * cm, 5 * cm]
final_vol.mother = "world"
final_vol.material = "G4_WATER"
sim.add_volume(final_vol)

The keyword arguments translation and rotation specify how the second shape is translated and rotated, respectively, with respect to the first shape prior to the boolean operation. The absolute placement in space in the simulation is irrelevant for this. On the other hand, the line final_vol.translation = [5 * cm, 5 * cm, 5 * cm] simply refers to the [common parameter](#Common parameters) which specifies the placement of the final volume in space with respect to its mother, in this case the world volume.

Only the finally resulting volume final_vol is actually added to the simulation while the others are only created as intermediate steps of the contruction.

Note that not all volumes are compatible with boolean operations. For example, image volumes cannot be combined. You will receive an error message when trying to apply booelan operations to incompatible volumes.

Boolean operations are a great tool to build complex shapes. The phantoms in opengate.contrib.phantoms are good examples. Also have a look at test016. Be aware, however, that the Geant4 user guide warns that very extensive use of boolean operations can slow down particle tracking speed.

1.3.11. Examples of complex geometries: Linac, SPECT, PET, phantoms¶

Examples of complex nested geometries, partly relying on boolean and repeat operations, can be found in the subpackages opengate.contrib.pet, opengate.contrib.spect, opengate.contrib.linacs, opengate.contrib.phantoms. Also have a look at some of the tests that use these geometries, e.g. test015 (iec phantom), test019 (linac Elekta), test028 (SPECT GE NM670), test037 (Philips Vereos PET).

1.3.12. Parallel worlds¶

TODO

1.4. Sources¶

Sources are the objects that create particles ex nihilo. The particles created from sources are called the Event in the Geant4 terminology, they got a EventID which is unique in a given Run.

Several sources can be defined and are managed at the same time. To add a source description to the simulation, you do:

source1 = sim.add_source('Generic', 'MySource')
source1.n = 100

Bq = gate.g4_units('Bq')
source2 = sim.add_source('Voxels', 'MySecondSource')
source2.activity = 10 * Bq

There are several source types, each one with different parameter. In this example, source1.n indicates that this source will generate 10 Events. The second source manages the time and will generate 10 Events per second, so according to the simulation run timing, a different number of Events will be generated.

Information about the sources may be displayed with:

# Print all types of source
print(sim.dump_source_types())

# Print information about all sources
print(sim.dump_sources())

1.4.1. Generic sources¶

The main type of source is called ‘GenericSource’ that can be used to describe a large range of simple source types. With ‘GenericSource’, user must describe 1) particle type, 2) position, 3) direction and 4) energy, see the following example:

from scipy.spatial.transform import Rotation  # used to describe a rotation matrix

MeV = gate.g4_units('MeV')
Bq = gate.g4_units('Bq')
source = sim.add_source('Generic', 'mysource')
source.mother = 'my_volume'
source.particle = 'proton'
source.activity = 10000 * Bq
source.position.type = 'box'
source.position.dimension = [4 * cm, 4 * cm, 4 * cm]
source.position.translation = [-3 * cm, -3 * cm, -3 * cm]
source.position.rotation = Rotation.from_euler('x', 45, degrees=True).as_matrix()
source.direction.type = 'iso'
source.energy.type = 'gauss'
source.energy.mono = 80 * MeV
source.energy.sigma_gauss = 1 * MeV

All parameters are stored into a dict-like structure (a Box). Particle can be ‘gamma’, ‘e+’, ‘e-’, ‘proton’ (all Geant4 names). The number of particles that will be generated by the source can be described by an activity source.activity = 10 * MBq or by a number of particle source.n = 100. The activity may be automatically decreased according to an exponential decay by setting the half-life source.half_life = 60 * sec. Alternatively, user can provide a TAC (Time Activity Curve) by means of two vectors (times and activities) :

starting_activity = 1000 * Bq
half_life = 2 * sec
times = np.linspace(0, 10, num=500, endpoint=True) * sec
decay = np.log(2) / half_life
activities = [starting_activity * np.exp(-decay * t) for t in times]
source.tac_times = times
source.tac_activities = activities

During the simulation, the activity of this source will be updated according to the current simulation time with a linear interpolation of this TAC. If the simulation time is before the first time or above the last one in the times vector, the activity is considered as zero. The number of elements in the times linspace (here 500) defined the accuracy of the TAC. See example test052.

The positions from were the particles will be generated are defined by a shape (‘box’, ‘sphere’, ‘point’, ‘disc’), defined by several parameters (‘size’, ‘radius’) and orientation (‘rotation’, ‘center’). The direction are defined with ‘iso’, ‘momentum’, ‘focused’. The energy can be defined by a single value (‘mono’) or Gaussian (‘gauss’).

The mother option indicate the coordinate system of the source. By default, it is the world, but it is possible to attach a source to any volume. In that case, the coordinate system of all emitted particles will follow the given volume.

It is possible to indicate a angle_acceptance_volume to the direction of a source. In that case, the particle will be created only if their position & direction make them intersect the given volume. This is for example useful for SPECT imaging in order to limit the particle creation to the ones that will have a chance to reach the detector. Note that the particles that will not intersect the volume will be created anyway but with a zero energy (so not tracked). This mechanism ensures to remain consistent with the required activity and timestamps of the particles, there is no need to scale with the solid angle. See for example test028 test files for more details.

Using direction.type = 'iso', the directions given to primary particles depends on 𝜃 and 𝜙 angles in a spherical coordinate system. By default, 𝜃 varies from 0° to 180° and 𝜙 varies from 0° to 360° (such that any direction is possible). One can define the 𝜃 and 𝜙 ranges (minimum and maximum values) like this:

source.direction.theta = [0, 10 * deg]
source.direction.phi = [0, 90 * deg]

Geant4 defines the direction as:

x = -sin𝜃 cos𝜙;
y = -sin𝜃 sin𝜙;
z = -cos𝜃.

So 𝜃 is the angle in XOZ plane, from -Z to -X; and 𝜙 is the angle in XOY plane from -X to -Y.

Source of ion can be set with the following (see test013)

source1 = sim.add_source('Generic', 'ion1')
source1.particle = 'ion 9 18'  # Fluorine18
source2 = sim.add_source('Generic', 'ion2')
source2.particle = 'ion 53 124'  # Iodine 124

There is some predefined energy spectrum of positron (e+):

source = sim.add_source('Generic', 'Default')
source.particle = 'e+'
source.energy.type = 'F18'  # F18 or Ga68 or C11 ...

It means the positrons will be generated following the (approximated) energy spectrum of the F18 ion. Source code is GateSPSEneDistribution.cpp. Energy spectrum for beta+ emitters are available : F18, Ga68, Zr89, Na22, C11, N13, O15, Rb82. See http://www.lnhb.fr/nuclear-data/module-lara. One example is available in test031.

There is a confine option that allows to generate particles only if their starting position is within a given volume. See phantom_nema_iec_body in the contrib folder. Note that the source volume MUST be larger than the volume it is confined in. Also, note that no particle source will be generated in the daughters of the confine volume.

All options have a default values and can be printed with print(source).

1.4.2. Voxelized sources¶

Voxelized sources can be described as follows:

source = sim.add_source('Voxels', 'vox')
source.particle = 'e-'
source.activity = 4000 * Bq
source.image = 'an_activity_image.mhd'
source.direction.type = 'iso'
source.energy.mono = 100 * keV
source.mother = 'my_volume_name'

This code create a voxelized source. The 3D activity distribution is read from the given image. This image is internally normalized such that the sum of all pixels values is 1, leading to a 3D probability distribution. Particles will be randomly located somewhere in the image according to this probability distribution. Note that once an activity voxel is chosen from this distribution, the location of the particle inside the voxel is performed uniformly. In the given example, 4 kBq of electrons of 140 keV will be generated.

Like all objects, by default, the source is located according to the coordinate system of its mother volume. For example, if the mother volume is a box, it will be the center of the box. If it is a voxelized volume (typically a CT image), it will the center of this image: the image own coordinate system (ITK’s origin) is not considered here. If you want to align a voxelized activity with a CT image that have the same coordinate system you should compute the correct translation. This is done by the function gate.image.get_translation_between_images_center. See the contrib example dose_rate.py.

1.4.3. Phase-Space sources¶

A phase-space source reads particles properties (position, direction, energy, etc.) from a root file and use them as events. Here is an example:

source = sim.add_source("PhaseSpaceSource", "phsp_source")
source.mother = plane.name
source.phsp_file = "input.root"
source.position_key = "PrePositionLocal"
source.direction_key = "PreDirectionLocal"
source.global_flag = False
source.particle = "gamma"
source.batch_size = 4000
source.n = 20000

In that case, the key “PrePositionLocal” in the root tree file will be used to define the position of all generated particles. The flag “global_flag” is False so the position will be relative to the mother volume (the plane here) ; otherwise, position is considered as global (in the world coordinate system).

Limitation: the particle timestamps is NOT read from the phsp and not considered (yet)

The particle type can be set by source.particle = "proton" option. Using this option all generated particles will be for example protons, overriding the particle type specified in the phase space.

Alternatively, by setting source.particle = None the particle type is read from the phase space file using the PDGCode. source.PDGCode_key = PDGCode specifies the name of the entry in the phase space file. Full listing:

source.PDGCode_key = "PDGCode"
source.particle = None

The PDGCode is defined by the particle data group (see https://pdg.lbl.gov/2023/web/viewer.html?file=../reviews/rpp2022-rev-monte-carlo-numbering.pdf). Here is a short overview of common particle types and its corresponding PDG Code

proton: 2212
neutron: 2211
electron: 11
gamma: 22
carbon ion C12: 1000060120

The naming of the phsp file entries generated by e.g. a GATE phase space actor changed over time, most notably from GATE v9 to GATE v10. Setting source.position_key = "PrePositionLocal" will cause the phsp source to look for particle positions in PrePositionLocal_X, PrePositionLocal_Y, PrePositionLocal_Z. source.direction_key = "PreDirectionLocal" will do the corresponding for the particle direction vector components in PreDirectionLocal_X, PreDirectionLocal_y, PreDirectionLocal_Z.

Alternatively, it is possible to directly set the individual components:

source.position_key = None"PrePositionLocal"
source.position_key_x = Position_X"
source.position_key_y = Position_X
source.position_key_z = Position_X
source.direction_key = None
source.direction_key_x = Direction_X
source.direction_key_y = Direction_X
source.direction_key_z = Direction_X

In case of simulating a realistic particle mix, for example the output after a linac, a phsp file could contain a mixture of particles. Typically, one would be interested in simulating a given number of primary particles (e.g. protons), simulating, but not counting as secondary particles (e.g. secondary electrons) in the number of particles to simulate. This can be acchieved by setting generate_until_next_primary = True. Furthermore, the PDG code of the primary particle needs to be specified, as well as a lower energy threshold (in order to identify secondary particles of the same type as the primary particle). The example below will consider protons above 90 MeV as primary particles. Every primary particle found in the phsp file will increase the number of particles simulated counter, while secondary particles (e.g. all other particles in the phsp file) will be simulated, but not be considered in (e.g. not increasing) the number of particles simulated.

source.generate_until_next_primary = True
source.primary_lower_energy_threshold = 90.0 * MeV
source.primary_PDGCode = 2212

For multithread: you need to indicate the entry_start for all threads, as an array, so that each thread starts in the phsp file at a different position. This done for example as follows (see test019_linac_phsp_source_MT.py). Warning, if the phsp reach its end, it will cycle and start back at the beginning.

total_nb_of_particle = 1e6
nb_of_threads = 4
source.entry_start = [total_nb_of_particle * p for p in range(nb_of_threads)]

See all test019 and test060 as examples.

1.4.4. GAN sources (Generative Adversarial Network)¶

(documentation TODO)

test034 : GAN for linac
test038 : GAN for SPECT
test040 : GAN for PET

generate training dataset
train GAN
use GAN as source ; compare to reference

1.4.5. PHID source (Photon from Ion Decay)¶

PHID (Photon from Ion Decay) is a virtual source model that generates photons emitted in the complex decay chain process of alpha-emitter radionuclides, typically for use during simulation of SPECT image acquisition. Given an alpha-emitter radionuclide, the model extracts from Geant4 databases the photon emission lines from all decaying daughters for both isometric transition and atomic relaxation processes. According to a given time range, abundances and activities in the decay chain are considered thanks to the Bateman equations, taking into account the decay rates and the initial abundances. It generates photons with the correct energy and temporal distribution, avoiding the costly Monte Carlo simulation of the complete decay chain. Photons emitted from Bremsstrahlung are ignored, but are not significant for SPECT imaging. Also, the model is not expected to be correct for gammas below 20-30 keV.

See Sarrut et al 2024 Phys. Med. Biol. https://doi.org/10.1088/1361-6560/ad3881

To use such a source, declare a “PhotonFromIonDecaySource” with an ion as particle name, like the “GenericSource”. Only the gammas emitted by atomic relaxation and isomeric transition will be created and tracked. The timing is taken into account by using a TAC (Time Activity Curve) automatically computed from the start and end time of the simulation. The TAC is then binned and the number of bins can be modified. See tests 053.

source = sim.add_source("PhotonFromIonDecaySource", "my_source")
source.particle = f"ion 89 225"
source.position.type = "sphere"
source.position.radius = 1 * nm
source.direction.type = "iso"
source.activity = 10 * kBq
source.atomic_relaxation_flag = True
source.isomeric_transition_flag = True
source.tac_bins = 200
source.dump_log = "phid_log.txt"
source.verbose = True

Also, several command lines tools are provided :

# print information about a radionuclide bi213, pb212, etc.
phid_info ac225

# plot time activity curve of a radionuclide. Options may by set to adapt the timing
phid_tac

# plot gammas lines from a radionuclide (whatever the time)
phid_gammas ac225
phid_atomic_relaxation ac225
phid_isomeric_transition ac225

1.4.6. Pencil Beam sources¶

The Pencil Beam source inherits from the Generic source, and retains therefore the same settings. The main difference consists in the sampling of the position and direction of the particles, which are not sampled independently, but are correlated. In fact, the Pencil Beam source is meant to describe a beam that can converge or diverge. This behaviour is modeled according to the Fermi-Eyges theory (Techniques of Proton Radiotherapy: Transport Theory B. Gottschalk May 1, 2012), that describes the correlated momentum spread of the particle with 4 parameters (each for x and y direction, assuming a beam directed as z):

spot size 𝜎
divergence 𝜃
emittance 𝜀
convergence flag [1,0] The parameters must satisfy the condition:

pi * sigma * theta >= epsilon

The user can set the beam parameters as shown in the example below, for a 120 MeV/n carbon ion beam.

source = sim.add_source("IonPencilBeamSource", "mysource")
source.energy.mono = 1440 * MeV
source.particle = "ion 6 12"  # carbon
source.position.translation = [100 * mm, 0 * mm, 0 * cm]
source.n = 20000
source.direction.partPhSp_x = [
    2.3335754 * mm,
    2.3335754 * mrad,
    0.00078728 * mm * mrad,
    0,
]
source.direction.partPhSp_y = [
    1.96433431 * mm,
    0.00079118 * mrad,
    0.00249161 * mm * mrad,
    0,
]

NOTE: the Pencil Beam source is created by default directed as the positive z axis. To rotate the source, use the source.position.rotation option.

Check all test044 for usage examples.

1.5. Physics¶

The managements of the physic in Geant4 is rich and complex, with hundreds of options. OPENGATE proposes a subset of available options.

1.5.1. Physics list and decay¶

First, the user needs to select a physics list. A physics list contains a large set of predefined physics options, adapted to different problems. Please refer to the Geant4 guide for a detailed explanation. The user can select the physics list with the following:

# Assume that sim is a simulation
phys = sim.get_physics_info()
phys.name = 'QGSP_BERT_EMZ'

The default physics list is QGSP_BERT_EMV. The Geant4 standard physics list are composed of a first part:

FTFP_BERT
FTFP_BERT_TRV
FTFP_BERT_ATL
FTFP_BERT_HP
FTFQGSP_BERT
FTFP_INCLXX
FTFP_INCLXX_HP
FTF_BIC
LBE
QBBC
QGSP_BERT
QGSP_BERT_HP
QGSP_BIC
QGSP_BIC_HP
QGSP_BIC_AllHP
QGSP_FTFP_BERT
QGSP_INCLXX
QGSP_INCLXX_HP
QGS_BIC
Shielding
ShieldingLEND
ShieldingM
NuBeam

And a second part with the electromagnetic interactions:

_EMV
_EMX
_EMY
_EMZ
_LIV
_PEN
__GS
__SS
_EM0
_WVI
__LE

The lists can change according to the Geant4 version (this list is for 10.7).

Moreover, additional physics list are available:

G4EmStandardPhysics_option1
G4EmStandardPhysics_option2
G4EmStandardPhysics_option3
G4EmStandardPhysics_option4
G4EmStandardPhysicsGS
G4EmLowEPPhysics
G4EmLivermorePhysics
G4EmLivermorePolarizedPhysics
G4EmPenelopePhysics
G4EmDNAPhysics
G4OpticalPhysics

Note that EMV, EMX, EMY, EMZ corresponds to option1, 2, 3, 4 (don’t ask us why).

WARNING The decay process, if needed, must be added explicitly. This is done with:

sim.enable_decay(True)
# or
sim.physics_manager = True

Under the hood, this will add two processed to the Geant4 list of processes, G4DecayPhysics and G4RadioactiveDecayPhysics. Those processes are required in particular if decaying generic ion (such as F18) is used as source. Additional information can be found in the following:

1.5.2. Optical Physics Processes¶

1.5.2.1. G4OpticalPhysics physics list¶

To include optical processes in the simulation, explicitly enable them with the following code:

sim.physics_manager.special_physics_constructors.G4OpticalPhysics = True

When G4OpticalPhysics is set to True, the following process are automatically added:

Cerenkov effect
Scintillation
Absorption
Rayleigh scattering
Mie scattering
Wave-length shifting
Boundary Scattering

WARNING: It’s important to note that merely including the G4OpticalPhysics physics list does not automatically activate the Cherenkov process. To generate Cherenkov photons, it’s necessary to set an appropriate electron physics cut in the relevant volume. Currently, setting the electron physics cut to 0.1 mm has been found effective:

sim.physics_manager.set_production_cut("crystal", "electron", 0.1 * mm)

You can find additional details about the G4OpticalPhysics physics list at the following link: https://geant4-userdoc.web.cern.ch/UsersGuides/AllGuides/html/ForApplicationDevelopers/TrackingAndPhysics/physicsProcess.html?highlight=g4opticalphysics#optical-photon-processes

1.5.2.2. Optical Physics Properties¶

The material property table stores the optical properties of materials, where each property is labeled with a name. These properties are of two types: constant properties, which consist of a single value, and property vectors, which are properties varying with the energy of the optical photon. A property vector comprises a series of pairs, each linking a specific energy level with its corresponding value.

To enable Optical physics, material property tables must be stored separately from the material database. This separation allows for easier modification of properties without altering the material database itself. In Gate 10, a default file named OpticalProperties.xml is used, located in the opengate/data folder. Users can specify a custom file by using:

sim.physics_manager.optical_properties_file = PATH_TO_FILE

1.5.2.3. Scintillation¶

A scintillator’s properties are influenced by its photon emission spectrum, which is characterized by an exponential decay process with up to three time constants. The contribution of each component to the total scintillation yield is defined by the parameters SCINTILLATIONYIELD1, SCINTILLATIONYIELD2, and SCINTILLATIONYIELD3. The emission spectra for these decays are specified through the property vectors SCINTILLATIONCOMPONENT1, SCINTILLATIONCOMPONENT2, and SCINTILLATIONCOMPONENT3, in addition to the time constants SCINTILLATIONTIMECONSTANT1, SCINTILLATIONTIMECONSTANT2, and SCINTILLATIONTIMECONSTANT3. These vectors indicate the probability of emitting a photon at a particular energy, and their total should equal one.

To initiate scintillation in a material, the first parameter to be set is SCINTILLATIONYIELD (in units of 1/Mev, 1/keV). This parameter denotes the average number of photons emitted per unit of energy absorbed. The actual photon count follows a normal distribution, with the mean value expressed as:

$$\mu_N = E \cdot \text{SCINTILLATIONYIELD}$$

The standard deviation of this distribution is:

$$\sigma_N = RESOLUTIONSCALE \cdot \sqrt{E \cdot \text{SCINTILLATIONYIELD}}$$

The parameter RESOLUTIONSCALE is derived from the scintillator’s energy resolution, which should exclude any electronic noise influences to reflect the intrinsic energy resolution of the scintillator. It is computed using the following formula:

$$ \text{RESOLUTIONSCALE} = \frac{R}{2.35} \cdot \sqrt{E \cdot \text{SCINTILLATIONYIELD}} $$

In this equation, R stands for the energy resolution (FWHM - Full Width at Half Maximum) at energy E.

<material name="LSO">
  <propertiestable>
    <property name="SCINTILLATIONYIELD" value="26000" unit="1/MeV"/>
    <property name="RESOLUTIONSCALE" value="4.41"/>
    <property name="SCINTILLATIONTIMECONSTANT1" value="40" unit="ns"/>
    <property name="SCINTILLATIONYIELD1" value="1"/>
    <propertyvector name="SCINTILLATIONCOMPONENT1" energyunit="eV">
      <ve energy="2.95167" value="1"/>
    </propertyvector>
    <propertyvector name="ABSLENGTH" unit="m" energyunit="eV">
      <ve energy="1.84" value="50"/>
      <ve energy="4.08" value="50"/>
    </propertyvector>
    <propertyvector name="RINDEX" energyunit="eV">
      <ve energy="1.84" value="1.82"/>
      <ve energy="4.08" value="1.82"/>
    </propertyvector>
  </propertiestable>
</material>

1.5.2.4. Cerenkov photons¶

Cerenkov light emission occurs when a charged particle traverses a dispersive medium at a speed exceeding the medium’s group velocity of light. This emission forms a cone-shaped pattern of photons, with the cone’s opening angle narrowing as the particle decelerates. To simulate Cerenkov optical photon generation in a material, the refractive index must be defined using the RINDEX property of the material.

1.5.2.5. Absorption¶

This process kills the particle. It requires the OpticalProperties.xml properties filled by the user with the Absorption length ABSLENGTH (average distance traveled by a photon before being absorbed by the medium).

1.5.2.6. Mie/Rayleigh Scattering¶

Mie Scattering is a solution derived from Maxwell’s equations for the scattering of optical photons by spherical particles. This phenomenon becomes significant when the radius of the scattering particle is approximately equal to the photon’s wavelength. The formulas for Mie Scattering are complex, and a common simplification used, including in Geant4, is the Henyey-Greenstein (HG) approximation. In cases where the size parameter (diameter of the scattering particle) is small, Mie theory simplifies to the Rayleigh approximation.

For both Rayleigh and Mie scattering, it’s required that the final momentum, initial polarization, and final polarization all lie in the same plane. These processes need the material properties to be defined by the user with specific scattering length data for Mie/Rayleigh scattering, denoted as MIEHG/RAYLEIGH. This represents the average distance a photon travels in a medium before undergoing Mie/Rayleigh scattering. Additionally, for Mie scattering, users must input parameters for the HG approximation: MIEHG_FORWARD (forward anisotropy), MIEHG_BACKWARD (backward anisotropy), and MIEHG_FORWARD_RATIO (the ratio between forward and backward angles). In Geant4, the forward and backward angles can be addressed independently. If the material properties only provide a single value for anisotropy (i.e., the average cosine of the scattering angle), the Materials.xml file might look something like this:

<material name="Biomimic">
  <propertiestable>
   <propertyvector name="ABSLENGTH" unit="cm" energyunit="eV">
     <ve energy="1.97" value="0.926"/>
     <ve energy="2.34" value="0.847"/>
    </propertyvector>
    <propertyvector name="RINDEX" energyunit="eV">
      <ve energy="1.97" value="1.521"/>
      <ve energy="2.34" value="1.521"/>
    </propertyvector>
    <property name="MIEHG_FORWARD" value="0.62" />
    <property name="MIEHG_BACKWARD" value="0.62" />
    <property name="MIEHG_FORWARD_RATIO" value="1.0" />
    <propertyvector name="MIEHG" unit="cm" energyunit="eV">
      <ve energy="1.97" value="0.04"/>
      <ve energy="2.34" value="0.043"/>
    </propertyvector>
  </propertiestable>
</material>

1.5.2.7. Fluorescence¶

Fluorescence involves a three-stage process: Initially, the fluorophore reaches an excited state after absorbing an optical photon from an external source (like a laser or lamp). This excited state typically lasts between 1-10 ns, during which the fluorophore interacts with its surroundings, eventually transitioning to a relaxed-excited state. The final step involves emitting a fluorescent photon, whose energy/wavelength is lower (or wavelength longer) than the excitation photon.

Geant4 models the process of Wavelength Shifting (WLS) in fibers, which are used in high-energy physics experiments. For example, the CMS Hadronic EndCap calorimeter utilizes scintillator tiles integrated with WLS fibers. These fibers absorb the blue light generated in the tiles and re-emit green light to maximize the light reaching the Photomultiplier Tubes (PMTs).

Users of Gate need to specify four properties to define the fluorescent material: RINDEX, WLSABSLENGTH, WLSCOMPONENT, and WLSTIMECONSTANT. WLSABSLENGTH indicates the absorption length of fluorescence, representing the average distance a photon travels before being absorbed by the fluorophore. This distance is typically short, but not zero to prevent immediate photon absorption upon entering the fluorescent material, which would result in fluorescent photons emerging only from the surface. WLSCOMPONENT details the emission spectrum of the fluorescent material, showing the relative intensity at different photon energies, usually derived from experimental measurements. WLSTIMECONSTANT sets the delay between absorption and re-emission of light.

1.5.2.7.1. Simulation of the Fluorescein¶

We define the refractive index of the fluorophore’s environment (water or alcohol):
<material name="Fluorescein">
<propertiestable>
<propertyvector name="RINDEX" energyunit="eV">
<ve energy="1.0" value="1.4"/>
<ve energy="4.13" value="1.4"/>
</propertyvector>

The WLS process encompasses both absorption and emission spectra. If these spectra overlap, a WLS photon might be absorbed and re-emitted repeatedly. To avoid this, one must ensure there is no overlap between these spectra. In the WLS process, there’s no distinction between original photons and WLS photons.

We describe the fluorescein absorption length taken from measurements or literature as function of the photon energy:
<propertyvector name="WLSABSLENGTH" unit="cm" energyunit="eV">
<ve energy="3.19" value="2.81"/>
<ve energy="3.20" value="2.82"/>
<ve energy="3.21" value="2.81"/>
</propertyvector>

We describe the fluorescein Emission spectrum taken from measurements or literature as function of the photon energy:
<propertyvector name="WLSCOMPONENT" energyunit="eV">
    <ve energy="1.771"  value="0.016"/>
    <ve energy="1.850"  value="0.024"/>
    <ve energy="1.901"  value="0.040"/>
    <ve energy="2.003"  value="0.111"/>
    <ve energy="2.073"  value="0.206"/>
    <ve energy="2.141"  value="0.325"/>
    <ve energy="2.171"  value="0.413"/>
    <ve energy="2.210"  value="0.540"/>
    <ve energy="2.250"  value="0.683"/>
    <ve energy="2.343"  value="0.873"/>
    <ve energy="2.384"  value="0.968"/>
    <ve energy="2.484"  value="0.817"/>
    <ve energy="2.749"  value="0.008"/>
    <ve energy="3.099"  value="0.008"/>
</propertyvector>
<property name="WLSTIMECONSTANT" value="1.7" unit="ns"/>
</propertiestable>
</material>

1.5.2.8. Boundary Processes¶

When a photon reaches the boundary between two mediums, its behavior is determined by the characteristics of the materials forming the boundary. If the boundary is between two dielectric materials, the photon’s reaction – whether it undergoes total internal reflection, refraction, or reflection – depends on factors such as the photon’s wavelength, its angle of incidence, and the refractive indices of the materials on either side of the boundary. In contrast, at an interface between a dielectric material and a metal, the photon may either be absorbed by the metal or reflected back into the dielectric material. For simulating a perfectly smooth surface, it’s not necessary for the user to input a G4Surface; the only essential property is the refractive index (RINDEX) of the materials on both sides of the interface. In such cases, Geant4 uses Snell’s Law to compute the probabilities of refraction and reflection.

1.5.3. Defining Surfaces¶

The photon travels through the surface between the two volumes Volume1 and Volume2. To create an optical surface from Volume1 to Volume2, the following command should be used:

sim.physics_manager.add_optical_surface(
  volume_from = "name of volume 1",
  volume_to = "name of volume 2",
  g4_surface_name = "name of the surface between two volumes"
)

The surface between Volume1 and Volume2 is NOT the same surface as that between Volume2 and Volume1; the surface definition is directional. When there is optical transport in both directions, two surfaces should be created. Surface name can be any surface defined in the SurfaceProperties.xml file.

1.5.3.1. EXAMPLES -¶

For instance if the objective is to create a PolishedTeflon_LUT surface from Volume1 to Volume2, and create a different surface like RoughTeflon_LUT from Volume2 to Volume1 -

sim.physics_manager.add_optical_surface(
  volume_from = "name of volume 1",
  volume_to = "name of volume 2",
  g4_surface_name = "PolishedTeflon_LUT"
)

sim.physics_manager.add_optical_surface(
  volume_from = "name of volume 2",
  volume_to = "name of volume 1",
  g4_surface_name = "RoughTeflon_LUT"
)

The keywords volume_from, volume_to, g4_surface_name are just for clarity. You can also define surfaces without them -

sim.physics_manager.add_optical_surface("name of volume 1", "name of volume 2", "PolishedTeflon_LUT")
sim.physics_manager.add_optical_Surface("name of volume 2", "name of volume 1", "PolishedTeflon_LUT")

This creates same surface from Volume1 to Volume2 and from Volume2 to Volume1.

1.5.3.2. LUT Davis Model¶

Available from GATE V8.0 onwards is a model for optical transport called the LUT Davis model [Roncali& Cherry(2013)]. The model is based on measured surface data and allows the user to choose from a list of available surface finishes. Provided are a rough and a polished surface that can be used without reflector, or in combination with a specular reflector (e.g. ESR) or a Lambertian reflector (e.g. Teflon). The specular reflector can be coupled to the crystal with air or optical grease. Teflon tape is wrapped around the crystal with 4 layers.

Surface names of available LUTs -

	BARE	TEFLON	ESR AIR	ESR GREASE
POLISHED	Polished_LUT	PolishedTeflon_LUT	PolishedESR_LUT	PolishedESRGrease_LUT
ROUGH	Rough_LUT	RoughTeflon_LUT	RoughESR_LUT	RoughESRGrease_LUT

The user can extend the list of finishes with custom measured surface data. In GATE, this can be achieved by utilising this tool to calculate LUTs. In the LUT database, typical roughness parameters obtained from the measurements are provided to characterize the type of surface modelled:

ROUGH Ra=0.48 µm, σ=0.57 µm, Rpv=3.12 µm
POLISHED Ra=20.8 nm, σ=26.2 nm, Rpv=34.7 nm

with Ra = average roughness; σ = rms roughness, Rpv = peak-to-valley ratio.

The desired finish should be defined in Surfaces.xml (file available in https://github.com/OpenGATE/GateContrib/tree/master/imaging/LUTDavisModel):

 <surface model="DAVIS" name="RoughTeflon_LUT" type="dielectric_LUTDAVIS" finish="RoughTeflon_LUT">
 </surface>

The detector surface, called Detector_LUT, defines a polished surface coupled to a photodetector with optical grease or a glass interface (similar index of refraction 1.5). Any surface can be used as a detector surface when the Efficiency is set according to the following example:

 <surface model="DAVIS" name="**Detector_LUT**" type="dielectric_LUTDAVIS" finish="Detector_LUT">
     <propertiestable>
      <propertyvector name="**EFFICIENCY**" energyunit="eV">
        <ve energy="1.84" value="**1**"/>
        <ve energy="4.08" value="**1**"/>
      </propertyvector>
    </propertiestable>
  </surface>

Running the simulation produces an output in the terminal confirming that the LUT data is read in correctly. The user should check the presence of these lines in the terminal. For example:

===== XML PATH ====: ./Surfaces.xml
===== XML PATH ====: ...
LUT DAVIS - data file: .../Rough_LUT.dat read in!
Reflectivity LUT DAVIS - data file: .../Rough_LUTR.dat read in!
===== XML PATH ====: ./Surfaces.xml
===== XML PATH ====: ...
LUT DAVIS - data file: .../Detector_LUT.dat read in!
Reflectivity LUT DAVIS - data file: .../Detector_LUTR.dat read in!

1.5.3.2.1. Detection of Optical Photons¶

Once the simulation is finished, the optical photon data can be found in the Hits Tree in the ROOT output. The Hits Tree consists of events that ended their path in the geometry defined as the sensitive detector (SD). Thus, photons can either be detected or absorbed in the crystal material when set as SD. The user can identify the optical photons from other particles using the PDGEncoding (-22 for optical photons).

NOTE - From Geant4 10.7, PDG code for optical photon has changed from 0 (zero) to -22.

1.5.3.2.2. Example¶

The example (https://github.com/OpenGATE/GateContrib/tree/master/imaging/LUTDavisModel) includes a 3 mm x 3 mm x 20 mm scintillation crystal coupled to a 3 mm x 3 mm detector area. The source is positioned at the side of the crystal, irradiating it at 10 mm depth. The set surface is RoughTeflon_LUT in combination with the Detector_LUT as the photo detector surface.

1.5.3.2.3. Background¶

The crystal topography is obtained with atomic force microscopy (AFM). From the AFM data, the probability of reflection (1) and the reflection directions (2) are computationally determined, for incidence angles ranging from 0° to 90°. Each LUT is computed for a given surface and reflector configuration. The reflection probability in the LUT combines two cases: directly reflected photons from the crystal surface and photons that are transmitted to the reflector surface and later re-enter the crystal. The key operations of the reflection process are the following: The angle between the incident photon (Old Momentum) and the surface normal are calculated. The probability of reflection is extracted from the first LUT. A Bernoulli test determines whether the photon is reflected or transmitted. In case of reflection two angles are drawn from the reflection direction LUT.

Old Momentum to New Momentum. The old momentum is the unit vector that describes the incident photon. The reflected/transmitted photon is the New Momentum described by two angles φ, 𝛳.

1.5.3.3. UNIFIED Model¶

The UNIFIED model allows the user to control the radiant intensity of the surface: Specular lobe, Specular spike, Backscatter spike (enhanced on very rough surfaces) and Reflectivity (Lambertian or diffuse distribution). The sum of the four constants is constrained to unity. In that model, the micro-facet normal vectors follow a Gaussian distribution defined by sigmaalpha ($\sigma_a$) given in degrees. This parameter defines the standard deviation of the Gaussian distribution of micro-facets around the average surface normal. In the case of a perfectly polished surface, the normal used by the G4BoundaryProcess is the normal to the surface.

An example of a surface definition looks like:

<surface name="rough_teflon_wrapped" type="dielectric_dielectric" sigmaalpha="0.1" finish="groundbackpainted">
 <propertiestable>
   <propertyvector name="SPECULARLOBECONSTANT" energyunit="eV">
     <ve energy="4.08" value="1"/>
     <ve energy="1.84" value="1"/>
   </propertyvector>
   <propertyvector name="RINDEX" energyunit="eV">
     <ve energy="4.08" value="1"/>
     <ve energy="1.84" value="1"/>
   </propertyvector>
   <propertyvector name="REFLECTIVITY" energyunit="eV">
     <ve energy="1.84" value="0.95"/>
     <ve energy="4.08" value="0.95"/>
   </propertyvector>
   <propertyvector name="EFFICIENCY" energyunit="eV">
     <ve energy="1.84" value="0"/>
     <ve energy="4.08" value="0"/>
   </propertyvector>
 </propertiestable>
</surface>

The attribute type can be either dielectric_dielectric or dielectric_metal, to model either a surface between two dielectrica or between a dielectricum and a metal. The attribute sigma-alpha models the surface roughness and is discussed in the next section. The attribute finish can have one of the following values: ground, polished, ground-back-painted, polished-back-painted, ground-front-painted and polished-front-painted. It is therefore possible to cover the surface on the inside or outside with a coating that reflects optical photons using Lambertian reflection. In case the finish of the surface is polished, the surface normal is used to calculate the probability of reflection. In case the finish of the surface is ground, the surface is modeled as consisting of small micro-facets. When an optical photon reaches a surface, a random angle $\sigma$ is drawn for the micro facet that is hit by the optical photon. Using the angle of incidence of the optical photon with respect to this micro facet and the refractive indices of the two media, the probability of reflection is calculated.

In case the optical photon is reflected, four kinds of reflection are possible. The probabilities of the first three are given by the following three property vectors:

SPECULARSPIKECONSTANT gives the probability of specular reflection about the average surface normal
SPECULARLOBECONSTANT gives the probability of specular reflection about the surface normal of the micro facet
BACKSCATTERCONSTANT gives the probability of reflection in the direction the optical photon came from

LAMBERTIAN (diffuse) reflection occurs when none of the other three types of reflection happens. The probability of Lambertian reflection is thus given by one minus the sum of the other three constants.

When the photon is refracted, the angle of refraction is calculated from the surface normal (of the average surface for polished and of the micro facet for rough) and the refractive indices of the two media.

When an optical photon reaches a painted layer, the probability of reflection is given by the property vector REFLECTIVITY. In case the paint is on the inside of the surface, the refractive indices of the media are ignored, and when the photon is reflected, it undergoes Lambertian reflection.

When the paint is on the outside of the surface, whether the photon is reflected on the interface between the two media is calculated first, using the method described in the previous section. However, in this case the refractive index given by the property vector RINDEX of the surface is used. When the photon is refracted, it is reflected using Lambertian reflection with a probability REFLECTIVITY. It then again has to pass the boundary between the two media. For this, the method described in the previous section is used again and again, until the photon is eventually reflected back into the first medium or is absorbed by the paint.

A dielectric_dielectric surface may have a wavelength dependent property TRANSMITTANCE. If this is specified for a surface it overwrites the Snell’s law’s probability. This allows the simulation of anti-reflective coatings.

1.5.3.3.1. Detection of optical photons¶

Optical photons can be detected by using a dielectric-metal boundary. In that case, the probability of reflection should be given by the REFLECTIVITY property vector. When the optical photon is reflected, the UNIFIED model is used to determine the reflection angle. When it is absorbed, it is possible to detect it. The property vector EFFICIENCY gives the probability of detecting a photon given its energy and can therefore be considered to give the internal quantum efficiency. Note that many measurements of the quantum efficiency give the external quantum efficiency, which includes the reflection: external quantum efficiency = efficiency*(1-reflectivity).

The hits generated by the detection of the optical photons are generated in the volume from which the optical photons reached the surface. This volume should therefore be a sensitive detector.

1.5.4. Electromagnetic parameters¶

WARNING : this part is work in progress. DO NOT USE YET.

Electromagnetic parameters are managed by a specific Geant4 object called G4EmParameters. It is available with the following:

phys = sim.get_physics_info()
em = phys.g4_em_parameters
em.SetFluo(True)
em.SetAuger(True)
em.SetAugerCascade(True)
em.SetPixe(True)
em.SetDeexActiveRegion('world', True, True, True)

WARNING: it must be set after the initialization (after sim.initialize() and before output = sim.start()).

The complete description is available in this page: https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html

1.5.5. Managing the cuts and limits¶

WARNING : this part is work in progress. DO NOT USE YET.

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/thresholdVScut.html

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/cutsPerRegion.html

https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/userLimits.html

1.6. Actors and Filters¶

The “Actors” are scorers can store information during simulation such as dose map or phase-space (like a “tally” in MCNPX). They can also be used to modify the behavior of a simulation, such as the MotionActor that allows to move volumes, this is why they are called “actor”.

1.6.1. SimulationStatisticsActor¶

The SimulationStatisticsActor actor is a very basic tool that allow to count the number of runs, events, tracks and steps that have been created during a simulation. Most of the simulations should include this actor as it gives valuable information. Once the simulation ends, the user can retrieve the values as follows:

# during the initialisation:
stats = sim.add_actor('SimulationStatisticsActor', 'Stats')
stats.track_types_flag = True

# (...)
output = sim.start()

# after the end of the simulation
stats = output.get_actor('Stats')
print(stats)
stats.write('myfile.txt')

The stats object contains the counts dictionary that contains all numbers. In addition, if the flag track_types_flag is enabled, the stats.counts.track_types will contain a dictionary structure with all types of particles that have been created during the simulation. The start and end time of the whole simulation is also available. Speeds are also estimated (primary per sec, track per sec and step per sec). You can write all the data to a file like in the previous GATE, via stats.write. See source.

1.6.2. DoseActor¶

The DoseActor computes a 3D edep/dose map for deposited energy/absorbed dose in a given volume. The dose map is a 3D matrix parameterized with: dimension (number of voxels), spacing (voxel size), and translation (according to the coordinate system of the “attachedTo” volume). There is no possibility to rotate this 3D matrix for the moment. By default, the matrix is centered according to the volume center.

Like any image, the output dose map will have an origin. By default, it will consider the coordinate system of the volume it is attached to, so at the center of the image volume. The user can manually change the output origin, using the option output_origin of the DoseActor. Alternatively, if the option img_coord_system is set to True the final output origin will be automatically computed from the image the DoseActor is attached to. This option calls the function get_origin_wrt_images_g4_position to compute the origin. See the figure for details.

Several tests depict usage of DoseActor: test008, test009, test021, test035, etc.

dose = sim.add_actor("DoseActor", "dose")
dose.output = output_path / "test008-edep.mhd"
dose.mother = "waterbox"
dose.size = [99, 99, 99]
mm = gate.g4_units.mm
dose.spacing = [2 * mm, 2 * mm, 2 * mm]
dose.translation = [2 * mm, 3 * mm, -2 * mm]
dose.uncertainty = True
dose.hit_type = "random"

1.6.3. PhaseSpaceActor¶

A PhaseSpaceActor stores any set of particles reaching a given volume during the simulation. The list of attributes that are kept for each stored particle can be specified by the user.

phsp = sim.add_actor("PhaseSpaceActor", "PhaseSpace")
phsp.mother = plane.name
phsp.attributes = [
    "KineticEnergy",
    "Weight",
    "PostPosition",
    "PrePosition",
    "ParticleName",
    "PreDirection",
    "PostDirection",
    "TimeFromBeginOfEvent",
    "GlobalTime",
    "LocalTime",
    "EventPosition",
]
phsp.output = "test019_hits.root"
f = sim.add_filter("ParticleFilter", "f")
f.particle = "gamma"
phsp.filters.append(f)

In this example, the PhaseSpace will store all particles reaching the given plane. For each particle, some information will be stored, as shown in the attributes array: energy, position, name, time, etc. The list of available attributes name can be seen in the file : GateDigiAttributeList.cpp. Here is the current list:

TotalEnergyDeposit
PostKineticEnergy PreKineticEnergy KineticEnergy TrackVertexKineticEnergy EventKineticEnergy
LocalTime GlobalTime TimeFromBeginOfEvent TrackProperTime
Weight
TrackID ParentID EventID RunID ThreadID
TrackCreatorProcess ProcessDefinedStep
ParticleName
TrackVolumeName TrackVolumeCopyNo
PreStepVolumeCopyNo PostStepVolumeCopyNo TrackVolumeInstanceID
PreStepUniqueVolumeID PostStepUniqueVolumeID HitUniqueVolumeID
Position PostPosition PrePosition EventPosition TrackVertexPosition
Direction PostDirection PreDirection PreDirectionLocal TrackVertexMomentumDirection EventDirection

The output is a root file that contains a tree. It can be analysed for example with uproot.

1.6.4. Hits-related actors (digitizer)¶

In legacy Gate, the digitizer module is a set of tools used to simulate the behaviour of the scanner detectors and signal processing chain. The tools consider a list of interactions occurring in the detector (e.g. in the crystal), named as “hits collections”. Then, this collection of hits is processed and filtered by different modules to end up with a final digital value. To start a digitizer chain, we must start defining a HitsCollectionActor, explained in the next sections.

1.6.4.1. DigitizerHitsCollectionActor¶

The DigitizerHitsCollectionActor is an actor that collects hits occurring in a given volume (or one of its daughters). Every time a step occurs in the volume a list of attributes is recorded. The list of attributes is defined by the user as follows:

hc = sim.add_actor('DigitizerHitsCollectionActor', 'Hits')
hc.mother = ['crystal1', 'crystal2']
hc.output = 'test_hits.root'
hc.attributes = ['TotalEnergyDeposit', 'KineticEnergy', 'PostPosition',
                 'CreatorProcess', 'GlobalTime', 'VolumeName', 'RunID', 'ThreadID', 'TrackID']

In this example, the actor is attached (mother option) to several volumes (crystal1 and crystal2 ) but most of the time, one single volume is sufficient. This volume is important: every time an interaction (a step) is occurring in this volume, a hit will be created. The list of attributes is defined with the given array of attribute names. The names of the attributes are as close as possible to the Geant4 terminology. They can be of a few types: 3 (ThreeVector), D (double), S (string), I (int), U (unique volume ID, see DigitizerAdderActor section). The list of available attributes is defined in the file core/opengate_core/opengate_lib/GateDigiAttributeList.cpp and can be printed with:

import opengate_core as gate_core
am = gate_core.GateDigiAttributeManager.GetInstance()
print(am.GetAvailableDigiAttributeNames())

    Direction 3
    EventDirection 3
    EventID I
    EventKineticEnergy D
    EventPosition 3
    GlobalTime D
    HitUniqueVolumeID U
    KineticEnergy D
    LocalTime D
    ParticleName S
    Position 3
    PostDirection 3
    PostKineticEnergy D
    PostPosition 3
    PostStepUniqueVolumeID U
    PostStepVolumeCopyNo I
    PreDirection 3
    PreDirectionLocal 3
    PreKineticEnergy D
    PrePosition 3
    PreStepUniqueVolumeID U
    PreStepVolumeCopyNo I
    ProcessDefinedStep S
    RunID I
    ThreadID I
    TimeFromBeginOfEvent D
    TotalEnergyDeposit D
    TrackCreatorProcess S
    TrackID I
    TrackProperTime D
    TrackVertexKineticEnergy D
    TrackVertexMomentumDirection 3
    TrackVertexPosition 3
    TrackVolumeCopyNo I
    TrackVolumeInstanceID I
    TrackVolumeName S
    Weight D

Warning: KineticEnergy, Position and Direction are available for PreStep and for PostStep, and there is a “default” version corresponding to the legacy Gate (9.X).

Pre version	Post version	default version
PreKineticEnergy	PostKineticEnergy	KineticEnergy (Pre)
PrePosition	PostPosition	Position (Post)
PreDirection	PostDirection	Direction (Post)

Attributes correspondence with Gate 9.X for Hits and Singles:

Gate 9.X	Gate 10
edep or energy	TotalEnergyDeposit
posX/Y/Z of globalPosX/Y/Z	PostPosition_X/Y/Z
time	GlobalTime

At the end of the simulation, the list of hits can be written as a root file and/or used by subsequent digitizer modules (see next sections). The Root output is optional, if the output name is None nothing will be written. Note that, like in Gate, every hit with zero deposited energy is ignored. If you need them, you should probably use a PhaseSpaceActor. Several tests using DigitizerHitsCollectionActor are proposed: test025, test028, test035, etc.

The two actors used to convert some hits to one digi are “DigitizerHitsAdderActor” and “DigitizerReadoutActor” described in the next sections and illustrated in the figure:

1.6.4.2. DigitizerHitsAdderActor¶

This actor groups the hits per different volumes according to the option group_volume (by default, this is the deeper volume that contains the hit). All hits (in the same event) occurring in the same volume are gathered into one single digi according to one of the two available policies:

EnergyWeightedCentroidPosition:
- the final energy (“TotalEnergyDeposit”) is the sum of all deposited energy
- the position (“PostPosition”) is the energy-weighted centroid position
- the time (“GlobalTime”) is the time of the earliest hit
EnergyWinnerPosition
- the final energy (“TotalEnergyDeposit”) is the energy of the hit with the largest deposited energy
- the position (“PostPosition”) is the position of the hit with the largest deposited energy
- the time (“GlobalTime”) is the time of the earliest hit

sc = sim.add_actor("DigitizerAdderActor", "Singles")
sc.output = 'test_hits.root'
sc.input_digi_collection = "Hits"
sc.policy = "EnergyWeightedCentroidPosition"
# sc.policy = "EnergyWinnerPosition"
sc.group_volume = crystal.name

Note that this actor is only triggered at the end of an event, so the mother volume to which it is attached has no effect. Examples are available in test 037.

1.6.4.3. DigitizerReadoutActor¶

This actor is the same as the previous one (DigitizerHitsAdderActor) with one additional option: the resulting positions of the digi are set in the center of the defined volumes (discretized). We keep two different actors (Adder and Readout) to be close to the previous legacy GATE versions. The additional option discretize_volume indicates the volume name in which the discrete position will be taken.

sc = sim.add_actor("HitsReadoutActor", "Singles")
sc.input_digi_collection = "Hits"
sc.group_volume = stack.name
sc.discretize_volume = crystal.name
sc.policy = "EnergyWeightedCentroidPosition"

Examples are available in test 037.

1.6.4.4. DigitizerGaussianBlurringActor¶

Digitizer module for blurring an attribute such as the time or the energy (single value only, not a vector). The blurring method can be “Gaussian”, “InverseSquare” or “Linear” :

For Gaussian: the sigma or the FWHM should be given blur_sigma or .blur_fwhm options.

For InverseSquare: blur_reference_value and blur_reference_value EQUATION

For Linear: blur_reference_value, blur_reference_value and blur_slope EQUATION

bc = sim.add_actor("DigitizerBlurringActor", "Singles_with_blur")
bc.output = "output.root"
bc.input_digi_collection = "Singles_readout"
bc.blur_attribute = "GlobalTime"
bc.blur_method = "Gaussian"
bc.blur_fwhm = 100 * ns

1.6.4.5. DigitizerSpatialBlurringActor¶

(documentation TODO)

warning: if blur leads to point outside volume (keep_in_solid_limits option). Useful for monocrystal. Should probably not be used for pixelated crystal.

1.6.4.6. DigitizerEnergyWindowsActor¶

(documentation TODO) for spect, test028

1.6.4.7. DigitizerProjectionActor¶

(documentation TODO) for spect, test028

1.6.4.8. DigitizerEfficiencyActor¶

Digitizer module for simulating detection with lowered (non-100%) efficiency. It is initialized with a float number representing the probability of storing the hit. An efficiency of 0 means that absolutely nothing will be recorded, whereas a probability of 1 means that all digis will be stored.

For each digi the digitizer receives, a random float is uniformly shot between 0 and 1: if this random float is below the efficiency set when initializing the digitizer, then the digit is kept, otherwise, it is discarded.

As an example, a DigitizerEfficiencyActor with an efficiency of 0.3 can be initialized as follows:

ea = sim.add_actor("DigitizerEfficiencyActor", "Efficiency")
ea.input_digi_collection = "Hits"
ea.efficiency = 0.3

A more detailed example can be found in test 57.

1.6.5. Coincidences Sorter¶

Please, be aware that the current version of the Coincidence sorter is still work in progress. Coincidence sorter is only offline yet.

The Coincidence Sorter searches, into the singles list, for pairs of coincident singles. Whenever two or more singles are found within a coincidence time window, these singles are grouped to form a Coincidence event.

As an example, a Coincidence Sorter is shown here:

singles_tree = root_file["Singles_crystal"]
...
ns = gate.g4_units.nanosecond
time_window = 3 * ns
policy="keepAll"

minSecDiff=1 #NOT YET IMPLEMENTED
# apply coincidences sorter
coincidences = coincidences_sorter(singles_tree, time_window, minSecDiff, policy, chunk_size=1000000)

As parameters Coincidence Sorter expects as input:

Singles Tree
Defined coincidence time window
Minimum sector difference or minimal distance between the detectors triggered the coincidence needed for removing geometrically impossible coincidences (Not yet implemented),
Policy to process the multiple coincidences. When more than two singles are found in coincidence, several types of behavior could be implemented.
Chunk size important for very large root files to avoid loading everything in memory.

GATE allows to model so far 2 different policy rules that can be used in such a case:

Policy	Description	Equivalent in Gate 9.X
keepAll	Each good pairs are considered	takeAllGoods
removeMultiples	No multiple coincidences are accepted, no matter how many good pairs are present	killAll

A more detailed example can be found in test 72.

1.6.6. MotionVolumeActor¶

(documentation TODO) test029

1.6.7. ARFActor (and ARFTrainingDatasetActor)¶

(documentation TODO) test043 The detector MUST be oriented such that the depth is Z dimension

1.6.8. LETActor¶

(documentation TODO) test050

1.6.9. ComptonSplittingActor¶

The Compton splitting actor generates N particles, each with a weight equal to the initial track weight divided by N, whenever a Compton process occurs. To tailor the splitting process to your specific application, you can use various options, as presented in test 71 :

compt_splitting_actor = sim.add_actor("ComptSplittingActor", "ComptSplitting")
compt_splitting_actor.mother = W_tubs.name
compt_splitting_actor.splitting_factor = nb_split
compt_splitting_actor.russian_roulette = True
compt_splitting_actor.rotation_vector_director = True
compt_splitting_actor.vector_director = [0, 0, -1]

The options include:

the splitting Number: Specifies the number of splits to create.
A Russian Roulette to activate : Enables selective elimination based on a user-defined angle, with a probability of 1/N.
A Minimum Track Weight: Determines the minimum weight a track must possess before undergoing subsequent Compton splitting. To mitigate variance fluctuations or too low-weight particles, I recommend to set the minimum weight to the average weight of your track multiplied by 1/N², with N depending on your application.

1.7. Additional functionalities¶

The command line tool opengate_info print information about the current installation (Geant4 version, ITK version etc).

The command line tool opengate_tests run all GATE tests. With the option -r only the last 10 tests and 1/4 of the remaining tests are run. With the option -i XX, it run the tests from XX. Each test dump log in the tests/log folder.

The command line tool opengate_user_info allows you to print all default and possible parameters for Volumes, Sources, Physics and Actors elements. This is verbose, but it allows you to have a dynamic documentation of everything currently available in the installed gate version.

1.8. Contributions¶

The “contribution folder” contains additional useful functions that do not belong to the core of gate.

1.8.1. Dose rate computation¶

(documentation TODO), test035

1.8.2. Phantoms¶

1.8.2.1. Phantom: IEC 6 spheres NEMA phantom¶

An analytical model of the 6 spheres IEC NEMA phantom is provided. It can be used as follows:

import opengate as gate
import opengate.contrib.phantoms.nemaiec as gate_iec

sim = gate.Simulation()
iec_phantom = gate_iec.add_iec_phantom(sim, 'iec_phantom')
activities = [3 * BqmL, 4 * BqmL, 5 * BqmL, 6 * BqmL, 9 * BqmL, 12 * BqmL]
iec_source = gate_iec.add_spheres_sources(sim, 'iec_phantom', 'iec_source', 'all', activities)
iec_bg_source = gate_iec.add_background_source(sim, 'iec_phantom', 'iec_bg_source', 0.1 * BqmL)

The rotation should be adapted according to your need. The order of the 6 spheres can be changed with the parameter sphere_starting_angle of the add_iec_phantom command.

Example can be found in test015 (and others).

1.8.2.2. Phantom: cylinder phantom for PET NECR¶

An analytical model of the simple NECR phantom (cylinder and linear source) is provided. It can be used as follows:

import opengate as gate
import opengate.contrib.phantoms.necr as gate_necr

sim = gate.Simulation()
necr_phantom = gate_necr.add_necr_phantom(sim, 'necr_phantom')
necr_source = gate_necr.add_necr_source(sim, 'necr_phantom')
necr_source.activity = 1000 * Bq

Example can be found in test049 (and others).

1.8.3. Radiation therapy linac¶

Important Notice: Please be aware that the models provided within the OpenGate toolkit are based on approximate simulations. Users are strongly encouraged to independently verify these models against empirical data to ensure their applicability and accuracy for specific use cases.

The following models are available:

Elekta Synergy, without patient specific collimation systems
Elekta Versa HD, with Agility multileaf collimator (160 leaves) and back-up jaws.

import opengate as gate
import opengate.contrib.linacs.elektasynergy as synergy
import opengate.contrib.linacs.elektaversa as versa

sim = gate.Simulation()
linac1 = synergy.add_linac(sim)
linac2 = versa.add_linac(sim)

1.8.3.1. LINACs reference systems :¶

Each LINAC head is simulated with a z-axis translation relative to the world center. This translation aligns the machine’s isocenter with the world center, with a user-defined Source-Axis Distance (SAD). The “translation_from_sad” function (example in test019_linac_elekta_versa.py) can be used to move the LINAC head with a translation relative to the SAD.

The “rotation_around_user_point” function enables LINAC head rotation around either the world center (i.e the isocenter) or a user-defined point. Axis and angle lists for each axis must be defined in a way consistent with Rotation_from_euler. An example illustrating how to use this function is available in test019_elekta_versa.py.

1.8.3.2. LINACs source types :¶

There are two way to simulate the LINAC photon beam: from an electron source impinging the tungsten target, and from a phase space source provided by the user (only .root extension)

1.8.3.2.1. Electron source¶

The function to create an electron is add_electron_source (test019_linac_elekta_versa.py). The electron source replicates the processes undergone by electrons to generate the photon beam. It’s positioned at the beginnin of the linac target and characterized as follows:

Electron momentum is strictly aligned along the z-axis.
Electron beam kinetic energy is 6.7 MeV for the LINAC Synergy and 6.4 MeV for the LINAC Versa.
Initial X-Y positions of electrons are defined as a 2 mm radius cylinder for the LINAC Synergy and as a 2D Gaussian distribution with a σ of 0.47 mm for the LINAC Versa. These parameters delineate the profile of the photon beam derived from the electron beam.

1.8.3.2.2. Phase space source¶

When simulating electron sources becomes too time-consuming, you can switch to using a phase space source. This source have to be put before your custom collimation system for patients. To get started, you’ll need a rootfile. The user is free to use its own rootfile, tailored to his photon beam setup.

The prerequisite parameters to include into your rootfile are :

Position in X, Y, Z
Direction in dX, dY, dZ
Energy
Weight (important if you’ve used any biasing methods)

The photon polarization is not available with the phase space source. If a user wants to build a phase space source encompassing a LINAC-specific photon beam, He can run the simulation with the electron source and use the phase_space_plan function to generate the required phase space.

1.8.3.3. Phase space plan¶

The phase space is generated through two sequential steps: firstly, employing the function “add_phase_space_plane” to establish a plane, and secondly, attaching a phase space actor to this designated plane using “add_phase_space.” Examples illustrating the utilization of these functions are provided in test019_linac_elekta_versa.py. This plane is positioned a few centimeters beyond the mirror, before the collimation system. These functions offer a convenient way to generate a phase space, subsequently usable as a source.

1.8.3.4. Empty LINAC box¶

In case users want to simulate a LINAC box without shaping and/or collimation systems, the function add_empty_linac_box is available and the created object behaves the same way as the normal LINAC head.

1.8.3.5. Collimation systems¶

1.8.3.5.1. Multi-Leaf Collimator¶

The LINAC Versa module includes a simulation model of the Agility multileaf collimator, featuring two banks of 80 leaves each. In the simulation, each leaf tip is rounded to resemble real leaf characteristics, and leaves operate independently. To use this collimator system, users can simply employ the “add_mlc” function, which sets up the multileaf collimator in a closed position. An example is provided in test019_linac_elekta_versa_with_mlc.py.

For users who wish to manually control specific leaf movements, the “field” function (elektaversa.py) is available. An example within the “set_rectangular_field” function demonstrates its usage. This function requires parameters including the mlc and jaw objects, as well as vectors specifying the real (not the projected positions at the isocenter) x-axis positions for the leaves within the multileaf collimator and y-axis positions for the jaws.

1.8.3.5.2. Back-up Jaws¶

The LINAC Versa module also includes a complex model of the Elekta Versa Y-jaws. Users can simply employ the “add_jaws” function, setting up the jaws in a closed position. The y-axis jaws position can also be modified using the “field” function.

1.8.3.5.3. Rectangular Fields¶

For a quick and easy establishment of a rectangular field, for example to compare measurements between experiments and simulations, users can utilize the “set_rectangular_field” function. An example is available in test019_linac_elekta_versa_with_mlc.py.

To use it, users just provide the MLC and jaws objects he has created, along with the desired x and y dimensions of the field at the isocenter. The LINAC head position relative to the isocenter needs to be specified to adjust the leaves and jaws accurately.

1.8.3.6. Reproduction of a radiotherapy plan irradiation¶

The LINAC Versa HD module, combined to the dicomrtplan module, enable the reproduction of an IMRT or VMAT irradiation sequence. By default, the LINAC head position is set-up at a source-axis distance of 100 cm. The rotation will be performed around the center of the world.

1.8.3.6.1. Lecture of the DICOM RT plan¶

The dicomrtplan.py module can be used to get useful data from DICOM radiotherapy plans for reproducing radiotherapy sessions. The “read” function returns a Python dictionary with lists of leaf positions, y-jaw positions, normalized monitor units, and isocenter for each LINAC head rotation:

import opengate as gate
import opengate.contrib.linacs.elektaversa as versa
import opengate.contrib.linacs.dicomrtplan as rtplan

paths = utility.get_default_test_paths(__file__, output_folder="test019_linac")
sim = gate.Simulation()
linac = versa.add_linac(sim)

rt_plan_parameters = rtplan.read(str(paths.data / "DICOM_RT_plan.dcm"))
leaves_positions = rt_plan_parameters["leaves"]
jaws_1_positions = rt_plan_parameters["jaws 1"]
jaws_2_positions = rt_plan_parameters["jaws 2"]
linac_head_angles = rt_plan_parameters["gantry angle"]
monitor_unit_weights = rt_plan_parameters["weight"]
isocenter = rt_plan_parameters["isocenter"]

All parameters related to modifying the LINAC head geometry are linearly interpolated between two subsequent control points. The deposited dose between these control points is assumed to be delivered from a single rotation angle, also interpolated between them. Note that the isocenter remains the same for all LINAC head rotations, but it’s returned as a list for convenience.

1.8.3.6.2. Simulation according to the DICOM RT plan¶

1.8.3.6.2.1. LINAC head rotation¶

To simulate a VMAT or IMRT session, we opted to realise one GEANT4 run for each LINAC head rotation. Motion actor are used to adjust the geometry based on the DICOM RT plan between each run. The “set_linac_head_motion” function links relevant parameters extracted from the DICOM RT plan with the LINAC head and the collimation system to be moved.

It’s important to note that the rotation is performed around the center of the simulation. Consequently, the patient geometry must be relocated to align the isocenter with the world center. An illustrative example is provided in test019_linac_elekta_versa_rt_plan_isocenter.py.

1.8.3.6.2.2. Modulation of the number of particle to send¶

To tailor the number of particles sent for each LINAC head rotation, a weight is applied to the “activity” specified by the user. This weight is extracted from the DICOM RT plan and is normalized relative to the median weight of the specific therapy.

The function “set_time_intervals_from_rtplan” (with an example provided in test019_linac_elekta_versa_with_mlc_rt_plan.py) converts the weight into irradiation time. Consequently, this time becomes directly proportional to the number of particles to be sent.

1.8.4. SPECT imaging systems¶

Important Notice: Please be aware that the models provided within the OpenGate toolkit are based on approximate simulations. Users are strongly encouraged to independently verify these models against empirical data to ensure their applicability and accuracy for specific use cases.

The following models are available:

GE Discovery 670 SPECT
Siemens Symbia Intevo Bold SPECT

import opengate as gate
import opengate.contrib.spect.ge_discovery_nm670 as discovery
import opengate.contrib.spect.siemens_intevo as intevo

sim = gate.Simulation()

spect = discovery.add_spect_head(sim, "discovery1", collimator_type="melp")
crystal = sim.volume_manager.get_volume(f"{spect.name}_crystal")
discovery.add_digitizer_tc99m(sim, crystal.name, "digit_tc99m")

spect = discovery.add_spect_head(sim, "discovery12", collimator_type="lehr")
crystal = sim.volume_manager.get_volume(f"{spect.name}_crystal")
discovery.add_digitizer_lu177(sim, crystal.name, "digit_lu177")

spect = intevo.add_spect_head(sim, "intevo1", collimator_type="melp")
crystal = sim.volume_manager.get_volume(f"{spect.name}_crystal")
intevo.add_digitizer_tc99m(sim, crystal.name, "digit_tc99m")

spect = discovery.add_spect_head(sim, "intevo2", collimator_type="lehr")
crystal = sim.volume_manager.get_volume(f"{spect.name}_crystal")
intevo.add_digitizer_lu177(sim, crystal.name, "digit_lu177")

test028

1.8.5. PET imaging systems¶

Important Notice: Please be aware that the models provided within the OpenGate toolkit are based on approximate simulations. Users are strongly encouraged to independently verify these models against empirical data to ensure their applicability and accuracy for specific use cases.

The following models are available:

Philips Vereos Digital PET
Siemens Biograph Vision PET

test037