Modeling quadrupole mass analyzers with OmniTrak

Linear quadrupoles are widely used for ion mass spectrometry. They can separate ions by the ratio of mass number to charge number (M/Z) without a magnetic field, leading to compact, light and inexpensive systems. An assembly consists of a source of low-energy ions that enter an extended electrostatic  quadrupole. The voltage applied to the lens has both DC and AC components. The orbits of most ions are unstable. Depending on the choice of voltages and RF frequency, only ions within a small range of M/Z have stable orbits and pass to a downstream detector. This report demonstrates the application of the three-dimensional OmniTrak code to characterize quadrupole mass analyzers.

Cross-section of an electrostatic quadrupole.

Figure 1. Cross-section of an electrostatic quadrupole. a}) Lines of constant potential. b}) Electric field lines.

Assume the quadrupole has ideal hyperbolic electrodes with an electrode-tip radius r0 (i.e., the distance from the axis to the closest point on an electrode). The beam propagates in the z direction. Figure 1 shows a cross section of the geometry in the x-y plane. The voltage applied to the top and bottom electrodes has the form of Eq. 1 where U is the DC component and V is the magnitude of the RF voltage with angular frequency Ω. Voltages with opposite polarity are applied to the y and x electrodes.


The electrostatic potential in the gap has the form of Eq. 2.  The y component of electric field is given by Eq. 3.  Equation 4 describes ion motion in the y direction. With the introduction of the dimensionless variables of Eqs. 5, 6 and 7, the equation of motion takes the form of Eq. 8. Motion in the x direction is described by Eq. 9.

Equations 8 and 9 must be solved numerically. There are regions in (a,q) space where ion motion in x or y may be stable or unstable. Here, the term stable implies that the trajectories are bounded oscillations. In unstable regions, the displacement increases exponentially so that in a short time ions strike the electrodes. In particular, ions have stable trajectories in both the x and y directions within the shaded region shown in Figure 2. (Note that there is also a symmetric stability region below the q axis. Changing a to -a simply exchanges the roles of x and y).

Region of a-q space with stable motion in both x and y

Figure 2. Region of a-q space with stable motion in both x and y.

For reference, note that the ratio of the parameters is proportional to the DC voltage component divided by the AC component (Eq. 10). A fixed ratio corresponds to a straight line in (q,a) space like the dashed line in Fig. 2. Motion along the line starting at the origin is equivalent to raising the magnitude of U and V while keeping a constant ration. For a section of the scan, the curve in Fig. 2 passes through the stability region. The tip of the stability region has the coordinates q = 0.706 and a = 0.237.

The operation of the quadrupole as a mass filter is clear if we plot regions of stability in terms of voltage values U and V rather than the dimensionless parameters a and q. Figure 3 shows stability regions for two values of M/Z. The M/Z value for the green shaded region is half that of the tan region. The working line follows a fixed voltage ratio, U/V < (0.237)/(2 × 0.706). Note that as the magnitude of the applied voltages increases, the working line first passes through the green stability zone and then through the tan zone. A plot of the detector response versus a scan of voltage magnitude will show two peaks corresponding to the two values of (M/Z). Higher voltages correspond to higher mass-to-charge ratios. Absolute values of (M/Z) may be inferred from known values of the voltage, quadrupole geometry and RF frequency. Note that the ratio U/V = 0.168. gives ideal resolution but no throughput. Lowering the DC voltage allows more particles to pass to the detector at the expense of reduced mass resolution. With no DC offset (U = 0), all ions pass through the quadrupole.

Stability regions and a working line for ion transport as a function of the applied voltages U and V

Figure 3. Stability regions and a working line for ion transport as a function of the applied voltages U and V. The shaded regions show the stability zone for two ion species. The green region has a mass to charge ratio (M/Z) equal to half that of the tan region.

We now have the theoretical resources to design a practical system. Take the quadrupole length as L = 15.0 cm. Fabricating hyperbolic electrodes is expensive, so we’ll use stock electrodes with a circular cross section to approximate the field. The choice is acceptable as long as the ions are confined close to the axis. The electrode diameter is d = 1.0 cm. The electrodes are displaced 0.9 cm from the propagation axis, giving the minimum distance from the axis to an electrode of rc = 0.4 cm. We can identify two questions that can be answered by a 3D code:

  • How does the choice of rc versus d affect the linearity of electric fields near the axis?
  • What is the value of ro for equivalent hyperbolic electrodes? In other words, if we used hyperbolic instead of circular electrodes with a given voltage, what choice of ro would give the same electric field variation near the axis?

