Basic python geometry scripting

Precepts

• Units may be specified but default to Geant4 ones (e.g. mm, rad).

• Rotations are made using Tait-Bryan angles (rotation about reference x,y,z axes).

• A Registry object should be used to hold all things related in a model and passed into the constructors of most objects.

• GDML-like full lengths are used instead of typically half lengths

Units

• The default units are mm and rad for length and angle.

• Most constructors will take units as an optional key word argument (‘kwarg’).

• The kwargs are typically lunit for length unit and aunit for angle unit.

• Units are defined in pyg4ometry.gdml.Units.py.

The following units (as strings) are accepted:

Unit	Value
nm	1e-6
um	1e-3
mm	1
cm	10
m	1e3
km	1e6
deg	\(\pi/180\)
degree	\(\pi/180\)
rad	1
radian	1
mrad	1e-3
urad	1e-6
eV	1e-3
keV	1
MeV	1e+3
none	1
ns	1e-9
us	1e-6
ms	1e-3
s	1

Examples:

reg = pyg4ometry.geant4.Registry()
boxSolid = pyg4ometry.genat4.solid.Box("aBox", 10, 20, 30, reg)

This defines a box with the default units (none specifed), so mm. We can specify them:

boxSolid = pyg4ometry.genat4.solid.Box("aBox", 10, 20, 30, reg, "cm")

Geant4 Python Scripting

Making use of pyg4ometry requires the following modules

import pyg4ometry

To make a simple geometry of a box located at the origin

# load pyg4ometry
import pyg4ometry

# registry to store gdml data
reg  = pyg4ometry.geant4.Registry()

# world solid and logical
ws   = pyg4ometry.geant4.solid.Box("ws",50,50,50,reg)
wl   = pyg4ometry.geant4.LogicalVolume(ws,"G4_Galactic","wl",reg)
reg.setWorld(wl.name)

# box placed at origin
b1   = pyg4ometry.geant4.solid.Box("b1",10,10,10,reg)
b1_l = pyg4ometry.geant4.LogicalVolume(b1,"G4_Fe","b1_l",reg)
b1_p = pyg4ometry.geant4.PhysicalVolume([0,0,0],[0,0,0],b1_l,"b1_p",wl,reg)

# visualise geometry
v = pyg4ometry.visualisation.VtkViewer()
v.addLogicalVolume(wl)
v.addAxes(20)
v.view()

Here is the vtk visualiser output of the above example

GDML Defines

In GDML there are multiple define objects that can be used parameterise geometry, materials etc. These can be used as variables or definitions and mean that any equations used will be retained in GDML output. For example a GDML constant can be created in the following way

# registry to store gdml data
reg = pyg4ometry.geant4.Registry()

# constant called x
x = pyg4ometry.gdml.Constant("x",10,reg)

The normal set of mathematical operations in python can be performed and evaluated

y = 2*x + 10
y.eval()

>> 30

The constant x can of course be changed and y re-evaluated

x.setExpression(5)
y.eval()

>> 20

Note

Standard mathematical functions can be used with GDML defines (Constant, Variable, etc). So sin, cos, tan, exp and so on, but pyg4ometry functions have to be used

x  = pyg4ometry.gdml.Constant("x",10,reg)
cx = pyg4ometry.gdml.cos(x)

So the box example above can be rewritten using constants

# load pyg4ometry
import pyg4ometry

# registry to store gdml data
reg  = pyg4ometry.geant4.Registry()

bx = pyg4ometry.gdml.Constant("bx","10",reg,True)
by = pyg4ometry.gdml.Constant("by",2*bx,reg,True)
bz = pyg4ometry.gdml.Constant("bz",2*by,reg,True)

# world solid and logical
ws   = pyg4ometry.geant4.solid.Box("ws",50,50,50,reg)
wl   = pyg4ometry.geant4.LogicalVolume(ws,"G4_Galactic","wl",reg)

# box placed at origin
b1   = pyg4ometry.geant4.solid.Box("b1",bx,by,bz,reg)
b1_l = pyg4ometry.geant4.LogicalVolume(b1,"G4_Fe","b1_l",reg)
b1_p = pyg4ometry.geant4.PhysicalVolume([0,0,0],[0,0,0],b1_l,"b1_p",wl,reg)

# visualise geometry
v = pyg4ometry.visualisation.VtkViewer()
v.addLogicalVolume(wl)
v.addAxes(20)
v.view()

Note

