Geant4.jl
Julia bindings for the Geant4 particle transportation toolkit. It is using CxxWrap.jl package to wrap C++ types and functions to Julia. Since the Geant4 toolkit is rather large and complex, writing the wrapper code by hand is not really an option. For this we use the package WrapIt that automates the generation of the wrapper code making use of the clang library.
Documentation of the concepts and how to write applications with the Geant4 toolkit can be found with the Application Developer Guide or the Classes and Members reference guide for a detailed description of each C++ class. In this document we will only highlight the differences between the Julia and the C++ API. We will document the additional types that have been added on top of the C++ classes to make the user interface more Julia friendly. To distinguish these new types from the types coming directly from C++ via the CxxWrap wrappers, these types are prefixed with G4JL
.
Installation
The Geant4.jl package does no require any special installation. Stable releases are registered into the Julia general registry, and therefore can be deployed with the standard Pkg
Julia package manager.
julia> using Pkg
julia> Pkg.add("Geant4")
Examples are located in a separate repository G4Examples to minimize dependencies since they use additional functionality such as graphics, plotting, analysis tools, etc. To use and play with the examples, the user can clone the examples repository and setup a complete Julia environment with:
$ git clone https://github.com/JuliaHEP/G4Examples.jl.git
$ julia --project=G4Examples.jl -e 'import Pkg; Pkg.instantiate()'
Getting started
Import the Geant4
module. All the wrapped Geant4 classes are exported since they are prefixed by G4
minimizing the chances of a name clash with other Julia symbols.
julia> using Geant4
julia> runManager = G4RunManager()
**************************************************************
Geant4 version Name: geant4-11-01-patch-01 [MT] (10-February-2023)
Copyright : Geant4 Collaboration
References : NIM A 506 (2003), 250-303
: IEEE-TNS 53 (2006), 270-278
: NIM A 835 (2016), 186-225
WWW : http://geant4.org/
**************************************************************
Geant4.G4RunManagerAllocated(Ptr{Nothing} @0x00007f9fcb6f9c50)
julia> methodswith(G4RunManager, supertypes=true)
[1] convert(t::Type{G4RunManager}, x::T) where T<:G4RunManager in Geant4 at /Users/mato/.julia/packages/CxxWrap/IdOJa/src/CxxWrap.jl:676
[2] AbortEvent(arg1::Union{CxxWrap.CxxWrapCore.CxxRef{<:G4RunManager}, Union{CxxWrap.CxxWrapCore.SmartPointer{T2}, T2} where T2<:G4RunManager}) in Geant4 at /Users/mato/.julia/packages/CxxWrap/IdOJa/src/CxxWrap.jl:618
...
[94] rndmSaveThisRun(arg1::Union{CxxWrap.CxxWrapCore.CxxRef{<:G4RunManager}, Union{CxxWrap.CxxWrapCore.SmartPointer{T2}, T2} where T2<:G4RunManager}) in Geant4 at /Users/mato/.julia/packages/CxxWrap/IdOJa/src/CxxWrap.jl:618
julia> v = GetVersionString(runManager)
ConstCxxRef{G4String}(Ptr{G4String} @0x00007ffed34df2d8)
julia> String(v)
" Geant4 version Name: geant4-11-01-patch-01 [MT] (10-February-2023)"
Note that class methods are called with the object instance as first argument. In C++ the GetVersionString
method would be called as runManager->GetVersionString()
while in Julia it is called as GetVersionString(runManager)
. Thanks to the Julia multi-dispatch we do not need to prefix the methods with the module name Geant4.GetVersionString(runManager)
, even for very common function names such as mag
.
julia> v = G4ThreeVector(1,2,3)
Geant4.CLHEP!Hep3VectorAllocated(Ptr{Nothing} @0x00007f9fcaf2a710)
julia> mag(v)
3.7416573867739413
Geant4 Julia interface
The main goal for defining a Geant4 application in the Julia interface is to create an instance of the G4JLApplication
type, where all the needed elements for running a Geant4 application are declared, such as the detector geometry, the physics list, the primary particle generator, the type of run manager, the user actions, etc.
These are the needed elements:
- detector. An instance of a detector structure inheriting from the abstract type
G4JLDetector
, in which all the detector parameters are defined. The user should also provide a method specialization ofGeant4.getConstructor(::G4JLDetector)::Function
to return the Julia function that toolkit needs to call in order to construct the geometry and return the pointer of the 'world' physical volume. There is no default. - field. An instance of the magnetic field class. The
G4JLUniformMagField(...)
function provides a uniform magnetic field. See later how to define a custom one. - simdata. An instance of the simulation data structure that the program will need to collect during the simulation execution. This mutable structure needs to inherit from the abstract type
G4JLSimulationData
and is completely user defined with counters, data structures to collect the hits or doses, histograms, etc. The default is an instance of typeG4JLNoData
. - nthreads. Number of worker threads to be used. The default is 0, which means serial mode. Any number > 0 will use the MT functionality of Geant4, and therefore the user would need to pay attention to the user actions that are run concurrently to avoid data races (see Julia doc on multi-threading)
- verbose. Verbosity level (for physics list). The default is 0.
- physics_type. The physics list predefined type. Default is
FTFP_BERT
. - generator_type. The primary generator generator type. The default is
G4JLParticleGun
, which encapsulates aG4ParticleGun
. The underlyingG4ParticleGun
can be obtained by callingGetGun()
. - user actions. Julia methods defining the different possible user actions (e.g. stepping action, tracking action, run action, event action). The default is no action.
- sdetectors. List of sensitive detectors. This is given as a
Vector
of pairslv::String => sd::G4JLSensitiveDetector
to associate logical volumes by name to sensitive detector instances (see next section). - scorers. List of scoring meshes defined with the function
G4JLScoringMesh
.
Once the G4JLApplication
is instantiated (and implicitly an instance of the G4RunManager
created), the user can control the application with the following commands:
configure(::G4JLApplication)
. It associates the physics list, generator and user actions to the selectedG4RunManager
instance.initialize(::G4JLApplication)
. It basically calls theInitialize()
method of the run manager and associate the declared sensitive detectors.reinitialize(::G4JLApplication, ::G4JLDetector)
. It re-defines the declared detector geometry with a new instance ofG4JLDetector
.beamOn(::G4JLApplication, ::Int)
. Starts a run with a given number of events.
Constructing the detector
Parameters of the detector are collected in a user defined mutable data structure inheriting from G4JLDetector
. The user also needs to provide a Julia method for constructing the geometry. This method needs to have the signature
<User_Det_Constructor_Function>(::G4JLDetector)::CxxPtr{G4VPhysicalVolume}
The only argument of the function gives access to the user defined structure with all the detector parameters.
The type CxxPtr{G4VPhysicalVolume}
denotes a C++ pointer to the G4VPhysicalVolume
type.
The user can use the native G4 classes for constructing the geometry such as the different type of solids (e.g. G4Box, G4Tubs, etc.), G4LogicalVolume
, G4PVPlacement
, G4PVReplica
, etc. Alternatively can the type G4JLDetectorGDML
to construct a detector from a GDML file.
Note that by default constructed C++ objects from Julia would get automatically deleted by the Julia garbage collector (GC) since a finalizer
gets installed to the wrapper classes. This is particularly a problem when constructing the geometry.
Currently for the following classes have the finalizer
disabled in the wrapper: G4PVPlacement
, G4LogicalVolume
, G4PVReplica
, G4Material
, G4Isotope
, G4Element
. This means that instances of them will not be deleted by Julia to avoid double deletion (often a crash) when the geometry gets deleted at the finalization of the application from the C++ side.
A pointer to any of the G4Solid
needs to be passed to the G4LogicalVolume
using move!(objref)
function to transfer the ownership of the referenced object to the C++ side. See the following example:
trackerS = G4Tubs("tracker", 0, trackerSize, trackerSize, 0, 360deg)
trackerLV = G4LogicalVolume(move!(trackerS), m_air, "Tracker")
G4PVPlacement(nothing, # no rotation
G4ThreeVector(0,0,0), # at (0,0,0)
trackerLV, # its logical volume
"Tracker", # its name
worldLV, # its mother volume
false, # no boolean operations
0, # copy number
checkOverlaps) # checking overlaps
Physics List
The user can provide one of the pre-defined physics lists, such as QGS_BIC
, QBBC
or FTFP_BERT
. Alternatively, the user can define a Julia struct
as a subtype of
G4VUserPhysicsList` and modify some of the physcis in the constructor. For example:
struct ScintPhysicsList <: G4VUserPhysicsList
function ScintPhysicsList(verbose)
pl = FTFP_BERT(verbose)
# replace G4EmStandardPhysics
ReplacePhysics(pl, move!(G4EmStandardPhysics_option4(verbose)))
# register G4OpticalPhysics
RegisterPhysics(pl, move!(G4OpticalPhysics(verbose)))
# activate scintillation
optpar = G4OpticalParameters!Instance()
SetProcessActivation(optpar, "Scintillation", true)
# activate cherenkov radiation
SetProcessActivation(optpar, "Cerenkov", true)
return pl
end
end
Magnetic field
The user can define either an uniform magnetic field or a custom one. To define an custom one:
- define first a user structure for the parameters that inherits from the abstract type
G4JLFieldData
- then, define a function with the signature
(result::G4ThreeVector, position::G4ThreeVector, params::G4JLFieldData)::Nothing
- and finally, with all this, instantiate the magnetic field calling the function
G4JLMagneticField(<name>, <data>; getfield_method=<function>)
Primary Particle Generator
The user can define either a custom primary particle generator or use one of the two defined ones:
G4JLGunGenerator. Uses the
G4ParticleGun
class of Geant4 that generates a single particle type with a fix kinetic energy, position and direction.G4JLGunGenerator(particle = "proton", energy = 3GeV, direction = G4ThreeVector(0,0,1), position = G4ThreeVector(0,0,-2940.0))
G4JLGeneralParticleSource. Uses the
G4GeneralParticleSource
class of Geant4 to generate one of more sources of primary particles with predefined distributions for energy, position and direction. An example can be:src1 = (particle="e+", intensity=0.5, ene=(type="Exp", min=2MeV, max=10MeV, ezero=2.), pos=(type="Plane", shape="Circle", centre=G4ThreeVector(1cm,2cm,0cm), radius=3cm), ang=(type="cos", mintheta=0deg, maxtheta=180deg)) src2 = (particle="gamma", intensity=0.5, ene=(type="Brem", min=2MeV, max=10MeV, temp=2e12), pos=(type="Plane", shape="Ellipse", centre=G4ThreeVector(3cm,1cm,0cm), halfx=1cm, halfy=2cm), ang=(type="iso", mintheta=0deg, maxtheta=180deg)) gps = G4JLGeneralParticleSource(sources = [ src1, src2 ])
the standard particle gun parameters works as well:
G4JLGeneralParticleSource(particle = "proton", energy = 3GeV, direction = G4ThreeVector(0,0,1), position = G4ThreeVector(0,0,0))
Custom Generator. It is fairly simple to write a custom generator. It is needed to define a structure for the parameters to configure the generator, and a two functions to initialize and generate the primary particles called for each event. Here is an example:
# define the data structure with the generator parameters mutable struct PlaneSourceData <: G4JLGeneratorData particleName::String particlePtr::CxxPtr{G4ParticleDefinition} energy::Float64 halfx::Float64 halfy::Float64 position::G4ThreeVector direction::G4ThreeVector end # define the constructor with the default parameters function PlaneSource(;particle="gamma", energy=0.07MeV, halfx=7cm, halfy=7cm, position=G4ThreeVector(0,0,-14.9cm), direction=G4ThreeVector(0,0,1)) data = PlaneSourceData(particle, CxxPtr{G4ParticleDefinition}(C_NULL), energy, halfx, halfy, position, direction) function init(data:: PlaneSourceData, app::G4JLApplication) data.particlePtr = FindParticle(data.particleName) end function generate( evt::G4Event, data:: PlaneSourceData)::Nothing mass = data.particlePtr |> GetPDGMass momentum = √((mass + data.energy)^2 - mass^2) pvec = momentum * data.direction pos = data.position + G4ThreeVector( data.halfx * (rand() - 0.5), data.halfy * (rand() - 0.5), 0) primary = G4PrimaryParticle(data.particlePtr, pvec |> x, pvec |> y, pvec |> z ) vertex = G4PrimaryVertex(pos, 0ns) SetPrimary(vertex, move!(primary)) # note that we give up ownership of the objects just created AddPrimaryVertex(evt, move!(vertex)) # note that we give up ownership of the objects just created end G4JLPrimaryGenerator("PlaneSource", data; init_method=init, generate_method=generate) end
User Actions
User actions are native Julia functions that are callbacks of the Geant4 toolkit. They are declared in the constructor of G4JLApplication
, so they do not need to be associated to a specific function name. All user actions receive a reference to the G4JLApplication
from which the user can obtain details of the actual application, such as the current detector, the physics, the generator, or the running simulation data. There are the available attributes of the application instance:
runmanager::G4RunManager # The C++ G4RunManager instance
detector::DET # User defined detector structure with all parameters
simdata::Vector{DAT} # User defined simulation data structs (each worker has its own)
physics::G4VUserPhysicsList # Physics List
generator::G4JLPrimaryGenerator # Primary particle generator
field::Union{G4Field, G4JLMagneticField} # Magnetic field instance
evtdisplay::G4JLDisplay # Event display instance
nthreads::Int32 # number of worker threads
verbose::Int32 # verbosity level for physics lists
sdetectors::Dict{String,Vector{G4JLSensitiveDetector}} # dictionary of sensitive detectors
scorers::Vector{G4JLScoringMesh} # vector of scoring meshes
The following are the currently defined possible user actions:
- stepping action. Called on each simulation step. The signature is
(::G4Step, ::G4JLApplication)::Nothing
. Consult the G4Step reference manual to see what you can get from it. - pre-tracking action. Called at the creation of a new participle being tracked. The signature is
(::G4Track, ::G4JLApplication)::Nothing
. Consult the G4Step reference manual to see what you can get from it. - post-tracking action. Called at the end of the particle being tracked. The signature is
(::G4Track, ::G4JLApplication)::Nothing
. Consult the G4Track reference manual to see what you can get from it. - begin-event action. Called at the beginning of each event. The signature is
(::G4Event, ::G4JLApplication)::Nothing
. Consult the G4Event reference manual to see what you can get from it. - end-event action. Called at the end of each event. The signature is
(::G4Event, ::G4JLApplication)::Nothing
. Consult the G4Event reference manual to see what you can get from it. - begin-run action. Called at the beginning of a run. The signature is
(::G4Run, ::G4JLApplication)::Nothing
. Consult the G4Run reference manual to see what you can get from it. - end-run action. Called at the end of a run. The signature is
(::G4Run, ::G4JLApplication)::Nothing
. Consult the G4Run reference manual to see what you can get from it.
Sensitive Detectors
The user can define sensitive detectors by defining a data structure and 3 callback functions, which will initialize, fill and dispose the defined data structure. Later, the instantiated sensitive detector would be associated to one or more logical volumes of the detector setup. Instantiating a G4JLSensitiveDetector
will require the following arguments:
- name. A string to identify the sensitive detector. No default.
- sd data. A instance of a user defined
G4JLSDData
mutable data structure that will passed to each callback invocation. - initialize method. User method that is called at the beginning of each event. The signature is
(::B2aSDData, ::G4HCofThisEvent)::Nothing
. - endOfEvent method. User method that is called at the end of each event. The signature is
(::B2aSDData, ::G4HCofThisEvent)::Nothing
. - processHits method. User method that is called at simulation step that ends at the associated logical volume. The signature is
(::B2aSDData, ::G4Step, ::G4TouchableHistory)::Bool
. Consult the G4Step reference manual to see what you can get from the G4Step. It returns true if a true hit is generated.
The following is a example on how to define a sensitive detector
#--------------------------------------------------------------------------------------------------
#---Define Crystal Sensitive Detector--------------------------------------------------------------
#--------------------------------------------------------------------------------------------------
#---SD collected data------------------------------------------------------------------------------
struct CrystalSDData <: G4JLSDData
hitcollection::Vector{Hit} # define a hit collection
CrystalSDData() = new(Hit[])
end
#---Initialize method------------------------------------------------------------------------------
function crystal_initialize(::G4HCofThisEvent, data::CrystalSDData)::Nothing
empty!(data.hitcollection) # empty the hit collection at every event
return
end
#---Process Hit method-----------------------------------------------------------------------------
function crystal_processHits(step::G4Step, ::G4TouchableHistory, data::CrystalSDData)::Bool
part = step |> GetTrack |> GetParticleDefinition
part == G4OpticalPhoton && return false
edep = step |> GetTotalEnergyDeposit
edep < 0. && return false
pos = step |> GetPostStepPoint |> GetPosition
push!(data.hitcollection, Hit(0., pos, edep, ScintCryst)) # fill the collection with a new Hit
return true
end
#---Create SD instance-----------------------------------------------------------------------------
G4JLSensitiveDetector("Crystal_SD", CrystalSDData(); # name an associated data are mandatory
processhits_method=crystal_processHits, # process hist method (also mandatory)
initialize_method=crystal_initialize) # intialize method
Scoring meshes
The user can also specify scoring meshes to obtain quantities on the defined grid. In Geant4 this is achieved using a set of UI commands. In this Julia interface this functionality has been encapsulated in a number of data structures. The function to create a scoring mesh is G4JLScoringMesh
and receive as arguments the the type and dimensions of the mesh, the position, the rotation, the number of bins in each dimension, and the quantities to accumulate with eventually some filter conditions. See for example the scoring mesh from RE03:
sc1 = G4JLScoringMesh("boxMesh_1",
BoxMesh(1m,1m,1m),
bins = (30, 30, 30),
quantities = [ energyDeposit("eDep")
nOfStep("nOfStepGamma", filters=[ParticleFilter("gammafilter", "gamma")])
nOfStep("nOfStepEMinus", filters=[ParticleFilter("eMinusFilter", "e-")])
nOfStep("nOfStepEPlus", filters=[ParticleFilter("ePlusFilter", "e+")])
]
)
The scoring mesh is added into the 'scorersargument when constructing a [
G4JLApplication`](@ref). After a run hs been executed, the user can obtain the quantity values (sum, sum2,entries) on the 3D grid just calling by accessing the quantity as an attribute of the scoring mesh. The returned 3D Julia array is shaped to the declared bins.
julia> beamOn(app,10000)
julia> sum, sum2, entries = sc1.eDep
julia> typeof(sum)
Array{Float64, 3}
julia> typeof(entries)
Array{Int64, 3}
julia> size(sum)
(30, 30, 30)
Detector and Event Display
For visualization applications, the user can create an instance of G4JLEventDisplay([settings file])
and give it to the constructor of G4JLApplication
in the evtdisplay
attribute. The constructor accepts a visualization settings file that will overwrite the default settings in the file Geant4.jl/ext/G4Vis/settings.jl. The format of the settings is Julia NamedTuple
. Here is an example:
(
display = (
show_axis = false,
),
trajectories = (
color = :yellow,
),
)
Examples
basic/B1
This is most basic example. For this example we have kept the interface closer to the native C++ interface. The only difference with respect the C++, is that that we need to instantiate a G4JLDetectorConstruction
with the Julia function that will be callback to construct the detector geometry as argument. This is because we cannot inherit from Julia the C++ class G4VUserDetectorConstruction
, which is the way foreseen in the native Geant4 toolkit to provide the specific user detector geometry. Similarly, for the user actions and primary particle generator we need to instantiate a G4JLActionInitialization
. The interaction with the application is done with the G4UImanager
.
To run it execute
julia --project=examples examples/basic/B1/B1.j
or execute the notebook examples/basic/B1/B1.ipynb
basic/B2a
This example fills a vector of TrackerHit
that is stored in the B2aSDData
simulation data structure for each event. This is achieved with a sensitive detector associated to the Chamber
logical volume. The example is using the high-level Julia interface with the instantiation of a G4JLApplication
declaring all the elements of the application (geometry, physics, simulation data, user actions, etc.)
To run it execute
julia --project=examples examples/basic/B2/B2a.jl
extended/RE03
This example makes use of the built-in scoring capability of Geant4 with a new Julia API interface creating an instance of G4JLScoringMesh
, instead of using the native Geant4 UI. The user defines a scoring mesh, and quantities to be collect and gets the results after the run.
To run it execute
julia --project=examples examples/extended/RE03/RE03.jl
WaterPhantom
This example is similar to RE03 but it defines a custom primary particle generator (MedicalBeam
) instead of using the predefined particle gun generator (G4JLGunGenerator
). It is a notebook and produces plots after each run. To run it execute
jupyter notebook examples/WaterPhantom/WaterPhantom.ipynb
See the rendered notebook.
TestEm3
This example comes from extended/electromagnetic/TestEm3 example. Since it requires additional packages such as FHist
for histograms and Plots
for their visualization, it has its own Julia environment in the folder examples/TestEm3
. It uses the Julia high-level interface with the instantiation of a G4JLApplication
declaring all the elements of the application.
To run it, execute
julia --project=examples -i examples/TestEm3/TestEm3.jl
Scintillator
Example with optical photons and customized physics list, as well as with a couple of sensitive detectors (for the crystal and silicon) and some simple analysis of the results. To run it, execute
julia --project=examples -i examples/Scintillator/Scintillator.jl
Visualization examples
The Geant4.jl project includes additional functionality for visualization under the extension directory ext/G4Vis/examples
. This is done in a different directory to separate and minimise the dependencies. The julia ext/G4Vis/examples/Project.toml
file has the complete list of dependencies needed for running these examples. In order to load all the required dependencies the user can execute the first time:
`julia --project=Geant4.jl/ext/G4Vis/examples -e 'import Pkg; Pkg.instantiate()'`
Note that depending on the actual platform and the desired interactivity, the user may need to choose a different Makie.jl
backend among the existing ones (GLMakie
, CairoMakie
, WGLMakie
,...).
B1vis.jl
This example uses the GLMakie
backend (OpenGL) of Makie. The use may change to other backends depending on his/her setup. To visualize the B1 detector do:
julia --project=ext/G4Vis/examples -i ext/G4Vis/examples/B1vis.jl
Note the option -i
to keep the interactive session open
B2aVis.jl
This example to visualize the detector and (a very simplistic) visualization of one event.
julia --project=ext/G4Vis/examples -i ext/G4Vis/examples/B2avis.jl
Solids.ipynb
This notebook shows all the possible solids in Geant4. This is work in progress and some solids do not have graphical representation yet.
jupyter notebook ext/G4Vis/examples/Solids.ipynb
See the rendered notebook
HBC30
This example mimics the CERN 30cm liquid hydrogen bubble chamber. It demonstrates the use of a uniform magnetic field (G4JLUniformMagField
). It is useful for displaying the detector and the produced particles in a customizable manner.
To run it, execute
julia --project=ext/G4Vis/examples -i ext/G4Vis/examples/HBC30/HBC30.jl
It also exists in a notebook format.
Building the wrapper code
We use the Geant4 native binary libraries and data from the binary package Geant4_jll, which has been produced with the BinaryBuilder
recipe. The wrapper library is downloaded from the binary package Geant4_julia_jll.
We have the possibility during the development of this package to re-generate locally new C++ wrapper code. For this we need to have wrapit
installed, which itself requires libclang
to be installed. If the executable is not found (not in the PATH), we can use the wrapper code that is already pre-generated and distributed with this package.
- what C++ classes get wrapped is controlled by the file
gen/Geant4.wit.in
. See the documentation of WrapIt for the details. - run the
gen/build.jl
script generate the wrapper code (if wrapit found) and build the wrapper library.
Please note that compiling the single generated wrapper file takes very long. This is due to the current implementation of WrapIt that places all wrapped types in a single file. This may change in the future.
Once the wrapper code is stabilized we move the generated code to the repository Geant4_cxxwrap to regenerate the binary package Geant4_julia_jll
using the BinaryBuilder
.