We’ll concentrate on the second issue, calculating the electric field with HiPhi. The quadrupole assembly is centered in a solution volume that covers the range -3.0 cm ≤ x,y ≤ 3.0 cm and 0.0 cm ≤ z ≤ 20.0 cm. The wall of the solution volume is grounded. The element size in the quadrupole region is 0.05 cm. Figure 4 shows the resulting mesh.

Conformal representation of the quadrupole electrodes

Figure 4. Conformal representation of the quadrupole electrodes.

I used normalized voltages of +1.0 V on the y electrodes and -1.0 V on the x electrodes in the electrostatic solution. Figure 5 shows calculated values of Ex in a scan through the axis along x at the quadrupole center. Over the range -0.2 cm ≤ x ≤ 0.2 cm, the field varied linearly with x to within 0.01%. The calculated dependence is Ex(x) = (1.255 × 10^5)x V/m. A comparison with Eq. 3 for hyperbolic electrodes implies the relationship of Eq. 11. The equivalent hyperbolic pole tip distance is 0.3993 cm, quite close to the physical distance rc = 0.4000 cm. My initial guess of rc  = 0.4 d turned out to give excellent results.

Scan of Ex(x)

Figure 5. Scan of Ex(x) at position y = 0.0 cm, z = 10.0 cm for normalized electrode voltage magnitude of 1.0 V.

The electrostatic solution was used as a basis for OmniTrak calculations of ion trajectories. I concentrated on singly-charged carbon ions (M = 12, Z = 1) with an initial kinetic energy of 40 eV. The drift velocity was 2.5 × 10^4 m/s, so the transit time through the spectrometer was about 6 μs. There should be several RF periods over the transit time to ensure that ions experience an average focusing force. I picked f = 2.5 MHz, corresponding to Ω = 1.571 × 10^7 1/s. A time step choice of Δt = 10.0 ns gave an acceptable value of 40 integration steps per RF period. I picked the PList command to initiate ion trajectories near the z = 0.0 cm boundary with velocity directed along z. To sample positions in x-y space, I used a spreadsheet to create a circular distribution of 36 particles at radius 0.1 cm. (There was no point in representing an on-axis particle. It would traverse the system with no transverse motion in all circumstances). Equations 6 and 7 can be used to find the U and V values for a working line that passes through the tip of the stability region (Fig. 2). Taking a = 0.237, q = 0.706 and substituting values for the quadrupole leads to the Eqs. 12 and 13.

The values for carbon ions are U = 14.65 V and V = 87.28 V. The value of U must be lower to observe transmitted ions.

The full OmniTrak input script is listed at the end of the report. The critical commands to represent the RF quadrupole are:

EFIELD3D: QuadAnalyzer.HOU 0.9875
MODFUNC E > 14.0 + 87.28*sin(1.571E7*t)

The first command loads the normalized electric field. The adjustment factor (β = 0.9875) determines the voltage magnitude and hence the position along the working line. We can make small changes to study the transmission width of an ion type, or large changes to study a different type. For example, the adjustment factor for singly-ionized iron is about β = (56/12) = 4.67. The ModFunc command serves two functions:

  • Add a time variation of fields during the ion transit.
  • Set the slope of the working line (i.e., determine the acceptance of the system).

With regard to the first function, all values of potential in the electrostatic solution are multiplied by the current value of the algebraic function. This process is valid when the propagation distance of an electromagnetic signal over an RF period is much larger than the system.

My procedure was to set a value of U less than 14.65 V in the ModFunc command and then to do a series of runs with varying β to change the total voltage level. Figure 6 shows the results for U = 12.0 V and 14.0 V. As expected, the peak signal occurred near the tip of the stability region and a lower value of U gave an increased transmission width. Figure 7 shows the effect of small changes of β on ion trajectories.

Particle transmission factor

Figure 6. Particle transmission factor as a function of the relative voltage magnitude and U.

Projected carbon ion trajectories

Figure 7. Projected carbon ion trajectories in the z-y plane for V = 87.28 V and U = 14.0 V. A small change in the applied voltage magnitude causes transmission to drop from 100% to 0%. Note that the y scale has been expanded for clarity.

