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-02-patch-01 [MT] (16-February-2024)
<< 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)
G4WT3 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT0 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT1 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT2 > .... G4ScoringMessenger::MeshBinCommand - G4ScoringBox
G4WT0 > G4ScoringManager has 1 scoring meshes.
G4WT0 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT0 > Size (x, y, z): (5, 5, 20) [cm]
G4WT0 > # of segments: (20, 20, 40)
G4WT0 > displacement: (0, 0, 0) [cm]
G4WT0 > registered primitve scorers :
G4WT0 > 0 eDep
G4WT1 > G4ScoringManager has 1 scoring meshes.
G4WT1 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT0 > 1 nOfStepGamma with gammafilter
G4WT2 > G4ScoringManager has 1 scoring meshes.
G4WT0 > 2 nOfStepEMinus with eMinusFilter
G4WT2 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT3 > G4ScoringManager has 1 scoring meshes.
G4WT1 > Size (x, y, z): (5, 5, 20) [cm]
G4WT3 > G4ScoringBox : boxMesh_1 --- Shape: Box mesh
G4WT0 > 3 nOfStepEPlus with ePlusFilter
G4WT2 > Size (x, y, z): (5, 5, 20) [cm]
G4WT2 > # of segments: (20, 20, 40)
G4WT1 > # of segments: (20, 20, 40)
G4WT2 > displacement: (0, 0, 0) [cm]
G4WT3 > Size (x, y, z): (5, 5, 20) [cm]
G4WT1 > displacement: (0, 0, 0) [cm]
G4WT1 > registered primitve scorers :
G4WT3 > # of segments: (20, 20, 40)
G4WT1 > 0 eDep
G4WT3 > displacement: (0, 0, 0) [cm]
G4WT2 > registered primitve scorers :
G4WT2 > 0 eDep
G4WT2 > 1 nOfStepGamma with gammafilter
G4WT2 > 2 nOfStepEMinus with eMinusFilter
G4WT1 > 1 nOfStepGamma with gammafilter
G4WT3 > registered primitve scorers :
G4WT1 > 2 nOfStepEMinus with eMinusFilter
G4WT2 > 3 nOfStepEPlus with ePlusFilter
G4WT1 > 3 nOfStepEPlus with ePlusFilter
G4WT3 > 0 eDep
G4WT3 > 1 nOfStepGamma with gammafilter
G4WT3 > 2 nOfStepEMinus with eMinusFilter
G4WT3 > 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}:
22.6907 8.23238 9.53181 14.7324 … 12.409 7.39246 17.7415
11.2568 10.365 21.2517 24.7346 17.5209 15.3881 1.71277
12.9987 15.2373 14.4912 48.1607 22.5194 19.1845 11.8564
9.42243 19.6698 23.8529 37.7654 19.5384 21.3335 13.5281
22.8952 24.6838 37.5823 52.9955 30.5014 33.8069 13.4116
25.7023 25.9152 33.594 74.3156 … 26.7333 24.8169 28.8356
32.9647 38.8203 49.2023 69.869 34.747 17.5358 23.1117
26.5994 42.2314 62.3937 68.9426 74.4483 26.0788 33.6273
11.4288 38.1839 59.5371 116.112 54.4134 35.2306 21.319
16.6416 46.9006 65.8704 109.926 40.7341 46.4459 27.2181
19.2279 44.3896 75.6179 104.706 … 73.4129 59.2247 32.4233
20.4118 47.9795 59.2302 126.023 58.434 62.0383 41.5008
24.6917 46.6389 60.4266 95.9413 80.301 29.0629 28.9642
17.0803 32.8747 49.331 87.1838 37.8035 38.9716 39.1364
12.9969 36.0215 43.2987 53.7823 57.9662 18.7645 20.4736
13.24 24.9719 27.4121 42.8602 … 29.3471 21.7373 11.0561
9.799 21.3981 30.4914 31.8881 35.2386 31.7529 11.3309
6.4313 15.8679 37.8778 28.2031 33.2141 9.74859 10.3116
5.79865 13.9472 13.9546 27.2462 14.7916 8.13246 11.2221
6.29478 14.4728 7.10103 16.5264 7.39579 14.2824 9.63563
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{MakieCore.EndPoints{Float32}, MakieCore.EndPoints{Float32}, Matrix{Float32}}}.
The available plot attributes for Heatmap{Tuple{MakieCore.EndPoints{Float32}, MakieCore.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
Generic attributes are:
clip_planes dim_conversions model transformation yautolimits
cycle label rasterize xautolimits zautolimits
Stacktrace:
[1] validate_attribute_keys(plot::Heatmap{Tuple{MakieCore.EndPoints{Float32}, MakieCore.EndPoints{Float32}, Matrix{Float32}}})
@ MakieCore ~/.julia/packages/MakieCore/zO6H4/src/recipes.jl:758
[2] push!(scene::Scene, plot::Plot)
@ Makie ~/.julia/packages/Makie/pFPBw/src/scenes.jl:489
[3] plot!
@ ~/.julia/packages/Makie/pFPBw/src/interfaces.jl:413 [inlined]
[4] plot!(ax::Axis, plot::Heatmap{Tuple{MakieCore.EndPoints{Float32}, MakieCore.EndPoints{Float32}, Matrix{Float32}}})
@ Makie ~/.julia/packages/Makie/pFPBw/src/figureplotting.jl:412
[5] plot!(fa::Makie.FigureAxis, plot::Heatmap{Tuple{MakieCore.EndPoints{Float32}, MakieCore.EndPoints{Float32}, Matrix{Float32}}})
@ Makie ~/.julia/packages/Makie/pFPBw/src/figureplotting.jl:407
[6] _create_plot(F::Function, attributes::Dict{Symbol, Any}, args::Matrix{Float64})
@ Makie ~/.julia/packages/Makie/pFPBw/src/figureplotting.jl:318
[7] #heatmap#37
@ ~/.julia/packages/MakieCore/zO6H4/src/recipes.jl:501 [inlined]
[8] top-level scope
@ In[9]:1
img = heatmap(sum[:,10,:]', color=:thermal, title="Edep (XZ)")
display("image/png", img)