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.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

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:

Material
 Name TissueEquivalent
 Component O 1.0000
 Component C 0.1940
 Component H 2.1032
 Component N 0.0390
 Density 1.00
 Insulator
End
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.

Footnotes

[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 : techinfo@fieldp.com.

[5] Field Precision home page: www.fieldp.com.

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:

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

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.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.

Variance reduction techniques in GamBet

The GamBet software suite calculates the interactions of energetic photons, electrons and positrons with matter. In a previous sequence of three articles (Monte Carlo methods versus moment equations: Part A, Part B and Part C), I discussed how Monte Carlo calculations approximate the behavior of large collections of particles. The essence of the Monte Carlo method is to follow the interactions of a finite set of model particles with the assumption that they represent the average behavior of the large set. The advantage is that there is a close connection with the interaction physics, so that a Monte Carlo calculation is often called a simulation. The disadvantage is that it is a statistical process, so the accuracy of predictions is limited by random variations set by the number of model particles or events, N. The accuracy in estimating an average quantity is usually expressed as <x> ± σ, where <x> is the ideal average and σ is the standard deviation set by the finite number of samples. The standard deviation equals the square root of the variance, the average of the squared differences from the mean value. The standard deviation scales as 1/√N. Improving the accuracy of estimates by a factor of 10 requires 100 times the number of model particles. The implication is that high accuracy leads to long run times.

The term variance reduction[1] applies to methods to optimize Monte Carlo calculations to gain an edge on the 1/√N limit. Variance reduction is not a formalized mathematical method, but rather a set of common-sense fixes that can go a long way toward reducing run time. The essence is finding ways not to waste time on particles that will not contribute to critical results. The success and validity of variance reduction depends strongly on user judgements. GamBet uses the Penelope package for interaction physics. Penelope has built-in features for variance reduction that are implemented in GamBet with the following commands:

ENHANCE NReg NSplit [ElecP, PhotP, PosiP, ElecS, PhotS, PosiS]

GamBet is unique among Monte Carlo codes in its use of a finite-element conformal meshes to represent the solution volume. Users can divide the space into a number of regions to represent different materials or sections of an object. The Enhance command improves statistics in a critical region (such as a detector). Particles entering the region are split into NSplit particles with statistical weight 1/NSplit. With the optional string parameters, the operation may be limited to specific particle types (primary or secondary electrons, photons and positrons).

REDUCE NReg NKill [ElecP, PhotP, PosiP, ElecS, PhotS, PosiS]

The Reduce command is the inverse of Enhance. The number of particle entering a specified non-critical region are reduced while increasing their statistical weight.

FORCE [ELEC,PHOT,POSI] [HELAS,...]  Factor [NReg]

This command instructs the program to increase the probability of low-probability interactions like bremsstrahlung emission. The statistical weight of reaction products is decreased to preserve the correct energy balance between reactions.

Beyond the Penelope techniques, GamBet is structured to help achieve short run times.  The code was designed to encourage the division of calculations into manageable segments. For example, an initial segment could address a radioactive source with shielding and collimation, while a second segment could address the interaction of forward-directed radiation with tissue. The segments are connected by an escape file which records the set of model particles that reach the boundary of the first segment. The key to variance reduction in GamBet is filtering and transforming escape files for optimal performance in a following segment. The escape distribution can be modified within a GamBet run with the following command:

ESCAPEFILTER Condition01 Condition02 ...

The conditions are strings like X>0.15, T<5.0E6,…. Particles must meet the combined conditions to be included in the escape file. Conditions may apply to spatial locations, kinetic energy and particle type. The idea is to limit particles to those that will play a role in the following segment and to limit the size of the escape file. Recently we expanded the EscapeFilter conditions to include particle direction: Ux>0.1,Uz>0.95,Ur<0.25,… Here, the quantities are the components of a unit vector pointing along the direction of the velocity. One application is to limit particles to those that are aimed toward a detector or target.

GenDist as a component in a variance reduction strategy

Figure 1. GenDist as a component in a variance reduction strategy.

The GamBet package includes GenDist, a powerful utility to create or to modify large particle distributions. GenDist can act as an additional stage between calculation segments (Figure 1) to optimize particle properties for reduced variance. The basic sequence is to read an escape file, filter or transform the particle parameters and to write a modify file to be used as the source for a following calculation segment. In the past, the operations were controlled interactively by the user in a program window. We have recently added a script capability for autonomous production runs. Here is a summary of the new script commands:

READ FPrefix.SRC
Load an escape file

WRITE Fprefix.SRC
Write a file of transformed particle parameters, applying any filter conditions that have been set.

AXIS [X,Y,Z]
Set a reference axis for evaluating transverse velocity in transformations and filters.

FILTER Condition Value
Set any number of filter conditions. The set of conditions is the same as those in the GamBet EscapeFilter command.

XFORM GENERAL XShift YShift ZShift XRotate YRotate ZRotate
Move or rotate the particles to match the coordinate system of the next computational segment.

XFORM UNIDIST Dist
XFORM NORMPLANE Pos
XFORM CLOSETOLINE HLine YLine
Move particles in ballistic orbits following their velocity vectors. The options are shifts 1) a uniform distance backward or forward, 2) to a plane normal to the current axis or 3) to their positions closest to a line parallel to the current axis. One main application is to find the effective radius of a bremsstrahlung source for X-ray imaging applications by back-projection to the target.

