A half-space is defined by a surface f x , y , z = 0 , and a sense ϵ indicating the side of the surface where the half-space is located. The ϵ ∈ − 1 , 1 is either 1 or −1 when the half-space is on the positive side ( f x , y , z > 0 ) or negative side ( f x , y , z < 0 ) of the surface.

Considering a half-space is created from a surface f x , y , z = 0 , the position of a point x , y , z related to the half-space can be indicated by Equation 1 S x , y , z = ϵ · f x , y , z < 0 , if outside = 0 , if on the surface  > 0 , if inside  . (1)

A point is inside a half-space or on the surface only if it is within the semi-algebraic subset P ⊂ R 3 (Bochnak et al., 2013), which P is defined in Equation 2: P = x , y , z ∈ R 3 | S x , y , z ≥ 0   . (2)

A half-space solid is bounded by a finite number of half-spaces using the Boolean intersection operation. Therefore, a point is inside or on the boundaries of a half-space solid only if it is inside the closed subset Q ⊂ R 3 : Q = ⋂ i n P i = ⋂ i n x , y , z ∈ R 3 | S i x , y , z ≥ 0 , (3) where n denotes the number of half-spaces, and i = 1 , 2 , 3 , … n denotes the index of half-spaces.

Under a subsequence of decompositions, it is proved that a BRep CAD solid can be converted to a closed semi-algebraic subset R ⊂ R 3 (Tsige-Tamirat, 2001): R = ⋃ j m ⋂ i n x , y , z ∈ R 3 | S j , i x , y , z ≥ 0 , (4) where S j , i x , y , z is defined in (1), m denotes the number of decomposed solids, and j = 1 , 2 , 3 , … m . Substituting (3) into (4) proves that a BRep solid can be represented by unions of finite numbers of half-space solids, R, defined in Equation 5: R = ⋃ j m Q j (5)

The half-space solid is a specific type of half-space CSG, and it is the key to the conversion of CAD geometries to CSG. By implementing this new solid type in Geant4, the modeling of CAD geometry becomes feasible.

The half-space solid type has been developed as a new Geant4 class called G4HalfSpaceSolid by inheriting from the base class G4VSolid. As shown in Figure 1, it has a surface list containing the objects of the half-space class G4HalfSpaceSurface and is implemented with mandatory functions for geometry tracking. G4HalfSpaceSurface is a base class for the classes of specific surface types that restore surface parameters and implement the intersection calculation. Currently, six surface types are supported: plane, cylinder, sphere, cone, general quadric, and torus.

The G4HalfSpaceSolid implements the functions of the G4VSolid interface that are necessary for ray tracking and navigation. The function Inside checks the relative position of a point, that is, inside, outside, or on the boundaries, regarding the solid. The functions DistanceToIn and DistanceToOut compute the distance of a particle entering or exiting this solid along a given direction. The functions SafetyToIn and SafetyToOut estimate the safety (underestimated) distance from a point to enter or exit the solid in any direction.

In the implementation of the function Inside, the bounding box of the solid with a small margin on each side is adopted for a first check of the point position. A point is surely outside the half-space solid when it is outside the bounding box. The points inside the bounding box will be further confirmed by the half-spaces of the solid. Based on (3), a point is outside the half-space solid when it is inside the complement of Q defined in Equation 6: Q ¯ = ⋃ i n P ¯ i   . (6)

Therefore, if S i x , y , z < 0 is true for one of the half-spaces P i , the point x , y , z is outside the solid. If not, the point is on the boundary when one half-space has S i x , y , z = 0 , or it is inside the solid when all half-spaces have S i x , y , z > 0 .

