Defining Primary Particles#
A Geant4 application requires a source of primary particles. We have several possibilities to provide a primary particle generator.
Particle Gun - This is the simplest primary particle generator. It is a single type of particle with a fix energy and direction
General Particle Source - This is a generalization of a particle gun with many more parameters to configure its behavior. See the G4GeneralParticleSource documentation.
Full customization - The user can write code to decide the primary particles to create, as well as, their energies and directions for each event.
We are going to exercise the three cases in this notebook
using Geant4
using Geant4.SystemOfUnits
# Do some analysis of energy, positions and directions
using DataFrames, Plots
function generate(generator, N)
df = DataFrame(energy = Float32[],
pos_x = Float32[], pos_y = Float32[], pos_z = Float32[],
dir_x = Float32[], dir_y = Float32[], dir_z = Float32[])
for i in 1:N
evt = G4Event()
generator.gen_method(evt, generator.data)
pos = evt |> GetPrimaryVertex |> GetPosition
dir = evt |> GetPrimaryVertex |> GetPrimary |> GetMomentumDirection
ene = evt |> GetPrimaryVertex |> GetPrimary |> GetKineticEnergy
push!(df, (ene, x(pos), y(pos), z(pos), x(dir), y(dir), z(dir)))
end
return df
end
generate (generic function with 1 method)
Particle Gun#
We have introduced a new type G4JLGunGenerator to facilitate the task of defining a
particlegun = G4JLGunGenerator(particle = "pi+",
energy = 330MeV,
direction = G4ThreeVector(0, 0, 1),
position = G4ThreeVector(0, 0, 0));
The actual G4ParticleGun object (the data.gun field) has not been created yet. It will be done by the G4RunManager during initialization. If we do initialize by hand now we will get a error since particles have not been created yet.
Lets build an application with just the generator
app = G4JLApplication(
generator = particlegun,
)
configure(app)
initialize(app)
**************************************************************
Geant4 version Name: geant4-11-02-patch-01 [MT] (16-February-2024)
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/
**************************************************************
Run one event and see the printout of the initial particle
ui`/tracking/verbose 1`
beamOn(app,1)
*********************************************************************************************************
* G4Track Information: Particle = pi+, Track ID = 1, Parent ID = 0
*********************************************************************************************************
Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName
0 0 0 0 330 0 0 0 World initStep
1 0 0 1e+03 330 2.73e-23 1e+03 1e+03 OutOfWorld Transportation
We can easily change the parameters of the particle gun with the following methods
SetParticleByName(particlegun, "e-")
SetParticleEnergy(particlegun, 10GeV)
SetParticleMomentumDirection(particlegun, G4ThreeVector(1,0,0))
beamOn(app,1)
*********************************************************************************************************
* G4Track Information: Particle = e-, Track ID = 1, Parent ID = 0
*********************************************************************************************************
Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName
0 0 0 0 1e+04 0 0 0 World initStep
1 1e+03 0 0 1e+04 5.12e-23 1e+03 1e+03 OutOfWorld Transportation
General Particle Source#
You can create a G4JLGeneralParticleSource with the same parameters (and same behavior) as the Particle Gun.
gps = G4JLGeneralParticleSource(particle = "proton",
energy = 10MeV,
direction = G4ThreeVector(1,0,0),
position = G4ThreeVector(0,0,0));
app = G4JLApplication(
generator = gps,
)
configure(app)
initialize(app)
df = generate(gps, 1000)
histogram(df.energy)
Note that The exception is due to the fact that Geant4 does not like to re-initialize. Please ignore the exception and look at the created particle. Now we can change the parameters in two different ways:
using the UI commands
instantiating a new
G4JLGeneralParticleSourcewith the needed parameters
reinitialize(app.generator; particle="geantino",
energy=2MeV,
pos=(type="Point", centre=G4ThreeVector(1cm,2cm,1cm)),
ang=(type="iso"))
df = generate(gps, 1000)
histogram(df.energy)
reinitialize(gps; particle="geantino",
ene=(type="Lin", min=2MeV, max=10MeV, gradient=1, intercept=1),
pos=(type="Plane", shape="Square", centre=G4ThreeVector(1cm,2cm,1cm), halfx=2cm, halfy=2cm),
ang=(type="cos",mintheta=10deg, maxtheta=80deg))
df = generate(gps, 1000)
histogram(df.energy)
scatter(df.pos_x, df.pos_y, df.pos_z)
Custom Primary Generator#
The user can provide fully custom primary particle generator. This is done by defining custom structure to configure the generator and two functions to initialize and generate the primary particles that will be called by the G4RunManager for each event.
The user configuration data structure should inherit from G4JLGeneratorData abstract type. Here is an example of a generator of a for a rectangle shape origined a mono-energy particles.
# define the data structure with the generator parameters
mutable struct PlaneSourceData <: G4JLGeneratorData
particleName::String # particle type
particlePtr::CxxPtr{G4ParticleDefinition} # keep the pointer to definition for performance
energy::Float64 # kinetic energy
halfx::Float64 # rectangle dimensions
halfy::Float64
position::G4ThreeVector # rectangle origin
direction::G4ThreeVector # particle direction
end
We defile now a function PlaneSource that will create an instance of G4JLPrimaryGenerator with the configuration data (using some default values) and the two functions: init and generate. The init converts the particle name to a definition, which is then used for the generate function at each event. The generate generate random points en X and Y and creates a G4PrimaryParticle and G4PrimaryVertex which are added to the event given as argument.
# define the constructor with the default parameters
function PlaneSource(;particle="gamma", energy=0.07MeV, halfx=7cm, halfy=7cm,
position=G4ThreeVector(0,0,0),
direction=G4ThreeVector(0,0,1))
data = PlaneSourceData(particle, CxxPtr{G4ParticleDefinition}(C_NULL), energy, halfx, halfy, position, direction)
function init(data:: PlaneSourceData, app::G4JLApplication)
data.particlePtr = FindParticle(data.particleName)
end
function generate( evt::G4Event, data:: PlaneSourceData)::Nothing
mass = data.particlePtr |> GetPDGMass
momentum = √((mass + data.energy)^2 - mass^2)
pvec = momentum * data.direction
pos = data.position + G4ThreeVector( data.halfx * (rand() - 0.5), data.halfy * (rand() - 0.5), 0)
primary = G4PrimaryParticle(data.particlePtr, pvec |> x, pvec |> y, pvec |> z )
vertex = G4PrimaryVertex(pos, 0ns)
SetPrimary(vertex, move!(primary)) # note that we give up ownership of the objects just created
AddPrimaryVertex(evt, move!(vertex)) # note that we give up ownership of the objects just created
end
G4JLPrimaryGenerator("PlaneSource", data; init_method=init, generate_method=generate)
end
PlaneSource (generic function with 1 method)
planesource = PlaneSource(energy=10MeV, halfx=10cm, halfy=10cm)
G4JLPrimaryGenerator{PlaneSourceData}("PlaneSource", PlaneSourceData("gamma", CxxPtr{G4ParticleDefinition}(Ptr{G4ParticleDefinition} @0x0000000000000000), 10.0, 100.0, 100.0, G4ThreeVector(0.0,0.0,0.0), G4ThreeVector(0.0,0.0,1.0)), var"#init#2"(), var"#generate#3"(), G4JLGeneratorAction[])
app = G4JLApplication(
generator = planesource,
)
configure(app)
initialize(app)
df = generate(planesource, 1000)
histogram(df.energy)
scatter(df.pos_x, df.pos_y)
scatter(df.pos_x, df.pos_y, df.pos_z)