Skip to content

Commit 752f296

Browse files
perematoperemato
committed
Added B4a example
1 parent 6186989 commit 752f296

File tree

4 files changed

+315
-2
lines changed

4 files changed

+315
-2
lines changed

‎docs/make.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ function create_extras(extrafiles...)
3636
return extra_mds
3737
end
3838

39-
basic_mds = process_literate("B1", "B2a", "B2aVis", "B3a")
39+
basic_mds = process_literate("B1", "B2a", "B2aVis", "B3a", "B4a")
4040
extend_mds = process_literate("GPS", "RE03", "TestEm3", "Solids")
4141
advanced_mds = process_literate("TPCSim", "HBC30", "WaterPhantom", "UserLib", "JuliaAction")
4242
extra_mds = create_extras("B2aDetector.jl", "B2aVisSettings.jl", "B3Detector.jl", "GPSDetector.jl",

‎docs/src/examples/B4a.lit‎

Lines changed: 301 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,301 @@
1+
# # Basic/B4a Example
2+
#
3+
# This example simulates a simple Sampling Calorimeter setup.
4+
# To demonstrate several possible ways of data scoring, the example is provided in four variants: %B4a, %B4b, %B4c, %B4d.
5+
# See [README](https://github.com/Geant4/geant4/blob/master/examples/basic/B5/README.md)
6+
# file for the example.
7+
8+
#md # !!! note "Note that"
9+
#md # You can also download this example as a
10+
#md # [Jupyter notebook](B4a.ipynb) and a plain
11+
#md # [Julia source file](B4a.jl).
12+
#
13+
#md # #### Table of contents
14+
#md # ```@contents
15+
#md # Pages = ["B4a.md"]
16+
#md # Depth = 2:3
17+
#md # ```
18+
19+
# ## Loading the necessary Julia modules
20+
# Load the `Geant4` and `Geant4.SystemOfUnits` modules.
21+
using Geant4
22+
using Geant4.SystemOfUnits: cm, cm3, mm, pGy, eplus, keV, g, eV, MeV, mole
23+
using FHist, Plots
24+
25+
## to force loading G4Vis extension we need to load the following module
26+
#md using CairoMakie
27+
#nb using CairoMakie
28+
#jl using GLMakie
29+
#md import DisplayAs: PNG #hide
30+
31+
# ## Define Detector Parameters struct
32+
# The `B4Detector` structure is defined with the default detector parameters.
33+
34+
mutable struct B4aDetector <: G4JLDetector
35+
## main input parameters
36+
const nofLayers::Int
37+
const checkOverlaps::Bool
38+
const absoThickness::Float64
39+
const gapThickness::Float64
40+
const calorSizeXY::Float64
41+
42+
## constructor with defaults values for parameters
43+
function B4aDetector(; nofLayers::Int=10,
44+
checkOverlaps::Bool=true,
45+
absoThickness::Float64=10mm,
46+
gapThickness::Float64=5mm,
47+
calorSizeXY::Float64=10cm)
48+
self = new(nofLayers, checkOverlaps, absoThickness, gapThickness, calorSizeXY)
49+
return self
50+
end
51+
end
52+
_layerThickness(d::B4aDetector) = d.absoThickness + d.gapThickness
53+
_calorThickness(d::B4aDetector) = d.nofLayers * _layerThickness(d)
54+
_worldSizeZ(d::B4aDetector) = 1.2 * _calorThickness(d)
55+
_worldSizeXY(d::B4aDetector) = 1.2 * d.calorSizeXY
56+
57+
# ## Defining the geometry constructor
58+
function B4aConstruct(det::B4aDetector)::CxxPtr{G4VPhysicalVolume}
59+
(; nofLayers, checkOverlaps, absoThickness, gapThickness, calorSizeXY) = det
60+
nist = G4NistManager!Instance()
61+
62+
63+
layerThickness = _layerThickness(det)
64+
calorThickness = _calorThickness(det)
65+
worldSizeXY = _worldSizeXY(det)
66+
worldSizeZ = _worldSizeZ(det)
67+
68+
G4Material("liquidArgon", z = 18., a = 30.95g/mole, density = 1.390g/cm3)
69+
70+
defaultMaterial = FindOrBuildMaterial(nist, "G4_Galactic")
71+
absorberMaterial = FindOrBuildMaterial(nist, "G4_Pb")
72+
gapMaterial = FindOrBuildMaterial(nist, "liquidArgon")
73+
74+
75+
76+
worldS = G4Box("World", 0.5 * worldSizeXY, 0.5 * worldSizeXY, 0.5 * worldSizeZ)
77+
78+
worldLV = G4LogicalVolume(worldS, defaultMaterial, "WorldLV")
79+
80+
worldPV = G4PVPlacement(nothing, # no rotation
81+
G4ThreeVector(), # at (0,0,0)
82+
worldLV, # its logical volume
83+
"World", # its name
84+
nothing, # its mother volume
85+
false, # no boolean operation
86+
0, # copy number
87+
checkOverlaps) # overlaps checking
88+
89+
##---------------------------------------
90+
calorimeterS = G4Box("Calorimeter", 0.5 * calorSizeXY, 0.5 * calorSizeXY, 0.5 * calorThickness)
91+
92+
93+
calorLV = G4LogicalVolume(calorimeterS, defaultMaterial, "CalorLV")
94+
95+
println(calorLV)
96+
97+
G4PVPlacement(nothing,
98+
G4ThreeVector(),
99+
calorLV,
100+
"CalorPV",
101+
worldLV,
102+
false,
103+
0,
104+
checkOverlaps)
105+
##----------------------------------------------
106+
107+
layerS = G4Box("Layer",
108+
0.5*calorSizeXY, 0.5*calorSizeXY, 0.5*layerThickness)
109+
110+
layerLV = G4LogicalVolume(layerS,
111+
defaultMaterial,
112+
"Layer")
113+
114+
G4PVReplica("Layer",
115+
layerLV,
116+
calorLV,
117+
kZAxis,
118+
nofLayers,
119+
layerThickness)
120+
121+
##---------------------------------------------------
122+
absorberS = G4Box("Abso",
123+
0.5*calorSizeXY, 0.5*calorSizeXY, 0.5*absoThickness)
124+
125+
absorberLV = G4LogicalVolume(absorberS, absorberMaterial, "AbsoLV")
126+
127+
G4PVPlacement(nothing,
128+
G4ThreeVector(0., 0., -0.5*gapThickness),
129+
absorberLV,
130+
"AbsPV",
131+
layerLV,
132+
false,
133+
0,
134+
checkOverlaps)
135+
##--------------------------------------
136+
137+
gapS = G4Box("Gap", 0.5*calorSizeXY, 0.5*calorSizeXY, 0.5*gapThickness)
138+
139+
gapLV = G4LogicalVolume(gapS, gapMaterial, "Gap")
140+
141+
G4PVPlacement(nothing,
142+
G4ThreeVector(0., 0., 0.5*absoThickness),
143+
gapLV,
144+
"GapPV",
145+
layerLV,
146+
false,
147+
0,
148+
checkOverlaps)
149+
150+
151+
##------------------------------------------
152+
##Visualization attributes
153+
SetVisAttributes(worldLV, G4VisAttributes!GetInvisible())
154+
SetVisAttributes(calorLV, G4VisAttributes!GetInvisible())
155+
156+
return worldPV # return a pointer to the G4PhysicalVolume
157+
end
158+
159+
# Instantiate the detector
160+
detector = B4aDetector()
161+
162+
# ## Primary Particle Generator
163+
worldZHalfLength = _worldSizeZ(detector)/2
164+
particlegun = G4JLGunGenerator(particle = "e-",
165+
energy = 300MeV,
166+
direction = G4ThreeVector(0., 0., 1.),
167+
position = G4ThreeVector(0,0,-worldZHalfLength))
168+
169+
# ## Defining the simulation data
170+
const Hist1D64 = Hist1D{Float64}
171+
mutable struct B4aSimData <: G4JLSimulationData
172+
fEnergyDeposit_Abs::Float64
173+
fEnergyDeposit_Gap::Float64
174+
175+
fStepLength_Abs::Float64
176+
fStepLength_Gap::Float64
177+
178+
fEdepHist_Abs::Hist1D64
179+
fStepLenHist_Abs::Hist1D64
180+
181+
fEdepHist_Gap::Hist1D64
182+
fStepLenHist_Gap::Hist1D64
183+
184+
B4aSimData() = new()
185+
end
186+
187+
# ## User Actions
188+
# ## Begin Run Action
189+
function beginrun(run::G4Run, app::G4JLApplication)::Nothing
190+
data = getSIMdata(app)
191+
192+
data.fEdepHist_Abs = Hist1D(;binedges=100.:1.:300.)
193+
data.fStepLenHist_Abs = Hist1D(;binedges=0.:2.:500.)
194+
195+
data.fEdepHist_Gap = Hist1D(;binedges=0.:1.:200.)
196+
data.fStepLenHist_Gap = Hist1D(;binedges=0.:2.:500.)
197+
nothing
198+
end
199+
200+
# ### End Run Action
201+
function endrun(run::G4Run, app::G4JLApplication)::Nothing
202+
nothing
203+
end
204+
# ### Begin Event Action
205+
function beginevent(evt::G4Event, app::G4JLApplication)
206+
G4JL_println("===============started event $(evt |> GetEventID)")
207+
data = getSIMdata(app)
208+
data.fEnergyDeposit_Abs = 0.0
209+
data.fStepLength_Abs = 0.0
210+
211+
data.fEnergyDeposit_Gap = 0.0
212+
data.fStepLength_Gap = 0.0
213+
return
214+
end
215+
# ###End event action
216+
function endevent(evt::G4Event, app::G4JLApplication)
217+
data = getSIMdata(app)
218+
## G4JL_println("AbsEdep: $(data.fEnergyDeposit_Abs)")
219+
## println("AbsStepLength: ", data.fStepLength_Abs)
220+
221+
push!(data.fEdepHist_Abs, data.fEnergyDeposit_Abs)
222+
push!(data.fStepLenHist_Abs, data.fStepLength_Abs)
223+
224+
push!(data.fEdepHist_Gap, data.fEnergyDeposit_Gap)
225+
push!(data.fStepLenHist_Gap, data.fStepLength_Gap)
226+
227+
println("Total Energy deposited in Absorber: ",data.fEnergyDeposit_Abs)
228+
println("Total Energy deposited in Gaps: ", data.fEnergyDeposit_Gap)
229+
G4JL_println("================event ended $(evt |> GetEventID) \n")
230+
return
231+
end
232+
# ### Stepping Action
233+
function stepaction(step::G4Step, app::G4JLApplication)
234+
data = getSIMdata(app)
235+
236+
volume = step |> GetPreStepPoint |> GetPhysicalVolume |> GetName
237+
edep = step |> GetTotalEnergyDeposit
238+
239+
stepLength = 0
240+
if step |> GetTrack |> GetDefinition |> GetPDGCharge != 0.
241+
stepLength = step |> GetStepLength
242+
end
243+
244+
if volume[] == "AbsPV"
245+
data.fEnergyDeposit_Abs += edep
246+
data.fStepLength_Abs += stepLength
247+
end
248+
249+
if volume[] == "GapPV"
250+
data.fEnergyDeposit_Gap += edep
251+
data.fStepLength_Gap += stepLength
252+
end
253+
return
254+
end
255+
256+
# ## Making the Application
257+
evtdisplay = G4JLEventDisplay(joinpath(@__DIR__, "B4aVis.jl"))
258+
Geant4.getConstructor(::B4aDetector)::Function = B4aConstruct
259+
app = G4JLApplication(detector = detector, # detector with parameters
260+
simdata = B4aSimData(), # simulation data
261+
generator = particlegun, # primary particle generator
262+
nthreads = 0, # # of threads (0 = no MT)
263+
physics_type = FTFP_BERT, # what physics list to instantiate
264+
evtdisplay = evtdisplay, # event display
265+
stepaction_method = stepaction,
266+
beginrunaction_method = beginrun, # begin-run action (initialize counters and histograms)
267+
endrunaction_method = endrun,
268+
begineventaction_method=beginevent,
269+
endeventaction_method = endevent
270+
)
271+
272+
configure(app)
273+
initialize(app)
274+
275+
# ## Run the Example
276+
# Run a single event and display it
277+
## ui`/tracking/verbose 1`
278+
beamOn(app,1)
279+
#nb display(evtdisplay.figure)
280+
#md PNG(evtdisplay.figure) #hide
281+
#jl display(evtdisplay.figure)
282+
283+
# Run for 100 events
284+
beamOn(app, 100)
285+
286+
# ## Plotting the histograms
287+
data = getSIMdata(app)
288+
lay = @layout [°; °; °; °]
289+
img = Plots.plot(layout=lay, show=true, size=(800,1000))
290+
291+
Plots.plot!(subplot=1, data.fEdepHist_Abs, title="Total EDep in Absorber", xlabel="Edep (MeV)", label="Abs Edep", show=true)
292+
Plots.plot!(subplot=2, data.fEdepHist_Gap, title="Total EDep in Gaps", xlabel="Edep (MeV)", label="Gap_Edep", show=true)
293+
294+
Plots.plot!(subplot=3, data.fStepLenHist_Abs, title="Total StepLeng in Absorbers", xlabel="Step Length (mm)", label="Abs_StepLen", show=true)
295+
Plots.plot!(subplot=4, data.fStepLenHist_Gap, title="Total StepLeng in Gaps", xlabel="Step Length (mm)", label="Gap_StepLen", show=true)
296+
#nb PNG(img) #hide
297+
#md PNG(img) #hide
298+
#jl display("image/png", img)
299+
300+
301+

‎docs/src/examples/B4aVis.jl‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
(
2+
display = (
3+
backgroundcolor = :black, # Display background color
4+
show_axis = true,
5+
# resolution = (640, 360),
6+
camera_zoom = 0.1,
7+
),
8+
trajectories = (
9+
color = :yellow,
10+
),
11+
)

‎docs/src/releases.md‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,8 @@
44
- New Features:
55
- Added new G4 classes to the Wrapper: G4VTwistedFaceted, G4Hype, G4TessellatedSolid, G4TriangularFacet, G4Polyhedron, HepPolyhedron
66
- Change the way to create the meshes of Solids for G4Vis using the GetPolyhedron() virtual function. Remove all obsolete code.
7-
- When drawing recursive geometry hierarchies, merging meshes is done before drawing them.
7+
- When drawing recursive geometry hierarchies, merging meshes is done before drawing them.
8+
- Added B4a example (thanks to @Chetan-379)
89

910
### 0.2.3 - 29-Aug-2025
1011
- New Features:

0 commit comments

Comments
 (0)