The key functions of G4HalfSpaceSolid are the particle-geometry intersection calculation in DistanceToIn and DistanceToOut. A neutral particle going in a straight trajectory is equivalent to a ray, which has a starting point and a direction vector. For charged particles, the curve trajectory is approximated by chords (Allison et al., 2016), and the intersection calculation method is hence similar to that of a straight trajectory. The intersections of a ray with all boundary surfaces of a half-space solid are performed, and the intersection points outside the solid are discarded because they are not valid. Normally, the distance of the ray to the nearest intersection point is the actual intersection distance. However, some exceptions require special caution. For example, in Figure 2, ray a crosses exactly the intersection line of the cylinder and the plane at Point-1 and enters the solid at Point-2. Point-2 is the correct entry point, but Point-1 might be mistakenly taken because it is the nearest valid answer. To avoid this issue, a probe point with a tiny distance behind Point-1 must be calculated. If the particle is outside the solid after walking in this probe step, Point-2 is the correct answer instead of Point-1.

Another issue is caused by solution tolerance. For example, an intersection point x 0 , y 0 , z 0 is supposed to have f x 0 , y 0 , z 0 = 0 , whereas f x 0 , y 0 , z 0 might be very closed but not equal to 0. In G4HalfSpaceSolid, a tolerance ε (e.g., 10⁻⁹) is adopted, and a point that satisfies f x 0 , y 0 , z 0 < ε is considered on the boundary surfaces. The direction of the particle is important for dealing with the precision problem. For example, in Figure 2, ray b enters the gray solid at Point-3 from the adjacent solid, but it could be already inside the solid because of the precision problem. An invalid answer will be produced in the function DistanceToIn, and the particle will travel blindly through this solid. Therefore, special treatment is performed to handle this issue. Assuming the normal of the boundary surface is always pointing outward the solid, the cosine of the angle θ between the direction of the ray and the surface normal is computed. If cos ⁡ θ < 0 , Point-3 is accepted as the entering point of the ray b, even though it is not on the trajectory.

The safety distances computed by the functions SafetyToIn and SafetyToOut are used to accelerate the geometry tracking. If the distance of a particle to the next collision site is smaller than the safety distance, the collision simulation is invoked without knowing the exact intersection distance to solid boundaries. Therefore, the safety distance must be smaller than the actual intersection distance, and the computational effort of calculating the safety distance should be considerably smaller as well. In SafetyToIn, an underestimated distance from a particle to the bounding box is computed. In SafetyToOut, the minimum isotropic distance of the particle to the boundary surfaces of the half-space solid is calculated. However, calculating the isotropic distance is very computationally expensive for the general quadric and torus surfaces. In these cases, the safety distances are set to 0. This conservative estimation does not affect the correctness of the actual simulation because the intersection calculation will be invoked.

With all these implementations, the G4HalfSpaceSolid is intended to be a kind of generic-purpose solid type for use in neutral or charged particle simulation. To create a G4HalfSpaceSolid in Geant4, the attributes include the bounding box, volume size, and surface area. In addition, the polyhedron is computed and provided in the later modeling method so that all the solids can be visualized in the Geant4 native visualization system. Manually creating a half-space solid is not straightforward; thus, a modeling approach has been developed in this work to generate a half-space solid model automatically from CAD geometries.

GDML is an application-independent geometry description format based on XML, and it is used as a general format in Geant4 for geometry persistence. It has self-consistent definitions of syntax that are specified in the GDML schema. Therefore, GDML can be extended to describe new solid types as long as the GDML schema is extended. In this work, the description of the half-space solid has been integrated into the GDML schema. As defined in the schema, a GDML file consists of five blocks—Define, the Material, the Solids, the Structure, and the Setup blocks. Extensions have been implemented in the Define and Solids blocks, as shown in Figure 3.

In the Define block, the Polyhedron is a polygon mesh used for visualizing the half-space solid. It is built by a list of nodes defined by three coordinates and a list of triangle facets consisting of three-node indices. In the Solids block, a half-space solid consists of a Surfaces list, which provides detailed definitions of half-spaces by the surfaces. It has a Boundary box defined by an upper and a lower point, a Volume, an Area, and a reference to the Polyhedron. The names of the solids and the polyhedrons must be uniquely given to be used as a reference for other blocks.