Finally, note that the calculation procedure can be automated for extended data analyses in the background. For example, to make a scan in voltage magnitude, the EField3D command could be changed to

EFIELD3D: QuadAnalyzer.HOU %1

where %1 represents a command line parameter. Multiple OmniTrak calculations could be controlled with a batch file with a form like this:

REM Xenos batch file, Field Precision
START /B /WAIT C:\fieldp_pro\xenos/omnitrak.exe QuadAnalyzer 0.75
START /B /WAIT C:\fieldp_pro\xenos/omnitrak.exe QuadAnalyzer 0.80
START /B /WAIT C:\fieldp_pro\xenos/omnitrak.exe QuadAnalyzer 0.85
START /B /WAIT C:\fieldp_pro\xenos/omnitrak.exe QuadAnalyzer 0.90

The end result is a set of particle files that may be analyzed with a text editor or a user script to find the fraction of particles that reach z = 20.0 cm.

Example: OmniTrak input script

  EFIELD3D: QuadAnalyzer.HOU 0.9782
  DUNIT:   1.0000E+02
  MODFUNC E > 14.0 + 87.28*sin(1.571E7*t)
  DT  1.0E-8
    12.0  1.0  40.0  0.0996   0.0087 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0966   0.0259 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0906   0.0423 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0819   0.0574 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0707   0.0707 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0574   0.0819 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0574  -0.0819 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0707  -0.0707 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0819  -0.0574 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0906  -0.0423 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0966  -0.0259 0.01 0.00 0.00 1.00
    12.0  1.0  40.0  0.0996  -0.0087 0.01 0.00 0.00 1.00
  PARTFILE: QuadAnalyzer


[1] PDF version of this report:

[2] Information on HiPhi:

[3] Information on OmniTrak:

[4] Contact us :

[5] Field Precision home page:

Building solids in MetaMesh

The programs Geometer and MetaMesh are used to build three-dimensional meshes for our multi-physics solution programs using constructive solid geometry. The user has the option to assemble objects from native models (e.g., spheres, cylinders, boxes, turnings,…) or to import complex shapes from SolidWorks and other 3D CAD programs. This article reviews how to make objects with conformal surfaces from native models. In particular, I want to clarify some points about conformal fitting.

To start, it’s important to understand the concept of the initial foundation mesh versus the final conformal mesh. The foundation mesh is a regular, structured set of box elements that fill the solution volume. Regions in the solution volume are collections in elements that represent physical objects (such as electrodes and dielectrics in electrostatic solutions). In the volume-fitting process, the goal is to make an assignment of region numbers to the foundation-mesh elements for the best representation of the physical objects. The process of surface fitting occurs after the foundation mesh is complete. Here, elements on the object surfaces are shaped to follow the mathematical definitions of native parts or CAD models to produce the final conformal mesh. As a result, elements near the surface become generalized hexahedrons.

The two mesh-generation tools serve the following functions:

1) Geometer is an interactive, graphical environment to create or to edit parts in the foundation mesh. A part is a single native or CAD model. Objects may be composed of multiple parts. The output of Geometer is a text script that lists specifications of the foundation mesh and the included parts. The script (FName.MIN) acts as input to MetaMesh. A script may also be prepared directly with a text editor.

2) MetaMesh analyzes the script to produce a binary output file of the coordinates of mesh nodes and element identities used by the solution programs. If the user makes no changes to the script produced by Geometer, the MetaMesh output is a regular mesh of box elements, identical to the foundation mesh. The user adds or edits commands of the MetaMesh script to control conformal fitting. This process is the focus of this report.

Example geometry

Figure 1. Example geometry.

Figure 1 shows the test geometry, a two-electrode gap. The solution volume covers the space -3.0 ≤ x,y ≤ 3.0, 0.0 ≤ z ≤ 7.0. The upper electrode has a spherical depression. The lower electrode has a spherical extension and a circular extraction port. The system has cylindrical symmetry about the z axis. A quick approach is to represent the upper and lower electrodes as turnings that extend beyond the boundaries of the solution volume. Instead, we will work with basic shapes to illustrate the construction procedure.

