Scintillation Example

This example is adapted from the original at settwi/g4-basic-scintillation to the Geant4.jl interface. Low-energy gammas are shot onto a CeBr₃ scintillating crystal that is read out by a silicon detector. The example counts the scintillation photons and the hits in the silicon, and records the energy deposited in the crystal.

Note that

You can also download this example as a Jupyter notebook and a plain Julia source file.

Table of contents

Loading the necessary Julia modules

  • Geant4, Geant4.PhysicalConstants and Geant4.SystemOfUnits for the Geant4 simulation
  • FHist, Printf and Plots to fill and plot the histograms
using Geant4
using Geant4.PhysicalConstants
using Geant4.SystemOfUnits
using FHist, Printf, Plots

# using GLMakie, Rotations, LinearAlgebra, IGLWrap_jll  # to force loading G4Vis extension

Define the Detector Geometry

The materials and their optical properties are defined in Parameters.jl, and the ScintDetector geometry in Detector.jl. We include these files not to clutter the example.

include(joinpath(@__DIR__, "Parameters.jl"))
include(joinpath(@__DIR__, "Detector.jl"))

crysdet = ScintDetector()
# world = scintConstruct(crystal)
# draw(world[])
Main.var"##142".ScintDetector(50.0, true, true, 20.0, 10.0, 2.0, 1.0, 1.0, 10.0, true, CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x00006ffee909af50)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x000000000000996f)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x000000000000998e)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x00006ffeeb000008)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x0000000000000007)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x0000000000010000)), CxxWrap.CxxWrapCore.CxxPtr{Geant4.G4Material}(Ptr{Geant4.G4Material}(0x0000000000140000)))

Define the Physics List

We start from FTFP_BERT, replace the electromagnetic physics with G4EmStandardPhysics_option4 and add G4OpticalPhysics. The scintillation and Cherenkov processes are then activated.

struct ScintPhysicsList <: G4VUserPhysicsList
    function ScintPhysicsList(verbose)
        pl = FTFP_BERT(verbose)
        ReplacePhysics(pl, move!(G4EmStandardPhysics_option4(verbose)))
        RegisterPhysics(pl, move!(G4OpticalPhysics(verbose)))
        # need to enable scintillation
        optpar = G4OpticalParameters!Instance()
        SetProcessActivation(optpar, "Scintillation", true);
        # I have found Cherenkov radiation to be error-prone
        SetProcessActivation(optpar, "Cerenkov", true);
        return pl
    end
end

Define the Particle Gun

The primary particles are 30 keV gammas shot along -z towards the crystal.

particlegun = G4JLGunGenerator(particle = "gamma",
                               energy = 30keV,
                               direction = G4ThreeVector(0, 0, -1),
                               position = G4ThreeVector(0, 0, 2cm))
Geant4.G4JLGunGenerator("ParticleGun", Geant4.G4JLParticleGunData(nothing, "gamma", G4ThreeVector(0.0,0.0,-1.0), G4ThreeVector(0.0,0.0,20.0), 0.03), Geant4.var"#init#19"(), Geant4.var"#gen#20"(), Geant4.G4JLGeneratorAction[])

Define the Simulation Data

The data structure ScintSimData is filled by the user actions during the simulation. It holds the number of scintillation photons in the current event and three histograms.

const Hist1D64 = Hist1D{Float64}
mutable struct ScintSimData <: G4JLSimulationData
    scintPhotonsPerEvent::Int64
    scintPhotonsHisto::Hist1D64
    siHitsHisto::Hist1D64
    crysEdepHisto::Hist1D64
    ScintSimData() = new()
end

The add! function adds up the histograms of two ScintSimData instances.

function add!(x::ScintSimData, y::ScintSimData)
    x.scintPhotonsHisto += y.scintPhotonsHisto
    x.siHitsHisto += y.siHitsHisto
    x.crysEdepHisto += y.crysEdepHisto
end
add! (generic function with 1 method)

The do_plot function plots the simulation data.

function do_plot(data::ScintSimData)
    (;scintPhotonsHisto, siHitsHisto, crysEdepHisto) = data
    lay = @layout [° °; °]
    plot(layout=lay, show=true, size=(1400,1000))
    plot!(subplot=1, scintPhotonsHisto, title="# scintillating photons/event", xlabel="# photons", show=true)
    plot!(subplot=2, siHitsHisto, title="# hits in Silicon/event", xlabel="# hits")
    plot!(subplot=3, crysEdepHisto, title="energy deposited in crystal", xlabel="keV")
end
do_plot (generic function with 1 method)

Define the User Actions and Sensitive Detectors

The sensitive detectors crystal_SD and silicon_SD and the user actions (beginrun, endrun, beginevent, endevent and stepping) are defined in UserActions.jl.

include(joinpath(@__DIR__, "UserActions.jl"))
findOpticalBoundary (generic function with 1 method)

Create the Application

app = G4JLApplication(; detector = crysdet,                           # detector with parameters
                        simdata = ScintSimData(),                     # simulation data structure
                        generator = particlegun,                      # primary particle generator
                        nthreads = 0,                                 # # of threads (0 = no MT)
                        physics_type = ScintPhysicsList,              # what physics list to instantiate
                        stepaction_method = stepping,                 # step action method
                        begineventaction_method = beginevent,         # begin-event action (initialize per-event data)
                        endeventaction_method = endevent,             # end-event action (fill histograms per event)
                        beginrunaction_method = beginrun,             # begin run action
                        endrunaction_method = endrun,                 # end run action
                        sdetectors = ["si_log" => silicon_SD,
                                      "cebr3_log" => crystal_SD]      # mapping of LVs to SDs (+ means multiple LVs with same name)
                      );

**************************************************************
 Geant4 version Name: geant4-11-04-patch-01 [MT]   (13-March-2026)
                       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/
**************************************************************

Configure, Initialize and Run

configure(app)
initialize(app)

# ui`/tracking/verbose 1`
beamOn(app,1000)
Checking overlaps for volume cebr3_phys:0 (G4Box) ... OK! 
G4OpticalSurface: data file /home/runner/.julia/artifacts/ce912e6d316dde2fa4ffe0e3ba28061c14783ea3/share/Geant4/data/RealSurface2.2/Rough_LUT.z successfully read in.
LUT DAVIS - data file: Rough_LUT.z read in! 
G4OpticalSurface: data file /home/runner/.julia/artifacts/ce912e6d316dde2fa4ffe0e3ba28061c14783ea3/share/Geant4/data/RealSurface2.2/Rough_LUTR.z successfully read in.
LUT DAVIS - reflectivity data file: Rough_LUTR.z read in! 
Checking overlaps for volume ptfe_phys:0 (G4SubtractionSolid) ... OK! 
Checking overlaps for volume al_phys:0 (G4SubtractionSolid) ... OK! 
Checking overlaps for volume si_phys:0 (G4Box) ... OK! 
### Birks coefficients used in run time
Run ended with a total of 1000 events

Plot the results

do_plot(app.simdata[1])

This page was generated using Literate.jl.