BEAMSECT2D NThet
BEAMSECT3D X 0 X X
The first command converts a distribution from 2D cylindrical calculation to one suitable for a 3D calculation. The second command limits 3D beam distributions to specific transverse quadrants or mirrors a beam distribution.

SCALE Fact
Change the size of a particle beam or convert spatial units to match calculation segments. For example, a source calculation may use units of mm, but units of m may be more appropriate for the detector calculation.

For more details, use these links to download the GamBet and GenDist manuals. Free updates to the new programs are available to GamBet and Xenos users.

Footnotes

[1] The term standard deviation reduction would be a better choice, but it is a less compelling phrase.

[2] Contact us : techinfo@fieldp.com.

[3] Field Precision home page: www.fieldp.com.

 

Moving thermal sources in HeatWave

The generation of X-rays via bremsstrahlung is energetically inefficient at electron beam energies of interest for medical applications. Therefore, high-power electron beams are necessary for procedures such as CT scans. Dissipation of energy to prevent target melting is a major issue. One solution is a rapidly-rotating target. Alternatively, we recently worked with a customer using a liquid jet of gallium alloy as the X-target.

HeatWave, a component of our Xenos suite for X-ray source design, performs 3D thermal conduction calculations. Energy deposition profiles determined by the GamBet Monte Carlo code can be imported as thermal source distributions. The distributions can be modulated with user-specified functions to represent pulsed or periodic beams. This capability can provide, at best, an rough approximation to a rotating target.

The HeatWave solution is performed on a stationary mesh — it would be difficult and unwieldy to introduce a moving mesh. On the other hand, it is relatively easy to move the heat source through the mesh, achieving the same result. Accordingly, we have introduced a new capability in HeatWave to generate exact moving-target simulations. This article summarizes its operation.

Power distributions from the programs Aether (microwave fields), RFE3 (time-dependent electric fields) and GamBet (Monte Carlo X-rays and electrons in matter) can be imported into HeatWave with the SourceFile command. The only restriction is that the field and particle solutions must be performed on the same mesh as that used for the thermal solution. The result is that the elements of the HeatWave solution have assigned source power densities. File sources may be used inb both static and dynamic solutions. In the dynamic mode, a time variation may be associated with the file source. The variation may be defined by a table of [t,f(t)] values using the SourceMod command. Here, f(t) is a multiplication factor. The factor may also be defined by an algebraic expression.

In the new version of HeatWave, the element power distribution is initially copied to a reference mesh variable. A time-dependent displacement vector is defined with the new commands XDisp, YDisp and/or ZDisp. As with the modulation functions, the time-dependent vector components may be defined by either a table of values (e.g., [t, x(t)]) or an algebraic function. The element power densities from the reference mesh are periodically mapped to the computational mesh using the current value of the displacement vector [x(t), y(t), z(t)]. The mapping algorithm conserves energy. Depending on the displacement, some source energy may be mapped outside the source region. As an example, heating of a rotating target could be modeled using a sawtooth function for x(t).

Figure 1 shows the results of a demonstration calculation. A 4 MeV electron beam of radius 5.0 mm with current 23.9 μA impinges on an aluminum target of thickness 10.0 mm. The top section shows the temperature distribution in the center of the target (z = 5.0 mm) with a total irradiation time of 10 s and a fixed source distribution. In contrast, the lower section shows a solution where the source moves the right 20.0 mm during the 10 s interval. Note that the peak temperature drops from 505.1 ºC t0 338.2 ºC.

Electron beam heating demonstration

Figure 1. Electron beam heating demonstration. Top: Stationary source. Bottom: Source moved 20.0 mm to the right over the 10 s calculation period.

HeatWave with the expanded functionality will be included in all future distributions of Xenos. A free upgrade is available to current users.

Footnotes

[1] Contact us : techinfo@fieldp.com.

[2] Field Precision home page: www.fieldp.com.