7.16. [ T-WWG ] section

This tally gives parameters used for the [Weight Window] section. The tally serves as a Weight Window Generator (WWG) to automatically obtain effective settings of [Weight Window] for a user-defined virtual space. When using this function, a test calculation using an input file including the [t-wwg] tally must be performed. Long calculations can then be performed using the [Weight Window] section generated from this test calculation. In the case of mesh=reg, you have to give the volume of each cell in the [volume] section. When you perform the restart calculation or use the sumtally function with an input file including [t-wwg], set axis other than wwg and ireschk=1 in the [parameters] section.

This tally scores the fluence of particles in each history \(j\), particle \(k\), energy mesh \(l\), and geometrical mesh \(m\), denoted as \(\Phi_{jklm}\). Then, at the end of each batch, it determines the lower limit of their weight window, \(W_{klm,\rm low}\), as written below:

(7.16.1)\[W_{klm,\rm low} = C_{kl} \left[ \left( \Sigma_j \Phi_{jklm} \right)_{k,l,x{\rm \%ile}} + \Sigma_j \Phi_{jklm} \right],\]

where \(\left( \Sigma_j \Phi_{jklm} \right)_{k,l,x{\rm \%ile}}\) is the \(x\)-percentile value of \(\Sigma_j\Phi_{jklm}\) for particle \(k\) and energy mesh \(l\) among all geometrical mesh \(m\). The value of \(x\) can be determined by an input parameter named pedestal whose default value is 0.1 that indicates 10 percentiles. The basic idea of this equation is to make \(W_{klm,\rm low}\) proportional to \(\Sigma_j\Phi_{jklm}\), but the introduction of the pedestal term, \(\left( \Sigma_j \Phi_{jklm} \right)_{k,l,x{\rm \%ile}}\), is required to avoid extremely low \(W_{klm,\rm low}\) that could cause excessive particle splitting.

\(C_{kl}\) in Eq. (7.16.1) is the normalization constant which depends on an input parameter named normww as written below:

(7.16.2)\[\begin{split}C_{kl} = \begin{cases} \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } & {\rm for \, normww=0} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{lm,{\rm max}} } & {\rm for \, normww=1} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{klm,{\rm max}} } & {\rm for \, normww=2} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } \frac{ \left( \Sigma_{jm} \Phi_{jklm} \right) } {\left( \Sigma_{jm} \Phi_{jklm} \right)_{l,{\rm min}}} & {\rm for \, normww=-1} \\ \frac{W_{\rm max}}{ \left( \Sigma_j \Phi_{jklm} \right)_{m,{\rm max}} } \frac{ \left( \Sigma_{jm} \Phi_{jklm} \right) } {\left( \Sigma_{jm} \Phi_{jklm} \right)_{kl,{\rm min}}} & {\rm for \, normww=-2} \end{cases}\end{split}\]

Here, the subscripts \(m, \max\), \(lm, \max\), and \(klm, \max\) denote the maximum values over the respective indices, while \(l, \min\) and \(lm, \min\) denote the minimum values over the corresponding indices. For example, \(\left( \sum_j \Phi_{jklm} \right)_{lm, \mathrm{max}}\) represents the maximum value of \(\left( \sum_j \Phi_{jklm} \right)\) over all energy groups l and spatial meshes m, and thus takes different values for each particle k. When normww is 0, 1, or 2, normalization is performed based on the maximum fluence at the mesh points of m, lm, and klm, respectively. On the other hand, when normww is -1 or -2, a correction term is added based on the energy group (for -1) or particle-energy group (for -2) that gives the minimum value of the sum over m, i.e., the total fluence over the entire spatial domain. As a result, biases in the total fluence among energy groups (for -1) or among particle-energy groups (for -2) are corrected, yielding weight window values that lead to a more balanced spatial distribution of Monte Carlo particles. Thus, if you want to focus on high-energy neutron transport, you can narrow the energy-mesh width of high-energy neutrons and set normww=1, so that \(W_{klm,\rm low}\) for that high-energy group will be smaller overall.