If we start a new script in Geometer with the given solution-volume dimensions, the program automatically adds a box region that fills the volume with the name SolutionVolume. Using the command to edit region names, we change the name of the first region to Air and the second and third region names to Upper and Lower. We then add a part with the name UpperElectrode and the region Upper. The part is a Box with dimensions Lx = Ly = 6.0 and Lz = 2.0 with a displacement 6.0 in z. To cut the depression, we add a sphere with name UpperCut, region Air, radius 2.5 and displacement 2.5 in z. The action changes the region identity of overlapped elements and nodes in the upper electrode to Air. Next we add a Box with name LowerElectrode, region Lower, dimensions Lx = Ly = 3.0 and Lz = 2.0 and displacement 1.0 in z to represent the lower electrode. The extension is a Sphere with name LowerExtension, region Lower, radius 2.0 and displacement z = 1.0. Finally, we add a Cylinder with name Aperture, region Air, radius 0.5, and length 7.0 to cut the hole. Figure 2 shows the resulting mesh created by MetaMesh if no changes are made to the script. It shows the familiar stair-step pattern for curved surfaces characteristic of regular meshes.

Mesh with no conformal fitting

Figure 2. Mesh with no conformal fitting.

Before proceeding, we should note several significant points.

1) Some of the parts (like the Aperture) extended out of the solution volume. This is not a problem — external sections of parts are ignored.

2) Regions may be created by the actions of multiple parts. For example, the lower electrode was the sum of the LowerElectrode and LowerExtension parts minus the Aperture part.

3) Note how the superposition principle affects the ordering of parts. For example, we enter the UpperElectrode and then cut then spherical depression by overwriting part of the object with Air elements. Only then do we then enter the LowerElectrode. If we entered both electrodes followed by the air sphere, both electrodes would have cuts.

4) Finally, there are special considerations for cutting holes when the mesh will be used in a solution program that depends on the node identity.

MetaMesh input script with fitting commands

Figure 3. MetaMesh input script with fitting commands.

Figure 3 shows the input script produced by Geometer. The highlighted lines have been added to implement conformal fitting. The Global section (Lines 1-14) performs three functions: a) set the dimensions of the solution volume, b) specify the sizes of elements in the foundation mesh and c) define names for the regions. Part 1 (15-21) is the solution volume region defined by geometer with all elements set to region number 1 (Air). Part 2 (22-28) is the upper electrode that overwrites the air elements. Part 3 (29-37) is the spherical concavity in the upper electrode. In Part 3, the command

Surface Region Upper

states that elements of the part in contact with elements of type Upper should be shaped so that they conform to the spherical surface. The command

Coat Upper Upper

states that nodes shared with elements of region Upper should be marked as Upper. This operation ensures that the nodes of all elements in the upper electrode are marked as part of the fixed potential region. (Note: a valid solution is obtained even if we forget to add this command. Solution programs like HiPhi automatically set all nodes connected to fixed-potential electrodes to the fixed-potential region number.)

Moving on, Part 4 (38-44) is the lower electrode while Part 5 (45-52) is the spherical extension. The command

Surface Part UpperCut

states that elements adjacent to those marked as part UpperCut are shaped to lie on the spherical surface. The final part is the cylindrical aperture. The length of 7.0 ensures that the cylinder extends above the spherical extension while the section outside the solution volume is ignored. The commands

Surface Region Lower
Coat Lower Lower

move nodes shared with elements of the electrode to conform to the cylindrical shape and correct the region assignment of these nodes. Figure 4 shows a detail with nodes of the cylinder area.

Detail of the mesh around the Aperture

Figure 4. Detail of the mesh around the Aperture showing node assignments.

To conclude, I’ll clarify why a Surface Part command is used in Part 4 rather than a Surface Region command. Fitting occurs after volume assignment operations (i.e., the foundation mesh has been completed). In this case, there are elements of with region identity Air in the aperture as well as in the gap between electrodes. Clearly, we do not want nodes inside the aperture to be shifted to the spherical surface of the extension. Such moves would cause severe distortion of the mesh. Therefore, we use a more specific designation to narrow down the set of nodes that should be shifted.


[1] Information on MetaMesh:

[2] Information on Geometer:

[3] Contact us :

[4] Field Precision home page:


OmniTrak: superposition of magnetic field tables

OmniTrak is our software package for charged-particle guns and beam dynamics. The solution program can read output files from HiPhi and Magnum for 3D electric and magnetic field information. These complex files record values of potential on conformal meshes. Tables are another option for field information. Tables record [Ex,Ey,Ez] or [Bx,By,Bz] on a regular box mesh. Table files are easy to create and facilitate import of field information from other programs. Sometimes, it is useful to use an intermediate table instead of directly reading information from HiPhi and Magnum. In this case, the user uses the MATRIX command in PhiView or MagView to convert the values of potential on the conformal mesh to a table of electric field or magnetic flux density values. This reduces the amount of work performed by OmniTrak and can speed up many calculations.

