Simulation of the Aleph TPC
In this example we demonstrate how to connect the PYTHIA MC generator to a to a Geant4 application that implements a very simplistic model of the ALEPH TPC.
The TPC detector is a simple cylinder filled with a mixture of Argon and Methane gas inside an axial 1.5 tesla magnetic field. The primary particles are generated by PYTHIA8 with the the LEP collision settings.
You can also download this example as a Jupyter notebook and a plain Julia source file.
Table of contents
- Loading the necessary Julia modules
- Define the TPC Detector
- Magnetic Field initialization
- Primary Particle Generator
- Event Display
- Create the Application
- Run the Application
Loading the necessary Julia modules
Geant4andGeant4.SystemOfUnitsfor the Geant4 simulationParametersfor the parameter handling in the detector definitionPYTHIA8for the PYTHIA8 generator
using Geant4
using Geant4.SystemOfUnits
using Parameters
using PYTHIA8
# to force loading G4Vis extension we need to load the following modules
using CairoMakie
using PYTHIA8: px, py, pzDefine the TPC Detector
The TPC detector is a simple cylinder filled with a mixture of Ar and CH4 gas
using Geant4.SystemOfUnits: m, cm, mole,cm3The data structure AlephTPC is defined with the default detector parameters
@with_kw mutable struct AlephTPC <: G4JLDetector
totalL::Float64 = 4.4m # total length of the TPC
innerR::Float64 = 0.31m # inner radius of the TPC
outerR::Float64 = 1.8m # outer radius of the TPC
checkOverlaps::Bool = false # do check overlaps when creating the geometry
endMain.var"##277".AlephTPCThe function construct is defined to create the geometry of the TPC detector. It receives an instance of the AlephTPC structure and returns a pointer to the physical volume of the world.
function construct(det::AlephTPC)::CxxPtr{G4VPhysicalVolume}
(; totalL, innerR, outerR, checkOverlaps) = det
##---Materials----------------------------------------------------------------------------------
nist = G4NistManager!Instance()
m_Ar = FindOrBuildMaterial(nist, "G4_Ar")
m_CH₄ = FindOrBuildMaterial(nist, "G4_METHANE")
ArCH₄ = G4Material("ArCH₄", 3.491*mg/cm3, 2)
AddMaterial( ArCH₄, m_Ar, 0.91 )
AddMaterial( ArCH₄, m_CH₄,0.09 )
m_ArCH₄ = FindOrBuildMaterial(nist, "ArCH₄")
m_air = FindOrBuildMaterial(nist, "G4_AIR")
m_Al = FindOrBuildMaterial(nist, "G4_Al")
##---Volumes------------------------------------------------------------------------------------
worldZHalfLength = 1.1 * totalL/2
worldS = G4Box("world", worldZHalfLength, worldZHalfLength, worldZHalfLength)
worldLV = G4LogicalVolume(worldS, m_air, "World")
worldPV = G4PVPlacement(nothing, G4ThreeVector(), worldLV, "World", nothing, false, 0, checkOverlaps)
fcOuterS = G4Tubs("fcOuter", outerR, outerR + 1.8cm, totalL/2, 0, 2π)
fcOuterLV = G4LogicalVolume(fcOuterS, m_Al, "fcOuter")
fcInnerS = G4Tubs("fcInner", innerR - 0.8cm, innerR, totalL/2, 0, 2π)
fcInnerLV = G4LogicalVolume(fcInnerS, m_Al, "fcInner")
chamberS = G4Tubs("chamber", innerR, outerR, totalL/2, 0, 2π)
chamberLV = G4LogicalVolume(chamberS, m_ArCH₄, "chamber")
G4PVPlacement(nothing, G4ThreeVector(), fcOuterLV, "fcOuter", worldLV, false, 0, checkOverlaps)
G4PVPlacement(nothing, G4ThreeVector(), fcInnerLV, "fcInner", worldLV, false, 0, checkOverlaps)
G4PVPlacement(nothing, G4ThreeVector(), chamberLV, "chamber", worldLV, false, 0, checkOverlaps)
##SetUserLimits(chamberLV, G4UserLimits(5mm))
##---Visualization attributes------------------------------------------------------------------
boxVisAtt = G4VisAttributes(G4Colour(1.0, 1.0, 1.0, 0.0))
chamberVisAtt = G4VisAttributes(G4Colour(0.0, 1.0, 1.0, 0.2))
fcVisAtt = G4VisAttributes(G4Colour(0.0, 1.0, 0.0, 0.1))
SetVisAttributes(worldLV, boxVisAtt)
SetVisAttributes(chamberLV, chamberVisAtt)
SetVisAttributes(fcInnerLV, fcVisAtt)
SetVisAttributes(fcOuterLV, fcVisAtt)
##---Always return the physical world-------------------------------------------------------------
return worldPV
end
Geant4.getConstructor(::AlephTPC)::Function = constructFinally, the detector is created with the default parameters
tpc = AlephTPC() # create the detectorMain.var"##277".AlephTPC
totalL: Float64 4400.0
innerR: Float64 310.0
outerR: Float64 1800.0
checkOverlaps: Bool false
Magnetic Field initialization
bfield = G4JLUniformMagField(G4ThreeVector(0,0, 1.5tesla))Geant4.G4JLUniformMagField("UnifiormB", Geant4.G4JLUniformMagFieldData(G4ThreeVector(0.0,0.0,0.0015)), Geant4.var"#getfield!#getfield!##0"(), Geant4.G4JLMagField[])Primary Particle Generator
@with_kw mutable struct PythiaGenerator <: G4JLGeneratorData
eCM::Float64 = 91.2*GeV
pythia = nothing
end
function LEPCollision()
data = PythiaGenerator()
function LEPinit(data::PythiaGenerator, app::G4JLApplication)
data.pythia = PYTHIA8.Pythia("", false)
data.pythia << "Beams:idA = -11" # positron
data.pythia << "Beams:idB = 11" # electron
data.pythia << "Beams:eCM = $(data.eCM/GeV)" # CM energy
data.pythia << "PDF:lepton = off" # no lepton PDFs
# Process: hadronic decays of Z0. Initialize. Histogram.
data.pythia << "WeakSingleBoson:ffbar2gmZ = on"
data.pythia << "23:onMode = off"
data.pythia << "23:onIfAny = 1 2 3 4 5"
# If Pythia fails to initialize, exit with error.
data.pythia |> init || error("Pythia initialization failed")
end
function LEPgenerate(g4evt::G4Event, data::PythiaGenerator)::Nothing
# Generate event. Skip if error. List first one.
data.pythia |> next || return
evt = data.pythia |> event
# Define primary vertex
vertex = G4PrimaryVertex(G4ThreeVector(0,0,0), 0)
# Loop over the particles in the Pythia event
for p in evt
if isFinal(p)
primary = G4PrimaryParticle(p |> PYTHIA8.id, px(p)*GeV, py(p)*GeV, pz(p)*GeV)
SetPrimary(vertex, move!(primary)) ## note that we give up ownership of the objects just created
end
end
AddPrimaryVertex(g4evt, move!(vertex)) ## note that we give up ownership of the objects just created
end
G4JLPrimaryGenerator("LEPCollision", data; init_method=LEPinit, generate_method=LEPgenerate)
endLEPCollision (generic function with 1 method)Event Display
Overwriting the default drawEvent function defined in the G4Vis extension
const G4Vis = Base.get_extension(Geant4, :G4Vis)
function G4Vis.drawEvent(evtdisp::G4Vis.G4JLEventDisplay)
s = evtdisp.lscene.scene
settings = evtdisp.settings
trajectories = G4EventManager!GetEventManager() |> GetConstCurrentEvent |> GetTrajectoryContainer
if trajectories != C_NULL
points = Point3{Float64}[]
for t in trajectories
t |> GetCharge == 0 && continue # skip neutral particles
npoints = GetPointEntries(t)
npoints < 3 && continue # skip particles with less than 3 points
for i in 1:npoints
pos = GetPoint(t, i-1) |> GetPosition
r = sqrt(pos[0]^2 + pos[1]^2)
z = pos[2]
(r > 1.8m || abs(z) > 2.4m) && break # kill particle outside the TPC
if r > 0.38m # skip points inside the beam pipe
push!(points, convert(Point3{Float64}, pos))
else
push!(points, Point3{Float64}(NaN, NaN, NaN))
end
end
push!(points, Point3{Float64}(NaN, NaN, NaN))
end
lines!(s, points, color=settings.trajectories.color, transparency=false, overdraw=false)
end
##wait_for_key("Press any key to continue with next event")
end
evtdisplay = G4JLEventDisplay(joinpath(@__DIR__, "TPCSimVisSettings.jl"))G4Vis.G4JLEventDisplay((display = (backgroundcolor = :black, resolution = (1280, 720), show_axis = false, camera_rotation = (0, 0, -1.5707963267948966), camera_zoom = 0.1), trajectories = (color = :yellow,), detector = (show_detector = true,)), G4Vis.stateChange, G4Vis.initDisplay, #undef, #undef)Create the Application
app = G4JLApplication(detector = tpc, # detector defined above
generator = LEPCollision(), # primary generator to instantiate
field = bfield, # magnetic field
physics_type = FTFP_BERT, # what physics list to instantiate
evtdisplay = evtdisplay # event display
);
configure(app)
initialize(app)
**************************************************************
Geant4 version Name: geant4-11-03-patch-02 [MT] (25-April-2025)
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/
**************************************************************
*------- PYTHIA Process Initialization --------------------------*
| |
| We collide e+ with e- at a CM energy of 9.120e+01 GeV |
| |
|------------------------------------------------------------------|
| | |
| Subprocess Code | Estimated |
| | max (mb) |
| | |
|------------------------------------------------------------------|
| | |
| f fbar -> gamma*/Z0 221 | 4.352e-05 |
| |
*------- End PYTHIA Process Initialization -----------------------*
*------- PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings (changes only) ------------------*
| |
| Name | Now | Default Min Max |
| | | |
| Beams:eCM | 91.20000 | 14000.000 0.0 |
| Beams:idA | -11 | 2212 |
| Beams:idB | 11 | 2212 |
| PDF:lepton | off | on |
| WeakSingleBoson:ffbar2gmZ | on | off |
| |
*------- End PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings -----------------------------*
-------- PYTHIA Particle Data Table (changed only) ------------------------------------------------------------------------------
id name antiName spn chg col m0 mWidth mMin mMax tau0 res dec ext vis wid
no onMode bRatio meMode products
23 Z0 3 0 0 91.18760 2.50419 10.00000 0.00000 7.87987e-14 1 1 0 0 0
0 1 0.1540492 0 1 -1
1 1 0.1194935 0 2 -2
2 1 0.1540386 0 3 -3
3 1 0.1193325 0 4 -4
4 1 0.1523269 0 5 -5
5 0 0.0335480 0 11 -11
6 0 0.0667305 0 12 -12
7 0 0.0335477 0 13 -13
8 0 0.0667305 0 14 -14
9 0 0.0334720 0 15 -15
10 0 0.0667305 0 16 -16
-------- End PYTHIA Particle Data Table -----------------------------------------------------------------------------------------
G4ChordFinder: stepperDriverId: 2
Run the Application
beamOn(app, 1)
-------- PYTHIA Info Listing ----------------------------------------
Beam A: id = -11, pz = 4.560e+01, e = 4.560e+01, m = 5.110e-04.
Beam B: id = 11, pz = -4.560e+01, e = 4.560e+01, m = 5.110e-04.
In 1: id = -11, x = 1.000e+00, pdf = 1.000e+00 at Q2 = 8.317e+03.
In 2: id = 11, x = 1.000e+00, pdf = 1.000e+00 at same Q2.
Subprocess f fbar -> gamma*/Z0 with code 221 is 2 -> 1.
It has sHat = 8.317e+03.
alphaEM = 7.818e-03, alphaS = 1.300e-01 at Q2 = 8.317e+03.
Impact parameter b = 0.000e+00 gives enhancement factor = 1.000e+00.
Max pT scale for MPI = 9.120e+01, ISR = 9.120e+01, FSR = 9.120e+01.
Number of MPI = 1, ISR = 0, FSRproc = 0, FSRreson = 5.
-------- End PYTHIA Info Listing ------------------------------------
-------- PYTHIA Event Listing (hard process) -----------------------------------------------------------------------------------
no id name status mothers daughters colours p_x p_y p_z e m
0 90 (system) -11 0 0 0 0 0 0 0.000 0.000 0.000 91.200 91.200
1 -11 (e+) -12 0 0 3 0 0 0 0.000 0.000 45.600 45.600 0.001
2 11 (e-) -12 0 0 4 0 0 0 0.000 0.000 -45.600 45.600 0.001
3 -11 (e+) -21 1 0 5 0 0 0 0.000 0.000 45.600 45.600 0.000
4 11 (e-) -21 2 0 5 0 0 0 0.000 0.000 -45.600 45.600 0.000
5 23 (Z0) -22 3 4 6 7 0 0 0.000 0.000 0.000 91.200 91.200
6 4 c 23 5 0 0 0 101 0 2.496 9.144 -44.579 45.600 1.500
7 -4 cbar 23 5 0 0 0 0 101 -2.496 -9.144 44.579 45.600 1.500
Charge sum: 0.000 Momentum sum: 0.000 0.000 0.000 91.200 91.200
-------- End PYTHIA Event Listing -----------------------------------------------------------------------------------------------
-------- PYTHIA Event Listing (complete event) ---------------------------------------------------------------------------------
no id name status mothers daughters colours p_x p_y p_z e m
0 90 (system) -11 0 0 0 0 0 0 0.000 0.000 0.000 91.200 91.200
1 -11 (e+) -12 0 0 3 0 0 0 0.000 0.000 45.600 45.600 0.001
2 11 (e-) -12 0 0 4 0 0 0 0.000 0.000 -45.600 45.600 0.001
3 -11 (e+) -21 1 0 5 0 0 0 0.000 0.000 45.600 45.600 0.000
4 11 (e-) -21 2 0 5 0 0 0 0.000 0.000 -45.600 45.600 0.000
5 23 (Z0) -22 3 4 6 7 0 0 0.000 0.000 0.000 91.200 91.200
6 4 (c) -23 5 0 10 10 101 0 2.496 9.144 -44.579 45.600 1.500
7 -4 (cbar) -23 5 0 8 9 0 101 -2.496 -9.144 44.579 45.600 1.500
8 -4 (cbar) -51 7 0 11 12 0 102 -2.892 -8.826 44.585 45.567 1.500
9 21 (g) -51 7 0 13 13 102 101 0.501 0.066 -1.877 1.944 0.000
10 4 (c) -52 6 6 23 23 101 0 2.391 8.760 -42.708 43.689 1.500
11 -4 (cbar) -51 8 0 19 19 0 103 -1.553 -7.359 38.768 39.519 1.500
12 21 (g) -51 8 0 14 15 103 102 -1.326 -1.466 5.769 6.098 0.000
13 21 (g) -52 9 9 16 16 102 101 0.488 0.064 -1.829 1.894 0.000
14 21 (g) -51 12 0 17 18 103 104 -0.900 -0.284 4.281 4.384 0.000
15 21 (g) -51 12 0 22 22 104 102 -0.393 -1.178 1.363 1.844 0.000
16 21 (g) -52 13 13 20 21 102 101 0.454 0.060 -1.704 1.764 0.000
17 21 (g) -51 14 0 31 31 105 104 -1.413 -1.644 9.897 10.132 0.000
18 21 (g) -51 14 0 32 32 103 105 0.276 -0.472 3.535 3.577 0.000
19 -4 (cbar) -52 11 11 33 33 0 103 -1.316 -5.528 29.616 30.194 1.500
20 -1 (dbar) -51 16 0 24 24 0 101 0.016 -0.556 -0.490 0.811 0.330
21 1 (d) -51 16 0 29 29 102 0 0.381 0.446 -1.017 1.219 0.330
22 21 (g) -52 15 15 30 30 104 102 -0.336 -1.007 1.166 1.577 0.000
23 4 (c) -71 10 10 25 28 101 0 2.391 8.760 -42.708 43.689 1.500
24 -1 (dbar) -71 20 20 25 28 0 101 0.016 -0.556 -0.490 0.811 0.330
25 433 (D*_s+) -83 23 24 52 53 0 0 1.594 5.437 -25.689 26.392 2.112
26 -313 (K*bar0) -83 23 24 42 43 0 0 0.162 1.084 -8.760 8.880 0.953
27 113 (rho0) -83 23 24 44 45 0 0 0.441 1.908 -6.675 6.988 0.662
28 223 (omega) -84 23 24 54 56 0 0 0.210 -0.226 -2.074 2.241 0.791
29 1 (d) -71 21 21 34 41 102 0 0.381 0.446 -1.017 1.219 0.330
30 21 (g) -71 22 22 34 41 104 102 -0.336 -1.007 1.166 1.577 0.000
31 21 (g) -71 17 17 34 41 105 104 -1.413 -1.644 9.897 10.132 0.000
32 21 (g) -71 18 18 34 41 103 105 0.276 -0.472 3.535 3.577 0.000
33 -4 (cbar) -71 19 19 34 41 0 103 -1.316 -5.528 29.616 30.194 1.500
34 -211 pi- 83 29 33 0 0 0 0 0.255 -0.125 -0.555 0.639 0.140
35 213 (rho+) -83 29 33 46 47 0 0 -0.023 -0.204 0.296 0.893 0.818
36 -213 (rho-) -83 29 33 48 49 0 0 -0.066 0.012 0.942 1.240 0.804
37 211 pi+ 83 29 33 0 0 0 0 -0.689 0.162 1.443 1.613 0.140
38 -211 pi- 83 29 33 0 0 0 0 0.568 -0.525 1.223 1.454 0.140
39 211 pi+ 84 29 33 0 0 0 0 -0.896 -0.301 3.716 3.837 0.140
40 313 (K*0) -84 29 33 50 51 0 0 -0.363 -2.802 11.177 11.563 0.887
41 -433 (D*_s-) -84 29 33 57 58 0 0 -1.193 -4.421 24.956 25.460 2.112
42 -321 K- 91 26 0 0 0 0 0 -0.209 0.443 -4.499 4.552 0.494
43 211 pi+ 91 26 0 0 0 0 0 0.371 0.641 -4.261 4.328 0.140
44 211 pi+ 91 27 0 0 0 0 0 0.047 0.795 -1.860 2.028 0.140
45 -211 pi- 91 27 0 0 0 0 0 0.394 1.113 -4.815 4.960 0.140
46 211 pi+ 91 35 0 0 0 0 0 0.209 0.120 0.373 0.465 0.140
47 111 (pi0) -91 35 0 59 60 0 0 -0.232 -0.324 -0.077 0.428 0.135
48 -211 pi- 91 36 0 0 0 0 0 0.182 -0.293 0.582 0.691 0.140
49 111 (pi0) -91 36 0 61 62 0 0 -0.247 0.305 0.360 0.550 0.135
50 321 K+ 91 40 0 0 0 0 0 -0.222 -2.661 10.127 10.484 0.494
51 -211 pi- 91 40 0 0 0 0 0 -0.141 -0.141 1.051 1.079 0.140
52 431 (D_s+) -91 25 0 63 64 0 0 1.466 4.909 -22.785 23.436 1.968
53 22 gamma 91 25 0 0 0 0 0 0.129 0.528 -2.905 2.955 0.000
54 211 pi+ 91 28 0 0 0 0 0 0.172 -0.123 -0.906 0.941 0.140
55 -211 pi- 91 28 0 0 0 0 0 -0.095 -0.227 -0.927 0.969 0.140
56 111 (pi0) -91 28 0 65 66 0 0 0.134 0.125 -0.242 0.332 0.135
57 -431 (D_s-) -91 41 0 67 68 0 0 -1.224 -4.222 24.350 24.822 1.968
58 22 gamma 91 41 0 0 0 0 0 0.031 -0.199 0.606 0.638 0.000
59 22 gamma 91 47 0 0 0 0 0 -0.134 -0.231 0.009 0.268 0.000
60 22 gamma 91 47 0 0 0 0 0 -0.098 -0.093 -0.087 0.160 0.000
61 22 gamma 91 49 0 0 0 0 0 -0.012 0.102 0.070 0.125 0.000
62 22 gamma 91 49 0 0 0 0 0 -0.235 0.202 0.290 0.425 0.000
63 333 (phi) -91 52 0 69 70 0 0 0.819 2.554 -10.484 10.869 1.020
64 213 (rho+) -91 52 0 71 72 0 0 0.646 2.355 -12.301 12.567 0.808
65 22 gamma 91 56 0 0 0 0 0 -0.024 -0.002 -0.014 0.028 0.000
66 22 gamma 91 56 0 0 0 0 0 0.157 0.127 -0.227 0.304 0.000
67 -321 K- 91 57 0 0 0 0 0 -0.217 -1.445 11.682 11.783 0.494
68 313 (K*0) -91 57 0 73 74 0 0 -1.008 -2.776 12.669 13.039 0.893
69 130 K_L0 91 63 0 0 0 0 0 0.428 1.238 -5.537 5.711 0.498
70 310 (K_S0) -91 63 0 75 76 0 0 0.391 1.316 -4.947 5.158 0.498
71 211 pi+ 91 64 0 0 0 0 0 0.432 0.526 -2.557 2.649 0.140
72 111 (pi0) -91 64 0 77 78 0 0 0.215 1.829 -9.744 9.918 0.135
73 321 K+ 91 68 0 0 0 0 0 -1.046 -2.589 11.320 11.670 0.494
74 -211 pi- 91 68 0 0 0 0 0 0.038 -0.187 1.348 1.369 0.140
75 211 pi+ 91 70 0 0 0 0 0 0.361 1.030 -3.250 3.431 0.140
76 -211 pi- 91 70 0 0 0 0 0 0.030 0.286 -1.697 1.727 0.140
77 22 gamma 91 72 0 0 0 0 0 0.045 0.122 -0.779 0.790 0.000
78 22 gamma 91 72 0 0 0 0 0 0.170 1.707 -8.965 9.128 0.000
Charge sum: 0.000 Momentum sum: 0.000 0.000 -0.000 91.200 91.200
-------- End PYTHIA Event Listing -----------------------------------------------------------------------------------------------
display(evtdisplay.figure)
Another event
beamOn(app, 1)display(evtdisplay.figure)
This page was generated using Literate.jl.