The fundamental concept of [t-wwg] is to generate a [weight window] section that can distribute MC fluence as uniformly as possible to all locations within the regions defined in [t-wwg]. On the other hand, if you want to evaluate the dose in a certain region of interest (ROI) with high accuracy in a short time, you can generate a [weight window] that can navigate particles to the ROI by combining the history counter and the maximum biasing factor, \(B_{\rm max}\), defined by the chwei parameter. For activating this algorithm, you must define the [counter] section to increase a counter value when a particle enters the ROI and specify the chwei(i) parameter in [t-wwg], where \(i\) is the counter number defined in [counter]. Then, this tally determines the biasing factor for each mesh, \(B_{klm}\), as written below:

(7.16.3)\[B_{klm} = \frac{\Sigma_j H_{ij} \varphi_{jklm}} {\left( \Sigma_j H_{ij} \varphi_{jklm} \right)_{k,l,{\rm max}} } \left[ B_{\rm max} -1 \right] +1.\]

where \(H_{ij}\) is the \(i\)-th history-counter value of history \(j\), and \(\varphi_{jklm}\) is the MC fluence for the condition. The strategy behind this equation is to maximize the biasing factor up to \(B_{\rm max}=\) chwei(i) at the geometrical mesh \(m\) that gives the maximum value of \(\Sigma_j H_{ij} \varphi_{jklm}\), that is the location where the MC particles pass through most frequently before or after arriving at the ROI, while minimizing the factor to 1 where \(\Sigma_j H_{ij} \varphi_{jklm}=0\). The lower limit of the weight window, \(W_{klm,\rm low}\), is calculated by replacing \(\Sigma_j \varphi_{jklm}\) with \(\Sigma_j \Phi_{jklm}/B_{klm}\) in Eqs. (7.16.1) and (7.16.2). Please see the literature [1] and the lecture notes [2] for more details.

The [t-wwg] parameters are formatted as follows.

Table 7.16.1 mesh

value	explanation
reg, xyz, tet	Mesh type. Only reg, xyz, and tet can be set. A mesh type subsection is required below this option.

Table 7.16.2 part

value	explanation
all (default), particle name	Tally particle.

Table 7.16.3 material

value	explanation
(optional)	Specify materials for scoring.
all	all is the default, same as no definition.
number of materials	To set number of materials, define the material numbers in the next line. The number of materials can be set as a negative value, in which case the specified materials are not included for scoring.
(next line)	2 5 8 Material numbers. If the list does not fit on one line, it can be split across multiple lines without using the line-connection symbol.

Table 7.16.4 e-type

value	explanation
1, 2, 3, 4, 5	Energy mesh. See Section 6.6.1 for mesh subsection format.

Table 7.16.5 t-type

value	explanation
1, 2, 3, 4, 5 (optional)	Time mesh. See Section 6.6.1 for mesh subsection format.

Table 7.16.6 unit

value	explanation
1 (optional)	1: [ \(1/\mathrm{cm}^2/\mathrm{source}\) ]

Table 7.16.7 axis

value	explanation
eng, reg, tet, t	x axis value of output data.
xy, yz, xz	2-dimensional. They are the same as [t-track].
wwg	To obtain the [Weight Window] section, set axis=wwg. To perform restart calculation, set two or more axis parameters, and then set the first axis to reg, eng, or t.

Only reg, xyz, and tet can be set as mesh in this tally, because the parameters in [Weight Window] are defined only for these mesh types. Note that the surfaces of xyz mesh should be slightly shifted from surfaces of cells. Please be especially careful if you are using parameters to link geometry boundary and xyz mesh ranges. After version 3.34, if the maximum and minimum xyz mesh values for [weight window] and [t-wwg] are defined as integers, the minimum value is shifted by \(-1*\) deltxyz and the maximum value by \(+2*\) deltxyz from the integer, see Q2.9: I got many warning messages about lost particles. What should I do?. mesh=tet can be used only for geometry containing tetrahedral-mesh geometry. The weight window factors for each constituent tetrahedron element for the tetrahedral-mesh geometry can be obtained. In short, axis should be set to wwg to obtain the parameters for [Weight Window]. Although eng, reg, tet, xy, yz, xz, and t can also be set, their results are not related with [Weight Window]. Note that restart calculation cannot be performed when axis=xy, yz, or xz are selected because the _err file is not created by this tally.