A user recently suggested another application of tables that is helpful when the fields of a magnetic transport element can be expressed as a linear superposition of the fields produced by different control coils. For example, it may be useful to model the performance of a magnetic scanning system by calculating orbits for several values of current in the the X and Y drive coils. A quick way to do the calculation would be to prepare tables for the fields of the X and Y drive coils with a normalized current, and then to load a superposition of the fields into OmniTrak with different scaling factors. The process could be controlled by a batch file for a parameter search.

Recently, we make changes to OmniTrak to enable this procedure. The command to load a table has the form

BTable3D COIL01.MTX 1.5

where COIL01.MTX is the name of the table in the working directory and 1.5 represents an optional scaling factor. The factor multiples values of [Bx,By,Bz] — the number may be positive or negative. In the past, a second BTable3D command would simply cause the previous table to be overwritten. In the current version of OmniTrak, multiple BTable3D statements result in a superposition of values. There is no limit to number of commands — here is an example with two:

BTable3D COIL01.MTX 1.5
BTable3D COIL02.MTX 2.0

OmniTrak reads the first table in the normal way with scaling factor 1.5. When the second command is encountered, the program opens the file and checks that the map parameters (IMax, JMax, KMax, Xmin, Xmax, Ymin,…) are identical. If so, OmniTrak reads the second table and adds values of {Bx,By,Bz] to those already in memory with the scaling factor 2.0. In summary:

  1. There is no limit to the number of BTable3D commands.
  2. Multiple maps must use the same mesh.
  3. BTable3D commands may appear in any order in the FIELDS section of the OmniTrak script.
  4. Shifts and rotations are applied after all tables have been superimposed.
  5. The user must ensure that the superposition is physically valid.
  6. The ETable3D command has also been extended.

Figure 1 shows data from a test example: two identical shielded solenoid lenses at different positions along the z axis. The blue squares represent the tabular field of the first coil and the red squares show the field of the second downstream coil. The line is the total field calculated with the OmniTrak BSCAN command resulting from the command set listed above. The following parameters were used in the Magview MATRIX command to create the tables: IMax = 21, XMin: -10.0, XMax = 10.0, JMax = 21, YMin = -10.0, YMax = 10.0, KMax = 240, Zmin: -12.0 and ZMax = 12.0.

Test case - superposition of two solenoid lenses

Figure 1. Test case – superposition of two solenoid lenses.


[1] Contact us :

[2] Field Precision home page:

GamBet shielding calculation walkthrough

In a previous article, I discussed the basis of the H*(10) standard for radiation safety. Here, I’ll continue by demonstrating techniques for shielding calculations. I’ll use a published simple benchmark calculation to compare the results to several other Monte Carlo codes and to illustrate some of the advantages of GamBet. For this, I am grateful to the people at RayXpert who have made reports on a number of benchmark calculations available on the Internet.  In particular, I admire their patience and fortitude in learning four other codes besides their own to make the comparisons.

Geometry of the benchmark calculation

Figure 1. Geometry of the benchmark calculation taken from the RayXpert report.

Figure 1 shows the geometry taken from the report H*(10) dose rate evolution with some different shielding. Targets are located 100 cm from a Co60 point source with activity 10^9 Bq[1]. The dose at a depth of 1.0 cm was calculated for a bare target and with three shield configurations. The setup illustrates two points I made in the last article about H*(10) calculations:

  1. The results are not sensitive to the detector geometry. In this case, a 4.0 cm radius sphere of tissue equivalent material is used rather than the ICRU 15.0 cm radius sphere.
  2. The radiation weighting factor (Wr) is taken as unity for the photons and secondary electrons and there is no organ correction for the phantom (Wt = 1.0). In this case, the biological dose (Sv) equals the physical dose (Gy).

I limited my calculations to the bare target and to one behind a 2 cm lead shield. Here are the results reported by RayXpert. The values are given in units of μGy/hour.

             Dosimex-G  Mercurad  Microshield MCNPX RayXpert
No shielding    349        353         347      355   365
2 cm lead       143        152         141      146   149