All GDML defines (Constant, Variable, etc) can be used in the construction of other pyg4ometry classes interchangeably instead of floats or strings (where strings are either numbers or a GDML expression)

Warning

Avoid reassigning variables used as defines, this can have unexpected consequences so for example

b1   = pyg4ometry.geant4.solid.Box("b1",bx,by,bz,reg)
b1.pX = 20              # do not do this
b1.pX.setExpression(20) # rather do this

Solids

The python geant4 solids match the Geant4 constructors as much possible (different constructor signatures are not supported in python). For example looking at the G4Box class

pyg4ometry.geant4.solid.Box(name, pX, pY, pZ, registry, lunit)

G4Box(const G4String& pName, G4double  pX, G4double  pY, G4double pZ)

A full list of solids can be found in Geant4 solids.

Warning

The parameters stick to the GDML convention of full lengths opposed to half lengths.

Materials

As with solids materials are defined in a similar way to Geant4 C++. Python does not have overloaded constructors, so unique signatures are needed, in contrast to Geant4.

To define a material from the Geant4 predefined (e.g. NIST) materials

import pyg4ometry.geant4 as _g4
wm = _g4.MaterialPredefined("G4_Galactic")
bm = _g4.MaterialPredefined("G4_Fe")

To define a single element in terms of atomic number, atomic mass and density.

import pyg4ometry.geant4 as _g4
wm = _g4.MaterialSingleElement("galactic",1,1.008,1e-25,reg)   # low density hydrogen
bm = _g4.MaterialSingleElement("iron",26,55.8452,7.874,reg)    # iron at near room temp

To define a compound two elements using the mass fraction

import pyg4ometry.geant4 as _g4
wm = _g4.MaterialCompound("air",1.290e-3,2,reg)
ne = _g4.ElementSimple("nitrogen","N",7,14.01)
oe = _g4.ElementSimple("oxygen","O",8,16.0)
wm.add_element_massfraction(ne,0.7)
wm.add_element_massfraction(oe,0.3)
bm = _g4.MaterialSingleElement("iron",26,55.8452,7.874,reg)    # iron at near room temp

To define a compound using number of atoms

import pyg4ometry.geant4 as _g4
bm = _g4.MaterialCompound("plastic",1.38,3,reg)    # Generic PET C_10 H_8 O_4
he = _g4.ElementSimple("hydrogen","H",1,1.008)
ce = _g4.ElementSimple("carbon","C",6,12.0096)
oe = _g4.ElementSimple("oxygen","O",8,16.0)
bm.add_element_natoms(he,8)
bm.add_element_natoms(ce,10)
bm.add_element_natoms(oe,4)

Material as a mixture of materials

import pyg4ometry.geant4 as _g4
bm     = _g4.MaterialCompound("YellowBrass_C26800", 8.14, 2, reg)
copper = _g4.MaterialPredefined("G4_Cu")
zinc   = _g4.MaterialPredefined("G4_Zn")
bm.add_material(copper, 0.67)
bm.add_material(zinc, 0.33)

Example of elements formed by isotopes

import pyg4ometry.geant4 as _g4
u235 = _g4.Isotope("U235", 92, 235, 235.044)
u238 = _g4.Isotope("U238", 92, 238, 238.051)
uranium = _g4.ElementIsotopeMixture("uranium", "U", 2)
uranium.add_isotope(u235, 0.00716)
uranium.add_isotope(u238, 0.99284)
bm = _g4.MaterialCompound("natural_uranium", 19.1, 1, reg)
bm.add_element_massfraction(uranium, 1)

NIST Materials

Geant4 has many predefined materials according to the NIST database. Their name typically starts with G4_. These typically can be used with MaterialPredefined and we do not need to specify the full composition - Geant4 will find them at run time.

However, in the case of conversion to FLUKA, these are fully expanded according to their definition in Geant4 based on a cache in pyg4ometry of the material compositions generated using BDSIM from Geant4 (10.7.p01 as of writing). Should the user wish to use these, they can be accessed from the functions in the geant4 module.

import pyg4ometry
nistHydrogenElement = pyg4ometry.geant4.nist_element_2geant4Element('G4_H')

Note, an ‘element’ cannot be used as a ‘material’ in a logical volume. We must upgrade it to a material for that. The NIST elements contain the appropriate mixture of natural isotopes and can be used in MaterialCompound as demonstrated above.

Alternatively, we can access the NIST materials and materials of elements.