Table 7.16.8 file

value	explanation
file name	Define output file names. This is required by each setting of axis.

Table 7.16.9 resfile

value	explanation
(optional, D=file)	Define a file name of the past tally in the restart calculation. Even if several axis parameters were defined, specify only one resfile.

Table 7.16.10 factor

value	explanation
(optional, D=1.0)	Normalization factor.

Table 7.16.11 title

value	explanation
(optional)	Title.

Table 7.16.12 angel

value	explanation
(optional)	ANGEL parameters.

Table 7.16.13 rshow

value	explanation
0 (default), 1, 2, 3	When mesh=reg and axis=xy,yz,xz, region border (1), material name (2), and region name (3) are plotted using this option. A xyz mesh section must be added below this option.

Table 7.16.14 resol

value	explanation
1 (default)	This option multiplies the region line resolution by a factor of resol with the gshow or rshow option set to define the line thickness.

Table 7.16.15 width

value	explanation
0.5 (default)	The option defines the line thickness.

Table 7.16.16 volume

value	explanation
(optional)	This option defines the volume for each region for reg mesh. Volume definitions are required below this option. Values defined in [volume] are used in the case of no definition. If special description such as using (  ) is used to specify a region in reg=, its internally defined region number is output in input echo when this volume subsection is not defined.

Table 7.16.17 epsout

value	explanation
0 (default), 1, 2	When epsout=1, results are plotted into eps files. The eps file is named by replacing the extension with .eps. When epsout=2, error bars are also displayed in the eps file, except for the 2-dimensional type, axis=xy, yz, xz, rz.

Table 7.16.18 ph5out

value	explanation
0 (default), 1	When set to 1, the tally results are output in the PHITS HDF5 file format (ph5). The file name is the same as the output file, with the extension changed to .ph5.

Table 7.16.19 trcl

value	explanation
(optional)	Coordinate transformation number or definition for r-z or xyz mesh.

Table 7.16.20 gslat

value	explanation
2 (default)	Option to draw lattice or tetrahedral geometry boundary lines when gshow or rshow is specified.
0	Not drawing.
1	Drawing. Note that the line of voxel phantoms may not be drawn clearly.
2	Not drawing the lines in the same cell.
3	Not drawing the lines in the same material. In this case, the boundary of two adjacent cells filled with the same material is not also drawn, even if not lattice or tetrahedral geometry. When gshow>=3, cell and lattice numbers are now shown.

Table 7.16.21 stdcut

value	explanation
(optional, D=-1)	Threshold value of STD cut off.

When specifying stdcut, PHITS automatically stop the calculation depending on values of STD, standard deviation. This function is available when stdcut is positive and itall=0,1 is set in [parameters] section. When all relative values of STD of the tally result are larger than 0 and smaller than stdcut at the last of one batch, the calculation is stopped. If stdcut in two or more tally sections is set, all the results of the tally sections have to satisfy the conditions in order for the function to work.

Table 7.16.22 pedestal

value	explanation
(optional, D=0.1)	Percentile value, \(x\) in Eq. (7.16.1), divided by 100 for calculating the pedestal of \(W_{klm,{\rm low}}\).

Table 7.16.23 normww

value	explanation
(optional, D=-2)	Normalization option, see Eq. (7.16.2). 0: No adjustment. \(\pm 1\): Adjust \(C_{kl}\) so that the MC fluences in each energy group are equivalent. \(\pm 2\): Adjust \(C_{kl}\) so that the MC fluences in each particle and energy group are equivalent.

