.. _xyzmesh-example: Example: 3-step decay dose calculation (reg/xyz mesh) ----------------------------------------------------- This example seeks to illustrate how one would go about using DCHAIN and PHITS to perform transport calculations utilizing the time-dependent decay spectrum from activation products. In this example, the ambient dose equivalent rate :math:`\dot{H}^*(10)` distribution in a room containing an activated target is evaluated. This idea could be expanded to a variety of real-world assessments such as how irradiation schedules, cooling times, materials compositions, and shielding configurations affect the spatial dose distributions. As the title of this section implies, this is a multi-part calculation. The three steps are: #. Run PHITS with an active source (ion beam, neutron source, etc.) and [T-Dchain] tally section(s) to determine neutron fluxes and high-energy nuclear reaction rates. #. Run DCHAIN with the files produced by [T-Dchain] to simulate activation and time-evolution of nuclide inventories and their resulting decay emission spectra. #. In PHITS, insert the decay radiation [Source] outputted to the ``*.pht`` file(s) by DCHAIN (and turn off the original active source, unless performing assessments during operation) and rerun PHITS with new tallies for desired calculations. The scenario presented here is similar to the one of the previous example but with some changes to the geometry (rectangular rather than cylindrical target and with a surrounding room geometry). In full, the following is simulated: A 1 cm diameter pencil beam of 250 MeV protons is incident on a rectangular-prism-shaped target composed of three segments: tungsten, water, and iron. This box is surrounded by beryllium on the sides and back. This structure is irradiated with a beam current of 100 nA for 10 minutes and then left to cool for 50 minutes in total. The objective will be to determine the ambient dose equivalent rate :math:`\dot{H}^*(10)` distribution in the room attributable to activation of the tungsten portion of the target at a time of 5 minutes after the end of this 10-minute irradiation. The target and room geometries are shown in :numref:`Figures %s `–:numref:`%s `; the beam is starting at z = -10 cm and traveling in the +z direction toward the target, striking the tungsten surface first. .. container:: float figure-grid figure-grid-2 :name: ex_xyzexgeo .. figure:: assets/example_calc_xyz__images__target_geometry.png :name: ex_targetgeo :width: 338.8pt Target .. figure:: assets/example_calc_xyz__images__room_geometry_cross-section.png :name: ex_roomgeo :width: 415.8pt Room The PHITS input file used and modified throughout this example calculation can be found at ``/dchain-sp/sample/3-step_dose_xyz/phits_3-step.inp`` and is also shown in :ref:`Appendix 6.3 `. It contains some commented lines and disabled sections which will be enabled later in the example. Most notably, this is the lines above the [Source] section for inserting a file, the room’s air and concrete regions in the [Cell] section, and the last [T-Track] tally. The two [G-Show] tallies were used to produce the images in :numref:`Figures %s `–:numref:`%s `, and all output presented in this example will be produced automatically by PHITS or DCHAIN. Additionally, this example seeks to illustrate usage of ``xyz`` meshes with [T-Dchain] and how it can have notable impact on results in comparison to just ``reg`` meshes, depending on the scenario. :: $ Photon source from DCHAIN output $ infl:{W_reg_target_gamma_t3.pht} $ This needs to come before the "off" source section or it will be skipped. [Source] s-type = 1 # mono-energetic axial source proj = proton # kind of incident particle :math:`\vdots` :: $ Room (void) 10 0 -16:(5 -20 15) $ room + doorway air 11 0 -15:(20 -21 #10) $ outer + maze concrete walls $ Room (with materials) $ 10 10 -0.0012 -16:(5 -20 15) $ room + doorway air $ 11 11 -2.3 -15:(20 -21 #10) $ outer + maze concrete walls :math:`\vdots` :: [T-Track] off title = Ambient dose equivalent H*(10) [mSv/hr] room map For the sake of minimizing run time, the air and concrete cells (10 and 11) are filled with void instead of their respective materials for now; though it should be noted that the room geometry would have an impact on the neutron flux within the regions. Additionally, in a serious evaluation, one would want to simulate many more histories than 10,000 (``maxcas``\ :math:`\times`\ ``maxbch``\ =10,000 here), but is reduced here for speed. There are a number of other time-saving measures implemented in this example that will be pointed out. Importantly, the [Parameters] required for [T-Dchain]/DCHAIN to function properly have been set, and the [Volume] for the regions of interest have been calculated. As a reminder, the [Volume] section is only necessary with ``reg`` meshes in [T-Dchain]; volumes are automatically calculated for ``xyz`` meshes. Variables have been declared to allow input of the beam current in nA which is then automatically converted to particles/sec for the [T-Dchain] ``amp`` parameter. There are three [T-Dchain] tallies. The first two are tallying over the same region, the tungsten target, but using differing meshes (``reg`` and ``xyz``), and the third [T-Dchain] tally is an ``xyz`` mesh covering the entire target block. The first [T-Track] tally is just a flux tally for particle fluxes in the target. After running this file in PHITS, this [T-Track] tally produces the plots shown in :numref:`Figures %s `–:numref:`%s `. :: $ required options for DCHAIN file(21) = c:/phits/dchain-sp/data # (D=c:/phits/dchain-sp/data) DCHAIN data folder jmout = 1 # (D=0) Density echo, 0:input, 1:number density e-mode = 0 # (D=0) Event generator mode is not recommended for DCHAIN igamma = 3 # (D=2) 0:No, 1:Old, 2:EBITEM, 3:EBITEM+Isomer :math:`\vdots` :: [Volume] $ required section for DCHAIN reg vol 1 c1*(2*c5)**2 2 c2*(2*c5)**2 3 c3*(2*c5)**2 4 ((c1+c2+c3+c4)*(2*c6)**2)-((c1+c2+c3)*(2*c5)**2) .. container:: float figure-grid figure-grid-2 :name: ex_fluxdists2 .. figure:: assets/example_calc_xyz__images__yz-track.png :name: ex_pflux2 :width: 377.4pt Proton flux .. figure:: assets/example_calc_xyz__images__yz-track-p2.png :name: ex_nflux2 :width: 377.4pt Neutron flux Additionally, three sets of DCHAIN input files have been generated, adopting the filenames provided to the ``file`` parameter in the [T-Dchain] tallies. Note that due to the temporary files generated by DCHAIN, only one DCHAIN simulation can be ran at a time in a single folder/directory, so they must either be ran sequentially or separated into separate folders to be ran simultaneously. For simplicity they will be ran sequentially here. In this example, DCHAIN is being controlled entirely through parameters provided to the [T-Dchain] tally (some of which are just passed along from PHITS to DCHAIN). :: [T-DCHAIN] title = W target (reg) mesh = reg reg = 1 file = W_reg_target.out # file name of dchain input file timeevo = 2 # time evolution / irradiation schedule 10.0 m 1.0 50.0 m 0.0 outtime = 4 # output times 10.0 m -1 m -5 m -50 m amp = c12 # (D=1.0) Source Intensity(source/sec) idivs = 50 # (D=50) number of calculation substebs used in irradiation iphtout = 2 # (D=1) 2 = separate [Source] output to files for each time Looking at the first [T-Dchain] tally (above), the mesh, output filename, time schedule of the beam / active source, requested output times, and nominal beam intensity are provided as is done normally as a baseline. The only additional parameter specified is ``iphtout``\ ``= 2`` which corresponds to the DCHAIN input ``PHITSOUT``; setting it to 2 causes the PHITS [Source] cards outputted by DCHAIN to be written to separate files for each output time. The convenience of this will become apparent later. Run ``W_reg_target.out`` through DCHAIN, which should only take a few seconds to complete, and the usual output files will be generated. However, as specified, the ``*.pht`` output has been split into four files corresponding to each of the requested output times. :: [T-DCHAIN] title = W target (xyz) mesh = xyz x-type = 2 nx = 8 xmin = -c5 xmax = c5 y-type = 2 ny = 8 ymin = -c5 ymax = c5 z-type = 2 nz = 3 zmin = 0 zmax = c1 file = W_xyz_target.out # file name of dchain input file timeevo = 2 # time evolution / irradiation schedule 10.0 m 1.0 50.0 m 0.0 outtime = 4 # output times 10.0 m -1 m -5 m -50 m amp = c12 # (D=1.0) Source Intensity(source/sec) idivs = 4 # (D=50) number of calculation substebs used in irradiation iphtout = 2 # (D=1) 2 = separate [Source] output to files for each time ipltmode = 1 # (D=0) enable DCHAIN 2D xy activity plot The second [T-Dchain] tally (above) nominally contains the same information except the region mesh ``reg`` has been substituted for a three-dimensional grid mesh ``xyz``. The tungsten portion of the target has dimensions 8 cm :math:`\times` 8 cm :math:`\times` 3 cm, and the ``xyz`` mesh is set to divide the region into an 8\ :math:`\times`\ 8\ :math:`\times`\ 3 grid of 1 cm\ :math:`^3` voxels. In addition to ``iphtout``\ ``= 2``, two other extra parameters have also been set. The first is ``idivs``\ ``= 4`` which controls how many steps each irradiation time step is subdivided into for the DCHAIN calculation. It has a default value of 50 which is a conservative general recommendation, particularly for neutron activation scenarios as it allows for nuclides to more accurately “grow in,” but what value of this is necessary is quite dependent on the half-lives of any important nuclides and the sizes of the time steps inputted into ``timeevo``. For a problem whose activation is primarily dependent on output from [T-Dchain]’s [T-Yield] tally for high-energy reactions, such as is the case in this example, many divisions is not necessary. It has been reduced here to only 4 to speed up the calculation. As covered in Section :ref:`4.8 `, DCHAIN’s run time scales strongly with the number of regions (192 here) and ``IDIVS``. The other extra parameter provided to [T-Dchain] was ``ipltmode``\ ``= 1`` which tells DCHAIN to also output activity map plots sliced in the xy plane of the grid mesh. Run ``W_xyz_target.out`` through DCHAIN; despite reducing ``IDIVS``, the much higher number of regions means this calculation will still take a bit longer than the previous one. In addition to the output files already discussed, ``W_xyz_target_pxy.eps`` (created from ``W_xyz_target_pxy.out``) will also be produced. As set with ``ipltmode``, this file shows plots of the mesh in the xy plane at each slice in z and at each output time. The results are shown in the :ref:`activity plots below `. As expected, activity is concentrated in the central region where the beam passed through, activity decreases with increasing cooling time, and the level of peripheral activation (outside of the beam trajectory) increases with depth in the target. .. container:: float figure-grid figure-grid-4 :name: ex_W-activity .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy_firstpage-p2.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p4.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p7.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p10.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p2.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p5.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p8.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p11.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p3.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p6.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p9.png :width: 184.8pt .. figure:: assets/example_calc_xyz__images__W_xyz_target_pxy-p12.png :width: 184.8pt :: [T-DCHAIN] title = Whole target structure mesh = xyz x-type = 2 nx = 1 xmin = -c6 xmax = c6 y-type = 2 ny = 8 ymin = -c6 ymax = c6 z-type = 2 nz = 12 zmin = 0 zmax = c1+c2+c3+c4 file = whole_target_yz-view.out # file name of dchain input file timeevo = 2 # time evolution / irradiation schedule 10.0 m 1.0 50.0 m 0.0 outtime = 4 # output times 10.0 m -10 s -1 m -50 m amp = c12 # (D=1.0) Source Intensity(source/sec) idivs = 4 # (D=50) number of calculation substebs used in irradiation ipltmode = 4 # enable DCHAIN 2D yz activity plot The third [T-Dchain] tally (above) is structured very similarly to the second one. It has an ``xyz`` mesh encompassing the entire target block with dimensions 16 cm :math:`\times` 16 cm :math:`\times` 23 cm. Since a 1 cm resolution mesh would result in 5888 regions and take a substantially longer time to run in DCHAIN, some compromises are made. Since this tally’s purpose is primarily to show how activity changes with time in each region/material, the x direction is collapsed into a single bin, and the spatial resolution in the y and z directions is halved. Thus, the mesh used is 1 :math:`\times` 8 :math:`\times` 12 with each voxel having dimensions 16 cm :math:`\times` 2 cm :math:`\times` 1.9 cm for a total of only 96 voxels. The output times have also been slightly altered. Rather than being at 0, 1, 5, and 50 minutes after the end of irradiation, the output times for this plot are 0 sec, 10 sec, 1 min, and 50 min after the end of irradiation. The ``idivs`` parameter is still set to 4, but the ``ipltmode`` parameter is now changed from 1 to 4, telling DCHAIN to output its two-dimensional activity plots in the yz plane, slicing along the x dimension (which has been reduced to a single bin in this case). Run ``whole_target_yz-view.out`` through DCHAIN, and a familiar set of outputs will be generated. However, note that this time there is only one ``*.pht`` file generated since the third [T-Dchain] tally did not alter ``iphtout`` from its default value of 1. These plots, contained in ``whole_target_yz-view_pyz.eps``, are shown in :numref:`Figures %s `–:numref:`%s `. From them we can see how activity diminishes rapidly from the Be region and more slowly in the W region. The notable fluctuations from voxel to voxel, in these and the previous plots, stem from the high variance coming from running the PHITS simulation for only 10,000 histories. Additionally, one should note that the number of histories needed for decent statistics in all of the voxels increases the finer the ``xyz`` mesh resolution is. .. container:: float figure-grid figure-grid-4 :name: ex_target-activity .. figure:: assets/example_calc_xyz__images__whole_target_yz-view_pyz.png :name: ex-target-activity-0s :width: 184.8pt t\ :sub:`cooling` = 0 sec .. figure:: assets/example_calc_xyz__images__whole_target_yz-view_pyz-p2.png :name: ex-target-activity-10s :width: 184.8pt t\ :sub:`cooling` = 10 sec .. figure:: assets/example_calc_xyz__images__whole_target_yz-view_pyz-p3.png :name: ex-target-activity-1m :width: 184.8pt t\ :sub:`cooling` = 1 min .. figure:: assets/example_calc_xyz__images__whole_target_yz-view_pyz-p4.png :name: ex-target-activity-50m :width: 184.8pt t\ :sub:`cooling` = 50 min One note worth pointing out is that for the ``xyz`` mesh, an individual voxel within the mesh can contain multiple materials. The first ``xyz`` mesh was only over the tungsten portion of the target, but this one covered the entire target block with the boundaries of the grid not selected specifically to isolate individual materials to each voxel. This can be seen in the region definition portion of the DCHAIN input deck ``whole_target_yz-view.out``; the first portion of which is shown below. Note that DCHAIN does not subdivide each voxel to keep the materials within them isolated; it just combines all of the materials within a voxel, treating it like a single region filled with a homogeneous mixture. :: ! --- calculation region details --- imat = 4 mtscore = 0 mt-list = 5 W-180 7.5882E-05 W-182 1.6756E-02 W-183 9.0514E-03 W-184 1.9376E-02 W-186 1.7975E-02 mt-list = 3 H-1 6.6880E-02 H-2 1.0033E-05 O-16 3.3446E-02 mt-list = 4 Fe-54 4.9642E-03 Fe-56 7.7928E-02 Fe-57 1.7997E-03 Fe-58 2.3951E-04 mt-list = 1 Be-9 1.2352E-01 imesh = 96 non ix iy iz fluxs volume volmat001 volmat002 volmat003 volmat004 1 1 1 1 1.2681E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 2 1 2 1 2.8259E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 3 1 3 1 5.6903E+09 6.1333E+01 2.8395E+01 0.0000E+00 0.0000E+00 3.2938E+01 4 1 4 1 1.7153E+10 6.1333E+01 2.8395E+01 0.0000E+00 0.0000E+00 3.2938E+01 5 1 5 1 1.7210E+10 6.1333E+01 2.8395E+01 0.0000E+00 0.0000E+00 3.2938E+01 6 1 6 1 5.6811E+09 6.1333E+01 2.8395E+01 0.0000E+00 0.0000E+00 3.2938E+01 7 1 7 1 2.7261E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 8 1 8 1 1.2740E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 9 1 1 2 1.5936E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 10 1 2 2 3.7344E+09 6.1333E+01 0.0000E+00 0.0000E+00 0.0000E+00 6.1333E+01 11 1 3 2 6.9048E+09 6.1333E+01 1.5863E+01 1.3542E+01 0.0000E+00 3.1929E+01 12 1 4 2 1.5202E+10 6.1333E+01 1.5863E+01 1.3542E+01 0.0000E+00 3.1929E+01 Now that DCHAIN has been ran, the PHITS [Source] sections outputted by DCHAIN into the ``*.pht`` files can be reinserted into the PHITS input to simulate transport of the decay radiation. In this example, this process for both the ``reg`` and ``xyz`` meshes is demonstrated, and then the differences in the final results between the two will be compared. Before that, a few modifications will be made to the PHITS input. Since we will be simulating decay-energy photons which are much faster to transport than 250 MeV protons and are not generating any activation we are interested in, let’s: •\ :math:`\quad` increase the number of particle histories per batch from 1000 to 100,000, :: maxcas = 100000 # (D=10) number of particles per one batch maxbch = 10 # (D=10) number of batches •\ :math:`\quad` disable the original source section, :: [Source] off s-type = 1 # mono-energetic axial source proj = proton # kind of incident particle e0 = 250.00 # energy of beam [MeV/u] •\ :math:`\quad` replace the void in the room regions 10 and 11 with air and concrete, :: $ Room (void) $10 0 -16:(5 -20 15) $ room + doorway air $11 0 -15:(20 -21 #10) $ outer + maze concrete walls $ Room (with materials) 10 10 -0.0012 -16:(5 -20 15) $ room + doorway air 11 11 -2.3 -15:(20 -21 #10) $ outer + maze concrete walls •\ :math:`\quad` and disable the [T-Dchain] tallies. :: [T-DCHAIN] off :math:`\vdots` :: [T-DCHAIN] off :math:`\vdots` :: [T-DCHAIN] off Because surveying the differences in the photon flux profile within the target area will be interesting, let’s make the first [T-Track] tally only track photons and rename its output file to something we will be able to modify to make easily distinguishable between the ``reg`` and ``xyz`` [Source] sections from DCHAIN, in this case just adding “``_reg-src``” before the file extension. (We will look at the ``reg`` [Source] first.) :: file = yz-track_reg-src.out part = photon Since we are interested in calculating the ambient dose equivalent rate distribution within the room, we will need to enable the second [T-Track] tally which has this calculation prepared. :: $ Photon source from DCHAIN output infl:{W_reg_target_gamma_t3.pht} $ This needs to come before the "off" source section or it will be skipped. [Source] off s-type = 1 # mono-energetic axial source proj = proton # kind of incident particle e0 = 250.00 # energy of beam [MeV/u] There are a few things to note about this tally. From the y axis mesh and the geometry, we can see this tally is set up to tally the dose rate in the first two meters above the concrete floor coming from particles between 1 keV and 10 MeV. The ``factor`` is used to scale the output units from pSv/sec to mSv/hr, and the multiplier is set to multiply all fluxes in region 10 (the air) by :math:`H^*(10)` fluence-to-dose conversion coefficients and to multiply fluxes in other regions by 0. The ``angel`` parameters are setting the bounds of the color bar that will appear in the output plot. Normally these are set well automatically, but since we will be comparing results from two different simulations (``reg`` versus ``xyz``) as the primary result, it is best to ensure their scales will be exactly identical. Also note that the name provided to ``file`` is formatted similarly to the one set for the other [T-Track] tally in regards to easily distinguishing between the ``reg`` and ``xyz`` mesh cases. All that remains to address is the [Source] section containing the decay radiation spectrum being emitted from the tungsten portion of the target. Above the original [Source] is a line for inserting a file via the ``infl`` function. Uncommenting that line will cause the contents of ``W_reg_target_gamma_t3.pht`` to be inserted into the PHITS input. :: [T-Track] title = Ambient dose equivalent H*(10) [mSv/hr] room map mesh = xyz x-type = 2 nx = 140 xmin = -550 xmax = 150 y-type = 2 ny = 1 ymin = -100 ymax = 100 z-type = 2 nz = 100 zmin = -150 zmax = 350 unit = 1 axis = xz e-type = 3 ne = 1 emin = 0.001 emax = 10 file = room-dose_reg-src.out set:c20[3600/1.0E+09] factor = c20 $ convert pSv/sec to mSv/hr multiplier = 6 $ number of regions using multiplier mat mset1 1 ( 0 -200 ) $ Zero out regions where we don't 2 ( 0 -200 ) $ care about the dose (inside walls 3 ( 0 -200 ) $ and the target). 4 ( 0 -200 ) 10 ( 1 -200 ) $ We only want to tally dose rate in air. 11 ( 0 -200 ) 2D-type = 7 # 1:Cont, 2:Clust, 3:Color, 4:xyz, 5:mat, 6:Clust+Cont, 7:Col+Cont gshow = 1 # 0: no 1:bnd, 2:bnd+mat, 3:bnd+reg 4:bnd+lat epsout = 1 $ automatically generate eps plot $ Set bounds of ANGEL plot color bar and make axis labels bigger angel = cmin(1.0E-5) cmax(1.0E+1) sangel = 2 x: {\Large z [cm]} y: {\Large x [cm]} Recall that due to selecting ``iphtout``\ ``= 2`` in [T-Dchain] (``PHITSOUT``\ ``= 2`` in DCHAIN), the decay radiation spectra for all regions is separated into individual files for each output time step. As the “``t3``” in the filename suggests, this file is for the third output time requested, 5 minutes after the end of the irradiation time step. Had ``iphtout`` been left at its default value of 1, we would also need to specify in the ``infl`` function the line numbers in the single ``*.pht`` file corresponding to the [Source] section for the third output time. However, before we can run this, some modifications need to be made to the [Source] section in ``W_reg_target_gamma_t3.pht``. The top portion of this file is shown below. The [Source] uses ``s-type = 5`` (same as ``s-type = 2``) wherein a region is selected and the user must construct a box around that region. PHITS will then, for every history, sample points randomly within this box until one is found which lies within the specified region. While the cell number provided to the ``reg =`` parameter here should be automatically entered correctly by DCHAIN, it cannot automatically construct this box since it is not spatially aware of cells in a ``reg`` mesh. :: c <>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<> c <>-<> no. 1 regionwise calculation data <>-<> c <>-<> region label : DUMMY001 <>-<> c <>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<> c no. 3 region= 1 output time= 15 [m] (after the last shutdown= 5 [m]) c -------------------------------------------------------------------------------- [ S o u r c e ] totfact = 2.9969E+10 $ total number of gamma-rays (n/sec/cc) x volume (cc) s-type = 5 proj = photon reg = 1 $ This cell ID should be revised for lattice or combined cell case x0 = $ Input minimum x of the region here x1 = $ Input maximum x of the region here y0 = $ Input minimum y of the region here y1 = $ Input maximum y of the region here z0 = $ Input minimum z of the region here z1 = $ Input maximum z of the region here dir = all e-type = 4 ne = 28 $ Energy Flux (arbitrary unit) $ ^^^^ Energy spectra (n/sec/cc) based on activity concentration (Bq/cc) 0.0010 1.4968E+06 0.0100 2.2244E+05 0.0200 3.3342E+06 0.0300 3.7151E+06 Thus, we must add the boundaries of this box to the [Soruce] section manually. To maximize efficiency, it is best to make the box constrain the region fairly tightly (but still contain all of the region). Since in this case the specified region is actually in the shape of a box, we can provide it with the exact dimensions as shown below. :: reg = 1 $ This cell ID should be revised for lattice or combined cell case x0 = -4 $ Input minimum x of the region here x1 = 4 $ Input maximum x of the region here y0 = -4 $ Input minimum y of the region here y1 = 4 $ Input maximum y of the region here z0 = 0 $ Input minimum z of the region here z1 = 3 $ Input maximum z of the region here One more thing to note is that the value assigned to ``totfact`` will scale the PHITS results to the actual decay emission intensity. In this case, the whole tungsten region is emitting 2.9969E+10 photons per second. This value is calculated by just multiplying the total volumetric gamma flux, 1.5609E+08 here (written to the ``*.act`` file, as shown below), by the volume of the region: 8 cm :math:`\times` 8 cm :math:`\times` 3 cm = 192 cm\ :math:`^3`. And 1.5609E+08 :math:`\gamma`/(cm\ :math:`^3\cdot`\ sec) :math:`\times` 192 cm\ :math:`^3` = 2.9969E+10 :math:`\gamma`/sec. :: total gamma-ray flux 1.5609E+08 +/- 7.5486E+06 ( 4.84%) [n/s/cc] A brief aside to explain the units conversions involved in the :math:`\dot{H}^*(10)` [T-Track] tally is now warranted. The basic units of the [T-Track] tally with ``unit = 1`` are :math:`\#`/cm\ :math:`^2` per source particle. Due to this ``totfact`` set in the [Source], the results are scaled by the emission rate whose units are :math:`\gamma`/sec. Additionally, the fluence-to-dose conversion coefficients used by the ``multiplier`` section in the [T-Track] tally have units of pSv\ :math:`\cdot`\ cm\ :math:`^2`/particle (:math:`\gamma` here). And lastly, the tally’s ``factor`` is set to change these units of time from seconds to hours and units of dose from pSv to mSv. Thus, these are all combined as shown below to obtain the tally’s final output units of mSv/hr. .. math:: \begin{equation*} \frac{\#}{\text{cm}^2} \cdot \frac{\gamma}{\text{sec}} \cdot \frac{\text{pSv}\cdot\text{cm}^2}{\gamma} \cdot \frac{\text{mSv}/\text{pSv}}{\text{hr}/\text{sec}} = \frac{\text{mSv}}{\text{hr}} \end{equation*} Now that the [Source] section is prepared, now we can run the ``dchain_ex-dose.inp`` PHITS input file through PHITS again. This will generate the files specified in the two [T-Track] tallies which are shown in :numref:`Figures %s `–:numref:`%s `. Note that while the units printed to the color bars read “Flux [1/cm\ :math:`^2`/source]” the actual units are as previously discussed. In :numref:`Figure %s `, the photon flux has units of :math:`\gamma`/(cm\ :math:`^2\cdot`\ sec), and in :numref:`Figure %s ` the ambient dose equivalent rate has units of mSv/hr. These plots can be modified by changing labels, sizes, etc. in the ``*.out`` file and running it through ANGEL; the specification of this colorbar’s label is located near the very bottom of these files. .. container:: float figure-grid figure-grid-2 :name: ex_dose-reg .. figure:: assets/example_calc_xyz__images__yz-track_reg-src.png :name: ex_gflux-reg :width: 394.6pt Activated tungsten decay photon flux [:math:`\frac{\gamma}{\text{cm}^2\cdot\text{sec}}`] in the target block .. figure:: assets/example_calc_xyz__images__room-dose_reg-src.png :name: ex_rdose-reg :width: 415.8pt :math:`\dot{H}^*(10)` [mSv/hr] room map For comparison, let’s repeat this for the ``mesh = xyz`` case. We need to make only a few modifications to the PHITS input file: •\ :math:`\quad` change the filename provided to the ``infl`` function to the ``*.pht`` file corresponding to the ``xyz`` mesh run at the same output time :: $ Photon source from DCHAIN output infl:{W_xyz_target_gamma_t3.pht} $ This needs to come before the "off" source section or it will be skipped. •\ :math:`\quad` and change the output file names in the two [T-Track] tallies to now reflect our usage of the ``xyz`` mesh instead of the ``reg`` mesh. :: file = yz-track_xyz-src.out :: file = room-dose_xyz-src.out And that’s it! Unlike with the ``reg`` mesh, the ``*.pht`` files produced for ``xyz`` meshes by DCHAIN have the source box dimensions already specified since the mesh’s grid structure is included in the PHITS tally output and passed to DCHAIN. Still, it is instructive to take a closer look at how the ``*.pht`` files produced for ``xyz`` meshes differ from those for ``reg`` meshes; a portion of the first section of ``W_xyz_target_gamma_t3.pht`` is shown below. :: c <>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<> c <>-<> no. 1 regionwise calculation data <>-<> c <>-<> region label : xyz.000001 <>-<> c <>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<> c no. 3 region= 1 output time= 15 [m] (after the last shutdown= 5 [m]) c --------------------------------------------------------------------------------- [ S o u r c e ] = 2.8071E+05 $ total number of gamma-rays (n/sec/cc) x volume (cc) in this xyz voxel at this time totfact = 2.8071E+05 $ total number of gamma-rays (n/sec/cc) x volume (cc) summed over all preceeding regions at this time; only totfact from final will be used s-type = 5 proj = photon reg = all $ This cell ID should be revised for lattice or combined cell case x0 = -4.000000E+00 $ Minimum x of the region x1 = -3.000000E+00 $ Maximum x of the region y0 = -4.000000E+00 $ Minimum y of the region y1 = -3.000000E+00 $ Maximum y of the region z0 = 0.000000E+00 $ Minimum z of the region z1 = 1.000000E+00 $ Maximum z of the region dir = all e-type = 4 ne = 25 $ Energy Flux (arbitrary unit) $ ^^^^ Energy spectra (n/sec/cc) based on activity concentration (Bq/cc) 0.0010 1.2268E+01 0.0100 2.6172E+01 0.0200 3.0745E+01 0.0300 3.5958E+01 0.0450 5.2974E-03 0.0600 5.1652E+02 0.0700 3.4553E+04 First, note that ``reg = all`` is used rather than providing ``reg =`` with a specific cell number; this is done for ``xyz`` meshes because a single voxel may contain portions of multiple cells. And perhaps the most obvious differences are the presence of `` = 2.8071E+05`` and that the file is much longer now, containing ```` sections for every voxel of the ``xyz`` mesh’s grid. This is the syntax for the multi-source in PHITS which allows simultaneous use of multiple source sections. The values of the ```` entries just tell PHITS what the intensity of each ```` section is relative to the others, and these ```` values are normalized to unity upon being read into PHITS. The actual scaling of the output is controlled by the ``totfact`` parameter. While you may notice that every ```` section has been provided with a ``totfact`` value, PHITS in fact will only use the final ``totfact`` value provided in a given [Source] section. The ``totfact`` values listed in each region is a running sum of the total emission rate for that and all preceding voxels. The top portion of the ```` entry for the final voxel is shown below. :: c <>-<> no. 192 regionwise calculation data <>-<> c <>-<> region label : xyz.000192 <>-<> c <>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<>-<> c no. 3 region= 192 output time= 15 [m] (after the last shutdown= 5 [m]) c --------------------------------------------------------------------------------- = 1.3708E+06 $ total number of gamma-rays (n/sec/cc) x volume (cc) in this xyz voxel at this time totfact = 2.9969E+10 $ total number of gamma-rays (n/sec/cc) x volume (cc) summed over all preceeding regions at this time; only totfact from final will be used s-type = 5 proj = photon reg = all $ This cell ID should be revised for lattice or combined cell case x0 = 3.000000E+00 $ Minimum x of the region x1 = 4.000000E+00 $ Maximum x of the region y0 = 3.000000E+00 $ Minimum y of the region y1 = 4.000000E+00 $ Maximum y of the region z0 = 2.000000E+00 $ Minimum z of the region z1 = 3.000000E+00 $ Maximum z of the region dir = all Note that for this final voxel the ``totfact`` value, ``2.9969E+10``, is equal to the same total photon emission rate as was in the ``reg`` mesh [Source] section, as it should be. Since it is the final ``totfact`` value listed in this [Source] section, it is the one which will be applied to scale all results. We can run the updated ``dchain_ex-dose.inp`` input file through PHITS one last time, and the value of ``totfact`` used by PHITS can be verified from the input echo section of the ``phits.out`` file, shown below. :: [ Source ] totfact = 2.99690E+10 # (D=1.0) global factor The sign of ``totfact`` in a multi-source controls how the individual ```` sections are sampled. When ``totfact`` is positive, all source particles will be born with the same weight, and the rate at which each ```` section is sampled scales with its value relative to the values of other ```` sections. (```` sections with higher values are sampled more often.) When ``totfact`` is negative, all ```` sections are sampled equally but the weights of the particles spawned by each ```` section are adjusted to reflect the relative intensity of that ```` section against the other ```` sections. By default it is positive, but if you wish to change this since the activities in an ``xyz`` mesh can often vary multiple orders of magnitude and you may wish to ensure that every single voxel is sufficiently sampled, only the final ``totfact`` value needs to be made negative. It should also be noted that setting ``iphtout``\ ``= 2`` in [T-Dchain] (``PHITSOUT``\ ``= 2`` in DCHAIN) is essentially mandatory for ``xyz`` meshes if you have more than one output time requested in ``outtime`` in [T-Dchain] (``ITOUT`` in DCHAIN). This is because DCHAIN performs its calculations for one region at a time, meaning normally output for all time steps is written for a region before advancing to the next region. This is how the ``*.act`` file is structured and how the ``*.pht`` files are structured by default, if not split apart as done here. However, for this type of simulation, we want to isolate the decay spectra for all regions at a single time rather than for a single region at all times. If we had used ``iphtout``\ ``= 1`` here, we would have to tediously provide the ``infl`` function in PHITS with 192 different sets of line ranges containing all of the different ```` sections at our single desired time or just only include the single output time in ``outtime``. Anyways, now that PHITS has been ran with the ``xyz`` mesh multi-source, we can observe its outputs, shown in :numref:`Figures %s `–:numref:`%s `. .. container:: float figure-grid figure-grid-2 :name: ex_dose-xyz .. figure:: assets/example_calc_xyz__images__yz-track_xyz-src.png :name: ex_gflux-xyz :width: 394.6pt Activated tungsten decay photon flux [:math:`\frac{\gamma}{\text{cm}^2\cdot\text{sec}}`] in the target block .. figure:: assets/example_calc_xyz__images__room-dose_xyz-src.png :name: ex_rdose-xyz :width: 415.8pt :math:`\dot{H}^*(10)` [mSv/hr] room map For the ``reg`` mesh, the activation and resulting decay emissions were evenly distributed throughout the entire tungsten region. However, now with the ``xyz`` mesh, we can see that the emissions are indeed more narrowly concentrated along the path of the beam as had already been witnessed in the 2D ``xyz`` activity map plots produced by DCHAIN. This results in a noticeable self-shielding effect as a large fraction of photons are born toward the center (further from the exterior) of the tungsten target and are attenuated before escaping the target block. This has the ultimate effect of reducing the ambient dose equivalent rate in the room due to a greater fraction of emissions not escaping the target. One may note that the statistical clarity of these results seems slightly worse than those for the ``reg`` mesh despite the same number of original beam proton and decay photon histories being simulated in each case, and this is the reason. The same number of source photons are produced in both simulations, but for the ``xyz`` mesh case more are born toward the center of the target and are attenuated before escaping into the room. Of course, statistics can always be improved by just running PHITS for more particle histories. :numref:`Figures %s `–:numref:`%s ` show how these plots, for both ``reg`` and ``xyz`` cases, change when increasing the number of histories from 1 million to 50 million, making the effects of self-shielding on the photon flux around the target and resulting reduction in :math:`\dot{H}^*(10)` in the room more apparent. (The colorbar labels have been adjusted in ANGEL too.) .. container:: float figure-grid figure-grid-4 :name: ex_dose-xyz2 .. figure:: assets/example_calc_xyz__images__high-stats__fixed-label__yz-track_reg-src.png :name: ex_gflux-reg2 :width: 394.6pt ``reg`` - Activated tungsten decay photon flux [:math:`\frac{\gamma}{\text{cm}^2\cdot\text{sec}}`] in the target block .. figure:: assets/example_calc_xyz__images__high-stats__fixed-label__room-dose_reg-src.png :name: ex_rdose-reg2 :width: 415.8pt ``reg`` - :math:`\dot{H}^*(10)` [mSv/hr] room map .. figure:: assets/example_calc_xyz__images__high-stats__fixed-label__yz-track_xyz-src.png :name: ex_gflux-xyz2 :width: 394.6pt ``xyz`` - Activated tungsten decay photon flux [:math:`\frac{\gamma}{\text{cm}^2\cdot\text{sec}}`] in the target block .. figure:: assets/example_calc_xyz__images__high-stats__fixed-label__room-dose_xyz-src.png :name: ex_rdose-xyz2 :width: 415.8pt ``xyz`` - :math:`\dot{H}^*(10)` [mSv/hr] room map For the sake of illustration, :numref:`Figures %s `–:numref:`%s ` show how the ``xyz`` plots, still with 50 million histories, would change if the ``totfact`` value had been made negative instead (by changing the sign of the final ``totfact`` in ``W_xyz_target_gamma_t3.pht``). Perhaps counterintuitively, the statistical clarity actually gets worse despite the voxels nearest the edge of the target being sampled more. Studying the ``*.gso`` and ``*.act`` files, one will see that the photons contributing the most to dose rate and with the most penetrating power are ones from nuclides created through high-energy beam reactions in the tungsten, and in the ``xyz`` mesh case they are only being spawned a significant fraction of the time in the voxels closest to the center of the target. Additionally, due to setting ``totfact`` as a negative value in order to sufficiently sample all voxels, it does mean that these voxels of the highest activity are only sampled a fraction of the time. Referring back to the :ref:`activity plots above `, we see in the :math:`8\times8` grid that the central 4 voxels contain the vast majority of this beam reaction activation, meaning in PHITS they are only sampled for 6.25% of histories here, less often than when ``totfact`` was positive. .. container:: float figure-grid figure-grid-2 :name: ex_dose-xyz2_pt .. figure:: assets/example_calc_xyz__images__high-stats__neg-totfact__fixed-label__yz-track_xyz-src.png :name: ex_gflux-xyz2_pt :width: 394.6pt ``xyz`` - Activated tungsten decay photon flux [:math:`\frac{\gamma}{\text{cm}^2\cdot\text{sec}}`] in the target block .. figure:: assets/example_calc_xyz__images__high-stats__neg-totfact__fixed-label__room-dose_xyz-src.png :name: ex_rdose-xyz2_pt :width: 415.8pt ``xyz`` - :math:`\dot{H}^*(10)` [mSv/hr] room map Still, depending on your specific scenario, the mesh resolution, how concentrated or more evenly distributed the activation is, the mechanisms contributing to activation, what decay emissions are of concern, and what you wish to evaluate in a secondary PHITS simulation, setting ``totfact`` negative may yield improved results. For instance, in an object which is activated throughout but with increased activation in its core/center but is so large/dense that decay emissions would almost never escape from the center (perhaps a very thick cask which had previously held a strong neutron source), setting ``totfact`` negative would allow the more abundant voxels which are closer to the surface, but are less activated, to be sampled more often. If performing some dose assessment near this object, this would result in better statistics in the same amount of particle histories since fewer histories are “wasted” by spawning in the voxels in the smaller central region, which is the most activated but whose emissions almost never escape and contribute to the tally results. In general, the idea is to set the sign of ``totfact`` such that the number of particles contributing to your tallies of interest is maximized; though it is not always immediately obvious which sign results in more efficient use of computation time.