import pyg4ometry
nistHydrogenMaterial = pyg4ometry.geant4.nist_material_2geant4Material('G4_H')
nistConcreteMaterial = pyg4ometry.geant4.nist_material_2geant4Material('G4_CONCRETE')

Detector Construction

This largely proceeds in exactly the same way as in G4 or GDML. Hierarchy of solids, booleans, logical, physical (replica, division, param) volumes.

0. Create registry to hold everything

1. Create solids

2. Create logical volumes

3. Place logical volumes (construct physical volumes)

4. Visualise

5. Check

6. Export

Transformations & Physical Volumes

Transformations in 3D are essential for the easy placement of solids in a CSG tree or LV placement. There is not a specific transformation class like in Geant4. The matrices and vectors used for placements are here typically Numpy arrays or matrices.

Geant4 has two possible constructors for a physical volume. These provide active and passive transformations. In pyg4ometry, only one is provided.

• The transform in a physical volume first translates the placed logical volume with respect to the mother logical, then rotates it.

The physical volume class is documented here: Geant4 module, but an example is shown here.

import pyg4ometry
r = pyg4ometry.geant4.Registry()
vacuum = _g4.MaterialPredefined("G4_Galactic")
water = _g4.MaterialPredefined("G4_WATER")
worldSolid = pyg4ometry.geant4.solid.Box("world_solid", 100, 100, 100, reg)
boxSolid = pyg4ometry.geant4.solid.Box("box_solid", 10, 20, 40, reg)
worldLV = pyg4ometry.geant4.LogicalVolume(worldSolid, vacuum, "world_lv", reg)
boxLV = pyg4ometry.geant4.LogicalVolume(boxSolid, water, "box_lv", reg)

pyg4ometry.geant4.PhysicalVolume([0,0,0],
                                 [0,0,0],
                                 boxLV,
                                 "box_pv",
                                 worldLV,
                                 reg)

This creates a box of water inside a box of vacuum. The box of water is 10 x 20 x 50 mm long (note mm are the default length units), and it is placed with no offset and no rotation (i.e. at the centre) of the world volume. Alternatively:

import numpy as np
pyg4ometry.geant4.PhysicalVolume([0,np.pi/3.0,0],
                                 [0,0,0],
                                 boxLV,
                                 "box_pv",
                                 worldLV,
                                 reg)

In this case, the box is placed with no offset but with a rotation of \(\pi/3\) radians about the y axis of the world box.

Note

The rotations are Tait-Bryan angles, which are rotations about the reference x,y,z axes. i.e. if there is a rotation about both x and y, these are independent and it is not a compound frame that is rotated. These are commonly thought of like an aircraft and called pitch, yaw and tilt.

There are utility functions for translation between different transformations in pyg4ometry.transformation. See Transformation.

Optical Surfaces

Optical surfaces can be created in a similar way as in Geant4 C++. A pyg4ometry.geant4.solid.OpticalSurface instance holds all the needed properties of the surface (including extra properties, e.g. for optical processes). This is then assigned to the surface between either

• two physical volumes: pyg4ometry.geant4.BorderSurface, or

• a logical volume and all its neighbouring volumes: pyg4ometry.geant4.SkinSurface.

opa = _g4.solid.OpticalSurface("AirSurface", finish="polished", model="glisur", surf_type="dielectric_dielectric", value="1", registry=reg)
opw = _g4.solid.OpticalSurface("WaterSurface", finish="ground", model="unified", surf_type="dielectric_dielectric", value="0", registry=reg)

_g4.SkinSurface("AirSurface", air_lv, opa, reg)
_g4.BorderSurface("WaterSurface", water_phys, world_phys, opw, reg)

Properties of Materials and Optical Surfaces

Materials and optical surfaces support adding properties that can be used by Geant4 to influence processes, e.g. for scintillation, refraction or other optical processes.

In the GDML, a matrix is used to hold the value(s) of the property.

• addProperty(name, matrix) - Add a property based on an existing pyg4ometry.gdml.Matrix object.

• addVecProperty(name, e, v, eunit='eV', vunit='') - Add a property based on a energy vector and a value vector.

• addConstProperty(name, value, vunit='')- Add a property that has only one constant value.

Units can be specified by setting the parameters eunit for the energy vector and vunit for the values. The given vectors are expected to be homogeneous in their units.

Note

Optical properties can only use units (or combinations of units) that are also defined in pyg4ometry. If needed, additional units can be added: pyg4ometry.gdml.Units.units['ps'] = 1e-12.

