Scoring Meshes#
The user can 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.
using Geant4
using Geant4.SystemOfUnits
using CairoMakie
Let’s define a simple detector (from the previous notebook on sensitive detectors)
struct SimpleDetector <: G4JLDetector ; end
#---Materials----------------------------------------------------------------------------------
nist = G4NistManager!Instance()
m_air = FindOrBuildMaterial(nist, "G4_AIR")
m_bgo = FindOrBuildMaterial(nist, "G4_BGO")
#---Volumes------------------------------------------------------------------------------------
worldS = G4Box("world", 1m, 1m, 1m)
worldLV = G4LogicalVolume(worldS, m_air, "World")
worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, "World", nothing, false, 0, false)
crystalS = G4Box("world", 5cm, 5cm, 10cm)
crystalLV = G4LogicalVolume(crystalS, m_bgo, "Crystal")
crystalPV = G4PVPlacement(nothing, G4ThreeVector(), crystalLV, "Crystal", worldLV, false, 0, false)
#---define getConstructor
Geant4.getConstructor(::SimpleDetector)::Function = (::SimpleDetector) -> worldPV
We create a scoring mesh of the same size as the crystal and with bins 20, 20 and 40 in z direction to collect energy deposit, number of steps produced by gammas, electrons and positrons.
sc1 = G4JLScoringMesh("boxMesh_1",
BoxMesh(5cm,5cm,20cm),
bins = (20, 20, 40),
quantities = [ energyDeposit("eDep")
nOfStep("nOfStepGamma", filters=[ParticleFilter("gammafilter", "gamma")])
nOfStep("nOfStepEMinus", filters=[ParticleFilter("eMinusFilter", "e-")])
nOfStep("nOfStepEPlus", filters=[ParticleFilter("ePlusFilter", "e+")])
]
);
Let’s initialize now the application
particlegun = G4JLGunGenerator(particle = "e-",
energy = 3GeV,
direction = G4ThreeVector(0,0,1),
position = G4ThreeVector(0,0,-1m))
app = G4JLApplication(;detector = SimpleDetector(), # detector with parameters
generator = particlegun, # primary particle generator
nthreads = 4, # number of threads (MT)
physics_type = FTFP_BERT, # what physics list to instantiate
scorers = [sc1] # list of scorers
);
**************************************************************
Geant4 version Name: geant4-11-03-patch-02 [MT] (25-April-2025)
<< in Multi-threaded mode >>
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(app)
initialize(app)
--- G4CoupledTransportation is used
.... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4ScoringManager has 1 scoring meshes.
G4ScoringBox : boxMesh_1 --- Shape: Box mesh
Size (x, y, z): (5, 5, 20) [cm]
# of segments: (20, 20, 40)
displacement: (0, 0, 0) [cm]
registered primitve scorers :
0 eDep
1 nOfStepGamma with gammafilter
2 nOfStepEMinus with eMinusFilter
3 nOfStepEPlus with ePlusFilter
run for 1000 events
beamOn(app,1000)
G4WT1 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT2 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT3 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT0 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT1 > G4ScoringManager has 1 scoring meshes.
G4WT1 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT2 > G4ScoringManager has 1 scoring meshes.
G4WT1 > Size (x, y, z): (5, 5, 20) [cm]
G4WT2 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT1 > # of segments: (20, 20, 40)
G4WT2 > Size (x, y, z): (5, 5, 20) [cm]
G4WT1 > displacement: (0, 0, 0) [cm]
G4WT2 > # of segments: (20, 20, 40)
G4WT1 > registered primitve scorers :
G4WT1 > 0 eDep
G4WT1 > 1 nOfStepGamma with gammafilter
G4WT1 > 2 nOfStepEMinus with eMinusFilter
G4WT0 > G4ScoringManager has 1 scoring meshes.
G4WT2 > displacement: (0, 0, 0) [cm]
G4WT2 > registered primitve scorers :
G4WT2 > 0 eDep
G4WT2 > 1 nOfStepGamma with gammafilter
G4WT1 > 3 nOfStepEPlus with ePlusFilter
G4WT2 > 2 nOfStepEMinus with eMinusFilter
G4WT2 > 3 nOfStepEPlus with ePlusFilter
G4WT3 > G4ScoringManager has 1 scoring meshes.
G4WT3 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT0 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT3 > Size (x, y, z): (5, 5, 20) [cm]
G4WT3 > # of segments: (20, 20, 40)
G4WT3 > displacement: (0, 0, 0) [cm]
G4WT0 > Size (x, y, z): (5, 5, 20) [cm]
G4WT0 > # of segments: (20, 20, 40)
G4WT0 > displacement: (0, 0, 0) [cm]
G4WT3 > registered primitve scorers :
G4WT3 > 0 eDep
G4WT0 > registered primitve scorers :
G4WT3 > 1 nOfStepGamma with gammafilter
G4WT0 > 0 eDep
G4WT3 > 2 nOfStepEMinus with eMinusFilter
G4WT0 > 1 nOfStepGamma with gammafilter
G4WT3 > 3 nOfStepEPlus with ePlusFilter
G4WT0 > 2 nOfStepEMinus with eMinusFilter
G4WT0 > 3 nOfStepEPlus with ePlusFilter
Now the scoring mesh variable sc1 holds the collected information
sum, sum2, entries = sc1.eDep;
Each of these variables (sum, sum2 and entries) is a 20×20×40 Array corresponding to the deposit energy.
Let’s see a X-Y plane of it in the middle of the Z axis:
sum[:,:,20]
20×20 Matrix{Float64}:
13.7256 14.9938 7.92104 13.8171 … 7.27156 9.63841 19.785
12.5688 29.1449 14.7244 27.9177 13.7953 13.2012 15.1446
9.0862 26.3659 18.3457 35.1564 17.8024 13.9003 15.4984
12.3851 15.7975 41.7305 41.863 26.1884 19.0315 19.7014
14.3283 31.9743 38.3296 34.5613 34.2372 26.2778 20.1119
26.0273 36.8281 46.8928 63.5484 … 38.885 36.7931 43.1647
31.7248 24.7029 34.6437 66.118 38.5972 26.8317 24.2865
30.0158 42.4134 61.1366 115.708 71.6774 37.5478 34.0374
31.9256 30.7449 59.6449 132.869 76.5851 65.04 36.357
38.0952 28.9673 71.953 125.173 48.5948 55.2591 23.0177
31.6647 45.7067 90.3182 95.1163 … 49.6353 53.9987 23.6632
40.1828 45.9377 59.8547 69.8218 60.2548 38.2881 29.3659
32.4741 28.0979 54.3396 94.8621 52.9933 61.4301 31.4744
23.9163 44.4536 51.3471 90.3794 43.5916 42.0923 21.3182
13.7058 26.7813 56.5094 65.2599 23.2429 14.2696 24.3575
24.269 34.1277 34.0417 48.8177 … 26.33 43.6835 23.2543
21.4099 13.3081 13.878 24.8339 25.716 19.2284 15.9123
13.6347 7.63564 15.3051 30.4052 18.9042 17.8803 17.1869
13.7882 9.50164 14.0989 18.6022 16.4974 6.60949 7.72048
8.92128 6.3244 15.1491 7.80917 7.27676 5.53977 6.25419
A better way to to plot a heatmap with Makie.
img = heatmap(sum[:,:,20], color=:thermal, title="Edep (XY)")
display("image/png", img)
Invalid attributes color and title for plot type Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}}.
The available plot attributes for Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}} are:
alpha depth_shift inspector_hover nan_color transparency
clip_planes fxaa inspector_label overdraw visible
colormap highclip interpolate space
colorrange inspectable lowclip ssao
colorscale inspector_clear model transformation
Stacktrace:
[1] validate_attribute_keys(plot::Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}})
@ Makie ~/.julia/packages/Makie/4JW9B/src/recipes.jl:931
[2] push!(scene::Scene, plot::Plot)
@ Makie ~/.julia/packages/Makie/4JW9B/src/scenes.jl:512
[3] plot!
@ ~/.julia/packages/Makie/4JW9B/src/interfaces.jl:212 [inlined]
[4] plot!(ax::Axis, plot::Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}})
@ Makie ~/.julia/packages/Makie/4JW9B/src/figureplotting.jl:431
[5] plot!(fa::Makie.FigureAxis, plot::Heatmap{Tuple{Makie.EndPoints{Float32}, Makie.EndPoints{Float32}, Matrix{Float32}}})
@ Makie ~/.julia/packages/Makie/4JW9B/src/figureplotting.jl:427
[6] _create_plot(F::Function, attributes::Dict{Symbol, Any}, args::Matrix{Float64})
@ Makie ~/.julia/packages/Makie/4JW9B/src/figureplotting.jl:330
[7] #heatmap#47
@ ~/.julia/packages/Makie/4JW9B/src/recipes.jl:517 [inlined]
[8] top-level scope
@ In[9]:1
img = heatmap(sum[:,10,:]', color=:thermal, title="Edep (XZ)")
display("image/png", img)