A Nested Boolean Solid type has been introduced in the GDML schema. Boolean operators like union, intersection, and subtraction can be used in a nested structure. The description of this solid will be interpreted into G4BooleanSolid in the extended Geant4 GDML parser. In the end, the Geant4 GDML parser has been extended to process the half-space solid information in the extended GDML file and construct a model from it.

The conversion of CAD geometry half-space solid has been achieved using McCad (Wu, 1987-1992), which is an open-source MC geometry conversion program based on the Open CASCADE (OCC) library (Open CASCADE Technology, 2015). As the key to the conversion process, an automatic decomposition function that generates a set of splitting surfaces and decomposes the CAD solid using OCC Boolean operations has been implemented. Algorithms have been optimized to detect and sort splitting surfaces, as well as introduce new assistant splitting surfaces.

The conversion process starts with solids. A BRep CAD solid is composed of boundaries, which are also called Faces. The analytic representation of a face is called a Surface. All the surfaces of a solid must be checked in McCad to determine whether they are splitting surfaces. A novel algorithm has been developed in McCad to detect splitting surfaces. Using this algorithm, a BRep solid is meshed into triangle facets employing the OCC library function BRepMesh_ IncrementalMesh. A surface is a splitting surface if it has collisions with at least one facet of other faces or at least two facets of other faces located on different sides of this surface. Taking the two solids in Figure 4 for illustration, Surface-A collides with Facet-B of Face-B; thus, Surface-A is a splitting surface. Facet-E and Facet-F are located on different sides of Surface-C; hence, Surface-C is also a splitting surface.

In some special cases, the boundary surfaces are not suitable for splitting surfaces even if they pass the detection. For example, in Figure 5A, the cylinder Surface-A is an invalid splitting surface. If cutting the solid by this cylinder, the resultant solid outside the cylinder is not a half-space solid, but a “ghost” solid. An algorithm is provided in McCad to introduce an additional splitting surface. The solid in Figure 5A is cut with the assistant surface Surface-B and Surface-C, and the solid is split into three solids with regular shapes. In addition, optimizing the splitting order can significantly reduce the resultant solids. For example, in Figure 5B, using Surface-C to cut the solid would produce fewer solids than Surface-D. Computational time is reduced in the conversion process by fewer cutting operations. Computational time is reduced in a particle transportation simulation due to fewer solids in the model. More detail about the implementations of these algorithms in McCad is given in (Lu et al., 2017).

To make the creation of half-space solids possible for Geant4, a new GDML interface has been developed in McCad to export GDML files. Material compositions are defined by their names and densities and are assigned to a group of CAD solids as attributes. The linking information of a material and a solid is provided in the Structure block of a GDML file. The half-space solids generated from decomposing a CAD solid are exported in the Solids block, and they are united as one nested Boolean solid to represent the original CAD solid. The polyhedrons are produced from the CAD solids for visualization purposes. As a reference, the names of the polyhedrons, materials, and solids are presented uniquely in the GDML file.

As a Geant4 tessellated solid is produced from faceting the CAD solids, it is clear that they can be generated using the OCC function BRepMesh_IncrementalMesh. Therefore, McCad has been implemented with a conversion function for G4TessellatedSolid (Qiu et al., 2016) so that a hybrid model can be produced using together half-space solids and tessellated solids.

As shown in Figure 6, the test comparisons cover all the Geant4 primitive solids supported by the half-space solid. The volume size, the relative position, the safety distance, and the intersection distance are computed, and the results of the primitives are taken as references. Each function is evaluated by testing 10⁶ samples. The points and rays are randomly sampled inside the bounding box with a margin of factor 0.5 in each direction.

The test of the function Inside is considered passed if all the sample points yield consistent results from the two solids. The comparison results of safety and intersection distances are presented by the maximum absolute difference among all the samples. The results are shown in Table 1. Note that the methods for calculating the safety distance in some geometry shapes are different between half-space solids and Geant4 primitives. In this case, the comparisons are ignored, and the results are presented as N/A. The CPU time used for testing all the functions is summed and compared, and the value is presented in Table 1 by ratios taking Geant4 primitives as references.

The tests of function Inside are passed in all the shapes. For the results of safety and intersection calculations, the differences between half-space solids and primitive solids are mostly below 10⁻⁹ mm, which is at the same level as the geometry tolerance. The only exceptional shape is the trapezoid. The reason could be that the trapezoid is defined as a Geant4 primitive G4Trd by the locations of its eight points, while it is defined as a half-space solid by the surface parameters. Due to the limit of the precision, the sloped plane of the trapezoid would not be identical between the two representations, which results in slight differences in the results. Nevertheless, the differences are negligibly small. In general, more CPU time is used for half-space solids than for the Geant4 primitives, which means further optimization of the codes and algorithms is needed for G4HalfSpaceSolid.

The use of half-space solid type in the modeling of complex geometries is tested using the ITER Benchmark model (Wilson et al., 2008), which is a CAD model of the International Thermonuclear Experimental Reactor (ITER) facility. This CAD model has more than 900 solids, which consist of cylinder, cone, quadric, and torus surfaces. It is converted using McCad and exported to a GDML file. Approximately 3,000 half-space solids are produced in this step. In order to verify this model, a tessellated solid model has been generated using McCad from the same CAD model as half-space solids. Tessellated solids are generated one-to-one with the half-space solid on the same CAD solid so that direct comparison on each solid is allowed, for example, on volume. The precision of the tessellated solid is controlled by a parameter called deflection, which is the relative tolerance of an edge/facet to the original curve surface defined for the tessellation process. The deflection used in this tessellated model is 10⁻³. These two models are shown in Figure 7. A source has been set up on a sphere surface covering the entire model, and the Geantino particles are generated randomly on the surface. The direction of the particle is identical to the surface normal, pointing to the inside of the sphere. The Geantino particle is a virtual neutral particle that sees all solids as transparent and has no collision during particle transport. In total 1E7 event is simulated, and the track lengths are scored in all solids using the G4PSTrackLength scorer, and the CPU time has also been recorded.

Due to the large number of solids, the number of facets in the tessellated solid model has exceeded the GDML limitation. To avoid this problem, the tessellated solid model was separated into two models by dividing the solids into two halves. Using the identical source definition, the total Geantino track length in the same spatial region must be identical, regardless of the geometry model. Therefore, this treatment does not affect the track lengths scored on the solids, assuming the tessellated solids are identically defined in the separated model as the full model. However, the simulation time summed up from the two separated models will be underestimated compared to the full model because fewer solids are involved in the particle navigation.

The computational time is shown in Table 2, using one 3.40 GHz processor. The computational performance of the half-space solid model is at least 30% faster than that of the tessellated solid model. This performance could possibly be further improved, for example, by using the well-optimized G4MultiUnion to unite the half-space solids instead of using the nested Boolean solid. In addition, a certain number of particles are lost due to the geometry errors in the high-complexity model. The lost particles in 10⁷ samples are given in Table 2 as well, which shows that the half-space solid model has fewer geometry errors. Note that the lost particles reflect the qualities of the McCad-generated models on complex geometries using two different approaches and do not imply or reflect issues of the G4HalfSpaceSolid and G4TessellatedSolid implementations. The comparison of the track lengths is presented in Figure 8 as ratios, which take the results from the tessellated solid model as references. The ratios are, in general, within ±0.5%. A few solids have ratios of up to 2%. These larger variations are probably due to the geometry difference between the two representations because tessellated solids are facet approximations of CAD solids. The relative error of the half-space solid calculation is provided in Figure 8 as well, with a similar trend found in the calculation of the tessellated solid model. It should be noted that the relative error is higher than the ratio of the track lengths, which is unexpected. The reason why track lengths comparison shows a much lower deviation is likely because the same random number generator and seed are used for the simulation of two models, resulting in source particles generated in the same directions in the two runs.