scint = _g4.Material(...)
scint.addConstProperty('SCINTILLATIONTIMECONSTANT1', 2.5, vunit='ns')
scint.addConstProperty('SCINTILLATIONYIELD', 8000, vunit='/MeV')
scint.addVecProperty('RINDEX', [1, 10], [1.3, 1.05])

Registry and GDML Output

Strictly speaking a registry class to store all of the GDML is not required. As with normal Geant4 given a lv pointer it should possible to form an aggregation hierarchy that contains all necessary objects. Now GDML breaks this as the structure is built up using name tags. For example a placement requires a position. In Geant4 this would just be a pointer to an transformation object, but GDML has two mechanisms to represent this, firstly child nodes of a PhysicalVolume tag or secondly a position define, see below

The registry class is a storage class for a complete GDML file. At the construction stage of almost all objects a registry is required. If the object is added to the registry then it will appear explicitly in the GDML output

Visualisation

Any logical volume lv can be visualised using:

v = pyg4ometry.visualisation.VtkViewer()
v.addLogicalVolume(lv)
v.addAxes(20)
v.view()

which will open a Vtk render window. The render window now receives keyboard and mouse commands. To exit render window q, to restart interaction with the visualiser

v.start()

There are also convenience methods of pyg4ometry.visualisation.VtkViewer() the allow changing of the viewing parameters. So if the viewer is active then render window needs to be stopped q and then commands can be typed into the terminal, for example

v.setOpactity(0.1)
v.setWirefrace()
v.start()

Overlap Checking

“Overlaps” is a general term used to describe malformed geometry. Such geometry is unphysical and may causing particle tracking problems in simulations such as stuck particles, or particles completely missing certain volumes entirely. Such errors are rarely easy to spot from results or running the simulation.

Given all the PVs (daughters) of a LV (mother) should be bounded by the LV/mother solid. It is possible to check between all daughter solid meshes and between daughters and the mother solid mesh. Given an pyg4ometry.geant4.LogicalVolume instance (“lv”), this check can be performed by calling the following code:

lv.checkOverlaps()

This will check only the immediate daughters of this logical volume. To descend further into a geometry, the recursive flag can be used:

lv.checkOverlaps(recursive=True)

See Geant4 module : LogicalVolume.checkOverlaps() for full details. A more complete example is:

# cd pyg4ometry/pyg4ometry/test/pythonGeant4
import pyg4ometry
r  = pyg4ometry.freecad.Reader("./T103_overlap_copl.gdml")
l = r.getRegistry().getWorldVolume()
l.checkOverlaps(recursive=False,coplanar=True,debugIO=False)
v = pyg4ometry.visualisation.VtkViewer()
v.addLogicalVolume(l)
v.view()

Text is by default only printed out when an overlap is found. Any overlaps will be prepared for visualisation in a VtkViewer (must be constructed and given the LV after this).

The following overlap checks are performed:

1. daughter with other daughter overlap

2. co-planar daughter with other daughter overlap

3. protrusion of a daughter from the mother volume

4. co-planar daughter with mother volume

Colour Coding

In the visualiser, text will be overlaid saying “overlap” where some kind of overlap is detected. Additionally, the actual overlap itself will be visualised and colour coded according to:

• red: protrusion overlap

• green: daughter-daughter overlap

• blue: co-planar overlap

Limitations

1. The overlap detection is performed by checking for overlaps in the visualisation meshes generated for each volume. In the case of curved solids (e.g. a cylinder), the mesh is not truly curved but a polygon. Very closely spaced curved surfaces may produce false overlaps. By default, all curved solids will use the same number of points around a circle, so usually we can “get away” with this if the curved solids aren’t rotated about their axis.

2. Currently, division and parameterised volumes are not handled explicitly.

Assemblies

In the case of assembly volumes, and if an overlap is detected, a unique name is built up based on the parent PhysicalVolume, the assembly and the PhysicalVolume inside it. Furthermore, this is done recursively is assemblies of assemblies (etc) are used. The name is built up with an underscore “_” for padding and the user should decode this from their input.

As there is no ‘mother’ of an assembly, there is no mother protrusion directly. The contents of an assembly are compared to all other daughters and the mother at the higher level in which they are placed.

GDML Output

To write an GDML file file given a pyg4ometry.geant4.registy instance reg.

import pyg4ometry
w = p4gometry.gdml.Writer()
w.addDetector(reg)
w.write('file.gdml')
# make a quick bdsim job for the one component in a beam line
w.writeGmadTester('file.gmad', 'file.gdml')