The targets subtend only a small fraction of the solid angle surrounding the source, so it would be foolish to use an isotropic source. Instead, this is an ideal opportunity to use the GamBet Radiation Source Tool, illustrated in Fig. 2. After the user enters values in the dialog fields, the program creates a standard source file. The idea is to fill a circle of a specified radius in the direction of interest with a large, random distribution of particles. The flux assigned to each particle is calculated to represent the specified source activity. In this case, I created a distribution of 200,000 photons of energy 1.17 MeV and an equal number at 1.33 MeV. The distribution filled a circle of radius 4.0 cm at a distance of 92.0 from the source, ensuring complete irradiation of the detector. The energy deposition was symmetric about the line from the source through the detector, so I further reduced run time by using the option for cylindrical scoring in GamBet. The air region was represented by the standard Penelope dry air model and the tissue-equivalent target was represented by the custom material definition discussed in the previous article:

 Name TissueEquivalent
 Component O 1.0000
 Component C 0.1940
 Component H 2.1032
 Component N 0.0390
 Density 1.00
GamBet Radiation Source Tool dialog

Figure 2. GamBet Radioactive Source Tool dialog.

Figure 3 shows the dose distribution calculated by GamBet for the bare detector with photons entering from the left[2]. Even for this simple case, there is quite a lot going on. The dose at the surface is lower than the peak dose in the material. This is the result of the buildup of the secondary-electron distribution. The energy loss-rate for electrons is much higher than that of photons of the same energy. The depth to reach the equilibrium dose is about 0.6 cm, hence the choice of the H*(10) measurement depth as 1.0 cm for γ rays in the MeV range. The dose rate decreases gradually through the material because of γ ray attenuation. Note the difference in the air dose rate at the entrance and exit of the target. The enhanced downstream dose rate results from electrons driven out of the material. The H*(10) dose rate determined by GamBet was 338 μGy/hour, in good absolute agreement with the values reported by RayXpert.

Dose-rate disttribution for a bare target

Figure 3. Dose-rate distribution for a bare target.

Figure 4 shows the dose-rate distribution with 2.0 cm lead shield placed a short distance in front of the target. The shield extends to the outer radius — the portion irradiated by the 4.0 cm radius beam is clearly visible. The dose buildup depth in lead is shorter than the element size, so it is not visible[3]. Also, note the absence of a buildup region in the target. Most of the equilibrium secondary electron distribution from the lead crosses the small air gap and enters the target. The GamBet H*(10) prediction is 146 μGy/hour, again comparable to the values determined by the other codes. Finally, Figure 5 shows dose-rate scans along the diameter of the spherical target for the two cases. The dashed line is the H*(10) point.

Dose-rate distribution for a 2.0 cm lead shield

Figure 4. Dose-rate distribution for a 2.0 cm lead shield.

Dose-rate scans along a diameter of the spherical target

Figure 5. Dose-rate scans along a diameter of the spherical target.

One conclusion of this study is that all radiation-shielding codes give reasonable results. This is not surprising because the physics and mathematics base for Monte Carlo radiation codes is over half a century old. The question within this constraint is: why choose GamBet? Beyond its low price, GamBet has several advantages, some of which are apparent in this study:

  • Zoning in GamBet is performed with conformal finite-element meshes, allowing plots such as those of Figure 3 and 4. Complex processes may enter even the simplest radiation calculations, so detailed pictures of dose distributions are essential for evaluating results.
  • GamBet is fast for two reasons: 1) the code is designed to encourage division of solutions into stages and 2) the code supports efficient parallel processing. For this example, the run time with 1.6 million incident model photons (4 parallel processes with 400,000 photons each) was 7 minutes.
  • The element approach in GamBet is more effective than combinatorial solid geometry to represent complex systems. The advantage is not apparent in this calculation, but is overwhelming when representing something like a human body phantom.
  • Although all Monte Carlo physics engines are about equal for MeV γ rays, Penelope has far more detailed support for low-energy X-rays.
  • GamBet has several advanced features, such as integration of calculated 3D electric and magnetic fields and direct coupling to electron beam and thermal transport programs.


[1] This article reviews the decay scheme of Co60: Modeling radioactive sources with GamBet.

[2] GamBet uses standard units of Gy/s. The conversion factor to μGy/hour is 3.6E9.

[3] The element size in GamBet determines the display resolution and the accuracy for fitting shapes, but has no effect on the physics of the radiation interactions.