When normww is set to a positive value, the overall weight window values become lower, so the weight of each particle is controlled primarily by particle splitting. Conversely, if normww is set to a negative value, the overall weight window values become higher, leading to more frequent Russian roulette. For calculation conditions with large variance per history, such as the case of sources with a wide energy distribution, a negative value will give better statistics per unit computational time, since the calculation time per history will be shorter. This is the reason why normww=-2 is selected as the default value. However, care should be taken to define particle type and energy mesh in the case of normww=-2, as most particles are killed by Russian roulette if there is even one group with extremely low particle fluence in each particle and energy group. For the part parameter, it is recommended to set all or neutron photon in the case of neutron and photon only transport calculations. Regarding the energy groups, it is recommended to set the number of energy bins just a few, or less, with a slightly narrower group width for the energy bands that contribute the most to the tally. If you are not sure which setting is optimal, we encourage you to use part=all and ne=1. In that case, the calculation results will not depend on normww.

Table 7.16.24 wwmax

value	explanation
(optional, D=0.99)	Normalization constant corresponding to the maximum value of \(W_{klm,{\rm low}}\), see Eq. (7.16.2).

Table 7.16.25 updateww

value	explanation
1 (default)	Specifies whether to automatically update the lower limit of the weight window during the calculation. If a [weight window] section already exists in the input, this is automatically set to 0.
0	[t-wwg] only outputs the weight window file, and the weight window during the calculation does not change.
1	The weight window is automatically updated from the tally results of each batch (weight window optimization within a single calculation).

Table 7.16.26 wwgcut

value

explanation

(optional, D=0.1)

Threshold for terminating the calculation upon convergence of the weight window. It plays a role similar to that of stdcut in other tallies. When updateww=1 (the WW values are updated at the end of each batch), the calculation conditions change during the run, so in principle the statistical uncertainty cannot be evaluated. Instead, PHITS outputs the fraction of conditions whose WW value changes by wwgband or more between batches (the change fraction). When this value becomes smaller than wwgcut for two consecutive batches, PHITS judges that the calculation results have converged and terminates the calculation. If wwgcut is set in two or more [t-wwg] sections, this function operates when all of them satisfy the condition.

Table 7.16.27 wwgband

value	explanation
(optional, D=0.1)	In the convergence check by wwgcut, a condition is counted as a “changed condition” when the ratio of the WW lower limit before and after the update between batches is equal to or greater than this value.

Table 7.16.28 ctmin(i)

value	explanation
(optional, D=-9999)	Minimum value for the i-th counter.

Table 7.16.29 ctmax(i)

value	explanation
(optional, D=9999)	Maximum value for the i-th counter.

Table 7.16.30 chmin(i)

value	explanation
(optional, D=-9999)	Minimum value for the i-th history-counter. This parameter cannot be specified in the batch variance mode, istdev=1.

Table 7.16.31 chmax(i)

value	explanation
(optional, D=9999)	Maximum value for the i-th history-counter. This parameter cannot be specified in the batch variance mode, istdev=1.

Table 7.16.32 chwei(i)

value

explanation

(optional, D=0)

Magnitude of particle-navigation capability in the generated [weight windows] using the i-th history counter, that is \(B_{\rm max}\) in Eq. (7.16.3). Only one of these parameters can be defined per [t-wwg], even if their counter IDs are different. The restart calculation or use of the sumtally function is now allowed when this parameter is defined. The recommended value is roughly 20 to 200; the larger the value, the stronger the particle navigation becomes.

Table 7.16.33 chplane

value

explanation

(optional, D=all)

The plane, xy, yz, or xz, for calculating the maximum value of \(\left( \Sigma_j H_{ij} \varphi_{jklm} \right)_{k,l,{\rm max}}\) in Eq. (7.16.3). The default value is all, which indicates that the maximum value is calculated from all geometrical meshes. When using particle-navigation function for duct streaming calculations, specifying a plane perpendicular to the duct will reduce the computation time. This parameter is valid only when mesh=xyz and chwei is defined. Note that the traceback function (tracemem > 0) enables more effective particle navigation, so this parameter is no longer recommended for use.

Table 7.16.34 tracemem

value	explanation
(optional, D=0)	The number of elements of the traceback memory allocated per thread (in units of real(8) = 8 bytes), used by the traceback function that improves the accuracy of the particle-navigation function. If it is set to 0, the traceback function is not used. However, if chwei is defined and no other [t-wwg] defines tracemem, it is automatically set to 1000000.

In weight window generation using the particle-navigation function (chwei), the weight window is generated by weighting the histories of particles that reach the region targeted by [counter]. However, in the method up to version 3.36, if even a single particle reached the target region within one history, particles were navigated along the tracks of all particles generated in that history, so the function did not work effectively when the number of secondary particles was large. Therefore, from version 3.37, by recording all particle tracks in a temporary memory (buffer) and reanalyzing them after the history ends (traceback), particle navigation can now be limited to the tracks of the parent particles of the particles that reached the target region.

tracemem (> 0) can be used only when all of the following conditions are satisfied. Otherwise, the calculation stops with an error.

• The mesh type must be xyz or tet (reg is not allowed).

• It must be used together with chwei (the traceback assumes the chwei mode).

• Only one [t-wwg] tally per input can specify tracemem > 0.

The traceback buffer is allocated for the size specified by tracemem times the number of parallel threads (since each OpenMP thread has its own independent buffer). Therefore, the larger the number of parallel threads, the more memory is required. If the buffer becomes insufficient during recording, the recording continues while replenishing the memory as needed (reallocating and expanding it). However, if the replenishment (reallocation) fails due to insufficient memory, a warning is output and the subsequent traceback recording for that history is aborted (the tracks of the aborted history are not reflected in the weight window generation). Setting tracemem sufficiently large can avoid the occurrence of replenishment and the abortion of recording due to insufficient memory.

Information on the memory used for the traceback is output to the standard output and the console output file at the end of the calculation in the following format.

tracemem summary for [t-wwg]: max used =       <value>, allocated =       <value>
         reallocations =   <value>, allocation failures =   <value>

• max used: The maximum number of buffer elements used in any thread (peak usage).

• allocated: The finally allocated buffer size (the number of elements per thread; if it was expanded by replenishment, that value).

• reallocations: The number of times the buffer was replenished (reallocated and expanded).

• allocation failures: The number of times the replenishment (reallocation) failed.

If max used is much smaller than allocated (the initial tracemem), tracemem can be reduced, and if reallocations or allocation failures occur, it is recommended to increase tracemem.

Table 7.16.35 elowthre

value	explanation
(optional, D=0)	Threshold energy [MeV] for the low-energy unbiasing method. Below this energy, the lower limit of the weight window is multiplied by elowbias. If it is set to 0, the low-energy unbiasing method is not used. Note that the number of energy bins, ne, must be set to 1 when using this method.

Table 7.16.36 elowbias

value	explanation
(optional, D=5.0)	The multiplying factor for the lower limit of the weight window for energies below elowthre.

An example output of this tally is as follows.

Listing 7.16.1 Example output of [t-wwg]

[ Weight Window ]
     mesh =  reg

set:  c101[0.0]  c102[c101+2.13496E-08]  c103[1.25478E-04]
set:  c104[0.0]  c105[c104+9.29569E-09]  c106[2.08987E-04]

    part = neutron
     eng =  2
          1.00000E-03   1.00000E+03

    reg          ww1                    ww2
    1   (8.94092E-05+c101)/c103  (2.08987E-04+c104)/c106
    2   (1.25478E-04+c101)/c103  (1.26817E-04+c104)/c106
    3   (8.53835E-05+c101)/c103  (5.78131E-05+c104)/c106

In this sample, the [Weight Window] parameters for neutrons with two energy bins are output. In principle, it is not necessary to change these parameters, but c101 or c104 should be specified when adding a constant value for each Weight Window.

[1]

Sato et al. Nucl. Instr. and Meth. B, 2024; 557: 165535.

[2]

phits/lecture/advanced/weightB and phits/lecture/advanced/shielding