[4] Contact us :

[5] Field Precision home page:

H*(10) calculations with GamBet

I recently had an inquiry from a prospective customer. He had investigated several Monte Carlo options and was impressed with one package because the developer emphasized that the software could perform H*(10) calculations. In reality, every Monte Carlo X-ray code has this capability. The only requirement is that the user be aware of the definition. To remove any doubts, I will review the significance and technical background of H*(10) calculations in this article.

In applications to personnel shielding, it is not sufficient simply to know the radiation field at a location. By radiation field, I mean the fluxes of energetic electrons and photons along with their energy spectra. The critical issue is how do the particles interact with tissue and how does the dose (deposited energy divided by mass) build up inside the organism. A further complication is that the same dose may have different biological effects, depending on the radiation and tissue type. The standard unit of physical dose is the gray (Gy), equal to 1 joule/kg of deposited energy. The unit of biological dose is the sievert (Sv) which also has units of J/kg. The difference is that the biological dose includes weighting factors that depend on the radiation type (Wr) and the tissue type (Wt) to indicate the relative potential for biological harm. Dose values are related by

D(Sv) = D(Gy)*Wr*Wt

The radiation weighting factor is Wr = 1.0 for photons and electrons, so we need not worry about it in GamBet calculations. The tissue weighting factor Wt is difficult to estimate, and therefore there is considerable disagreement about the values. For whole-body irradiation, the convention is to take Wt = 1.0. The implication is that for most shielding calculations at energies of interest for medical applications, sieverts are equivalent to grays.

The critical issue is the dose buildup through interactions with tissue. To standardize measurements and calculations, the ICRU (International Commission on Radiation Units and Measurements) defined the Directional Dose Equivalent, H*(d). The quantity applies to radiation moving predominantly in one direction, such as X-rays outside a radiation shield. A different quantity, Hp (the personnel dose equivalent), applies to approximately isotropic radiation fields, like the time that Mr. Spock was trapped in the Propulsion Chamber.

The ideal calculation or measurement defined by the ICRU uses a 30-cm-diameter sphere of tissue-equivalent plastic with a density of 1 g/cm^3 and a mass composition of 76.2%  oxygen,  11.1% carbon, 10.1% hydrogen and 2.6% nitrogen. The quantity H*(d) is the biological dose at a depth d below the surface of the sphere in the direction of the radiation (Figure 1). The most common choice of d is 10.0 mm, hence the term H*(10). The reasoning is as follows. When a flux of gamma rays enters a material, the dose increases moving away from the surface because of the generation of secondary electrons (e.g., Compton electrons). The electrons deposit energy more rapidly than photons. The secondary electron density grows with distance until the production rate equals the absorption rate. At this equilibrium point, the dose reaches a maximum. For gamma rays in the 1 MeV range (typical of radioactive sources) in a material with density 1 gm/cm^3 (e.g., water), the equilibrium depth is a 10.0 mm or less. The following article describes benchmark calculations that illustrate this effect.

Geometry to define the Directional Dose Equivalent

Figure 1. Geometry to define the Directional Dose Equivalent.

To conclude, here are a couple practical suggestions for H*(10) calculations in GamBet. First, in defining custom materials, GamBet follows the Penelope convention using the stochiometric composition (the relative number of atoms in an equivalent molecule) rather than mass fractions. To make the conversion, assume that F0, FC, FH and FN are the stochiometric fractions for oxygen, carbon, hydrogen and nitrogen. Because the quantities are relative, we take F0 = 1.00. Using the atomic weights and mass fraction of 0.762 for oxygen, we can determine the other quantities from these equations:

FC*12.011/15.994 = 0.111/0.762
FH*1.00794/15.994 = 0.101/0.762
FN*14.00674/15.994 = 0.026/0.762

Here’s the result as it would appear in the Composition section of a GamBet script:

  Name TissueEquivalent
  Component O 1.0000
  Component C 0.1940
  Component H 2.1032
  Component N 0.0390
  Density 1.00

Finally, note it is not necessary to use a 30 cm sphere of the tissue equivalent material in calculations. Any block of material will give about the same answer, as long as the object has depth and width greater than d.

The next article describes a walkthrough example that gives some good insights into the physics of the H*(10) calculation. The example illustrates several GamBet techniques and new features and provides an opportunity to make comparisons with several other radiation-shielding codes.


[1] Contact us :

[2] Field Precision home page: