- Email: [email protected]
AMOS – An effective tool for adjoint Monte Carlo photon transport
NIM B Beam Interactions with Materials & Atoms
Nuclear Instruments and Methods in Physics Research B 251 (2006) 326–332 www.elsevier.com/locate/nimb
AMOS – An effective tool for adjoint Monte Carlo photon transport Dorothea Gabler a, Ju¨rgen Henniger a
a,*
, Uwe Reichelt
b
Radiation Physics Group, Institute of Nuclear and Particle Physics, Technische Universita¨t Dresden, IKTP, 01062 Dresden, Germany b Institute of Nuclear and Hadron Physics, Forschungszentrum Rossendorf, Germany Received 20 March 2006; received in revised form 30 May 2006 Available online 6 September 2006
Abstract In order to expand the photon version of the Monte Carlo transport program AMOS to an adjoint photon version, AMOS Pt, the theory of adjoint radiation transport is reviewed and evaluated in this regard. All relevant photon interactions, photoelectric effect, coherent scattering, incoherent scattering and pair production, are taken into account as proposed in the EPDL 97. In order to simulate pair production and to realise physical source terms with discrete energy levels, an energy point detector is used. To demonstrate the qualification of AMOS Pt a simple air-over-ground problem is simulated by both the forward and the adjoint programs. The results are compared and total agreement is shown. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Monte Carlo; Radiation transport; Photon transport; Transport equation; Adjoint; Air-over-ground
1. Introduction The theory of adjoint transport has been long since developed for photon and neutron transport [1–3]. Adjoint transport is already provided by various Monte Carlo transport codes which are often restricted to adjoint neutron transport and are mostly held in multigroup representation [4–11,13,14]. If the latter is the case adjoint transport is rather easily implemented by simply using the adjoint transition matrices instead of the forward ones. A totally different approach is necessary for continuous energy mode. Only recently a continuous energy adjoint photon transport code was developed using analytical forms for the different interaction cross sections [12]. But with this solution low energy problems cannot be handled satisfactorily, when for example the Klein–Nishina formula for incoherent scattering is not sufficient. AMOS [15] is a validated Monte Carlo tool especially designed for photon and electron transport providing a *
Corresponding author. Tel.: +49 351 463 32479; fax: +49 351 463 37040. E-mail address: [email protected] (J. Henniger). 0168-583X/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nimb.2006.07.005
continuous energy scale for particle energies while using tabulated interaction data grouped into up to 5000 energy intervals of variable width. This large number of energy groups leads to the benefit of speeding up the simulation, because an interpolation is not necessary and the interaction data can be used directly. Hence the here presented algorithm of adjoint photon transport does not originate from neutron or adjoint neutron transport. It is instead optimised for photons. Elementarily all relevant photon interactions, photoelectric effect, coherent and incoherent scattering, and pair production, must be treated upon the same fundamentals as in forward transport. For example, the cross sections, the coherent form factors and incoherent scattering functions are to be regarded as given in the EPDL 97 [16] which provides all photon interaction data used in AMOS. In order to inherit the AMOS Monte Carlo method, adjoint group data have to be generated in the same form as the forward analogon. This leads to the same profit of widely unrestricted geometries and fast simulation. However, to be able to realise pair production and discrete energy source terms an energy point detector similar to a spatial point detector had to be implemented slowing down
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
327
the simulation to a certain degree. Still, adjoint transport is always preferable for geometries with large source and small detector volumes in phase space.
wE,X multiplied with the appropriate response function fw!R Z Z R ¼ . . . d3 r dE d2 X dtf w!R ðr; E; X; tÞ wE;X ðr; E; X; tÞ:
2. Formalism of forward and adjoint Monte Carlo transport
ð6Þ
A detailed description of the forward and adjoint formalism for radiation transport can be found for example in [1]. The shortened version here presented aims for pure photon transport in the non-analogue Monte Carlo simulation program AMOS and it is therefore already an interpretation of any general formalism. The basis of forward Monte Carlo radiation transport is the integral Boltzmann equation in operator notation
In the case of the effect of interest being the fluence U, the response function equals 1=ðV det Rt ðr; EÞÞ r 2 V det ; w!U f ðr; E; V det Þ ¼ ð7Þ 0 r 62 V det
b Tb vE;X ðr0 ; E0 ; X0 ; t0 Þ vE;X ðr; E; X; tÞ ¼ sE;X ðr; E; X; tÞ þ C ð1Þ of the emergent particle density vE,X with the source term sE,X. The subscripts, in this case E and X, denote the varib and Tb ables the radiation quantities are differentiated to. C are integral operators acting to the right on a function b is the collision operator f(r, E, X, t). C b ¼ Cðr; b E0 ! E; X0 ! X; tÞf ðr; E0 ; X0 ; tÞ ¼ Cf
Rðr; E0 ! E; X0 ! XÞ f ðr; E0 ; X0 ; tÞ Rt ðr; E0 Þ
Z Z
2
0
Z
dX 4p
1
dE0
0
ð2Þ
and Tb the transport operator Tb f ¼ Tb ðr0 ! r; E0 ; X0 ; t0 ! tÞf ðr0 ; E0 ; X0 ; t0 Þ Z 1 h i 0 0 dL Rt ðr; E0 Þ ebðr;L;E ;X Þ f ðr0 ; E0 ; X0 ; t0 Þ ¼
ð3Þ
RL with bðr; L; E0 ; X0 Þ ¼ 0 dL0 Rt ðr L0 X0 ; E0 Þ and r 0 = r LX 0 , t 0 = t L/v. L and L 0 are path lengths to integrate over, v is the particle velocity. The differential cross section R(r, E 0 ! E, X 0 ! X) describes the scattering from (E 0 , X 0 ) into (E, X). Rt(r, E 0 ) is the total cross section for all interactions which may take place at r for photons of energy E 0 . As already noted, differentiated radiation field quantities are subscripted with the differentiating parameters. Instead, differentiated cross sections as d3 Rðr; E0 Þ dE d2 X
Rt ðr; E0 Þ ¼ Rs ðr; E0 Þ þ Ra ðr; E0 Þ;
ð4Þ
are marked by arrows. In contrast to this notation the arrows in transport and collision operators are there for reasons of illustration and do not denote differentiation. As resulting quantity an effect of interest R measured by a detector is desired. In general it is calculated via the collision density wE;X ðr; E; X; tÞ ¼ Tb ðr0 ! r; E; X; t0 ! tÞ vE;X ðr0 ; E; X; t0 Þ;
ð8Þ
where Ra is the absorption cross section and Rs is the scattering cross section excluding scattering into E = 0 keV. With this the collision operator can be written as Z 1 Z Z b ¼ Cf d2 X0 dE0 4p
0
Rðr; E0 ! E; X0 ! XÞ ¼
and yields the sensitive detector volume Vdet. Eq. (7) directly shows the advantage of using wE,X instead of vE,X for the integration to R. Because fw!U is nonzero only in the detector region, the integral is limited to the volume of the physical detector. Before further discussion, a few assumptions have to be made. First, absorption processes are taken to be scattering into E = 0 keV. Second, they are the only interactions scattering downward into E = 0 keV. And last it is presumed that only downscattering occurs, possible upscattering processes are not taken into account. Then the total cross section Rt(r,E 0 ) can be split into two terms
0
Rs ðr; E0 ! E; X0 ! XÞ þ Ra ðr; E0 Þ 0 0 f ðr; E ; X ; tÞ Rt ðr; E0 Þ ð9Þ and because particles of energy E = 0 keV do not contribute to any effect of interest the second summand can be disregarded. This results in Z Z Z 1 2 0 b Cf ¼ dX dE0 4p
0
Rs ðr; E0 Þ Rs ðr; E0 ! E; X0 ! XÞ 0 0 f ðr; E ; X ; tÞ Rt ðr; E0 Þ Rs ðr; E0 Þ ð10Þ and the first term in the right hand side integral can be formulated as the so-called non-absorption probability PNA(r,E 0 ) = 1 Ra(r,E 0 )/Rt(r,E 0 ). The second term can be expanded into a sum of products with two terms each " Z Z Z 1 X 2 0 0 b Cf ¼ dX dE P NA ðr; E0 Þ ðP j ðr; E0 Þ 0 8j # 4p fj ðr; E0 ! E; X0 ! XÞÞ f ðr; E0 ; X0 ; tÞ
ð11Þ
ð5Þ which results from applying Tb once more to the emergent particle density. Then, R results from the integration over
constituting all scattering interactions which change energy and direction of a particle from (E 0 , X 0 ) to (E, X) with
328
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
E > 0 keV. For each interaction j fj(r, E 0 E, X 0 X)/ 0 Rj(r, E ) is the probability density function (pdf) for (E, X) with E > 0 keV. Pj(r, E 0 ) = Rj(r, E 0 )/Rs(r, E 0 ) is the conditional probability for the respective scattering interaction under the condition that no absorption takes place. After having settled that upscattering is excluded, particles with energies E = 0 keV can be disregarded. From here on all particle density quantities and cross sections are defined to except energies E = 0 keV and Eq. (1) is applicable for these new quantities. The path of an emergent particle in AMOS is therefore followed by first choosing a free path length out of the transport kernel. Then a scattering process is compelled by multiplying PNA to the particle’s weight, whereby absorption is regarded implicitly. Again it has to be said, that this formulation is only for non-analogue transport, treating the absorption processes as pure leaks. Now, one of the possible scattering processes is chosen via the probabilities Pj. The new energy and direction are sampled from the corresponding fj. In AMOS a particle has energy and direction in a continuous scale (double precision floatingpoint numbers) and the change through scattering is calculated via the exact formulas. The scattering data, which is to say all PNA, Pj and fj, are instead provided for energy groups of E 0 . Between those no interpolation is carried out. After this short description of the forward formulation, the transition to the adjoint can now be carried out. Regarding Eq. (6) the adjoint collision density wþ E;X is fixed by Z Z R ¼ . . . d3 r dE d2 X dt½wþ ð12Þ E;X sE;X and therefore has to obey the adjoint Boltzmann equation w!R wþ þ C þ T þ wþ E;X ¼ f E;X :
ð13Þ
The problem with these adjoint quantities is that the collision and transport operators look totally different compared to the forward operators. Therefore the Monte Carlo simulation may not be derived from Eqs. (12) and (13). Instead, an adjoint integral Boltzmann equation similar to (1) can be set up ~ E;X ðr; E; X; tÞ ¼ f~ w!R ðr; E; X; tÞ þ C ~ E;X ðr0 ; E0 ; X0 ; t0 Þ; ~ T~ w w ð14Þ [4], using a modified adjoint collision density ~ E;X ðr; E; X; tÞ ¼ Rt ðr; EÞ wþ ðr; E; X; tÞ. The accordw E;X ingly modified response function f~ w!R ðr; E; X; tÞ ¼ Rt ðr; EÞ f w!R ðr; E; X; tÞ acts as the adjoint source term. Further, the modified adjoint collision and transport operators are contained, now similar to the forward ones. Here the integration to the effect of interest could be for~ E;X multiplied with mulated through the integration over w an appropriately modified source term of wE,X. But instead, as in the forward case, in practice R is calculated via a different adjoint radiation field quantity, the modified adjoint emergent particle density
~ E;X ðr0 ; E; X; t0 Þ w ~vE;X ðr; E; X; tÞ ¼ T~ ðr0 ! r; E; X; t0 ! tÞ ; Rt ðr0 ; EÞ ð15Þ and the modified physical source term is simply given by ~sE;X ðr; E; X; tÞ ¼ sE;X ðr; E; X; tÞ:
ð16Þ
So again the integration to R Z Z R ¼ . . . d3 r dE d2 X dt ~sE;X ðr; E; X; tÞ ~vE;X ðr; E; X; tÞ ð17Þ is limited to a definite volume by ~sE;X – the volume of the physical source instead of the physical detector. The modified adjoint transport operator has exactly the same form as the forward transport operator T~ ¼ Tb
ð18Þ
and is only labelled differently to keep consistency in the adjoint formulation. Hence the Monte Carlo simulation of the free path between two collision points is forward and adjoint exactly the same. In contrast, the modified adjoint collision operator ~ ¼ Cðr; ~ E Cf E0 ; X X0 tÞf ðr; E0 ; X0 tÞ Z Z Z 1 ¼ d2 X0 dE0 4p
0
Rs ðr; E ! E0 ; X ! X0 Þ 0 0 f ðr; E ; X ; tÞ ; Rt ðr; E0 Þ
ð19Þ
differs substantially from the forward one. In particular, ~ and the different arrow directions in the arguments of C Rs are to be pointed out. At this point it is substantial to discriminate between the different meanings of designations. The physical quantity Rs ðr; E ! E0 ; X ! X0 Þ ¼
d3 Rs ðr; EÞ dE0 d2 X0
ð20Þ
is the differential scattering cross section, in which the arrow directions are set by definition. The arrow directions ~ are put this way for reasons of illustration. Looking in C at the adjoint Monte Carlo particles the primed quantities, for example E 0 , still describe their properties before the interaction process. But in the adjoint everything happens ‘‘backward’’ compared to the physical reality, which is supposed to be pointed out by the backward arrow directions. On the same argumentation as in the forward case, the ~ can be expanded into a sum of products kernel of C ~ E Cðr;
E0 ; X
X0 ; tÞ f ðr; E0 ; X0 ; tÞ ¼
h P TNA ðr; E0 Þ Rj ðP Tj ðr; E0 Þ fjT ðr; E
Z Z
d2 X0
4p
E0 ; X
Z
1
dE0
0
i X0 ÞÞ f ðr; E0 ; X0 ; tÞ ;
ð21Þ
where the sum runs over all upscattering interactions which may change the particle’s energy and direction from (E 0 ,X 0 ) to (E,X) regarding E 0 > 0 keV. In this formalism the fjT are
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
the adjoint pdfs of (E,X) for the single scattering processes j fjT ðr; E
E0 ; X
X0 Þ ¼
Rj ðr; E ! E0 ; X ! X0 Þ RTj ðr; E0 Þ
with the adjoint cross sections Z Z Z 1 d2 X dERj ðr; E ! E0 ; X ! X0 Þ: RTj ðr; E0 Þ ¼ 4p
ð22Þ
ð23Þ
0
As the forward ones, in general they do not depend on the particle’s initial direction X 0 but on its energy E 0 and the scattering material at r. P Tj are the conditional probabilities for the single scattering processes P Tj ðr; E0 Þ ¼
RTj ðr; E0 Þ RTs ðr; E0 Þ
ð24Þ
because RTs is the sum of all adjoint scattering cross sections RTj , again excluding absorption. P TNA is the ‘‘adjoint nonabsorption probability’’ P TNA ðr; E0 Þ ¼
RTs ðr; E0 Þ ; Rt ðr; E0 Þ
ð25Þ
which as a physical quantity is not analogous to the forward one because it is not a probability at all. It must be pointed out that RTs is the adjoint scattering cross section but Rt is the forward total cross section. These two quantities are not directly connected to each other. But still P TNA is part of the adjoint collision operator and accounts for the loss of particles over absorption processes. Since P TNA is not a probability there is no way of playing an analogous adjoint Monte Carlo game. But as AMOS is a non-analogue Monte Carlo simulation program this is of no consequence, and the adjoint radiation transport may be carried out exactly the same way as the forward simulation. Therein P TNA is multiplied to the particle’s weight. Through P Tj the scattering process j and through fjT energy and direction are chosen. Again the scales of energy and direction are held continuous while group data is used for the interactions. The main difference between forward and adjoint transport are the interaction data which must be newly processed for the adjoint transport. One more remark about the adjoint non-absorption probability has to be stated at this point. As to be seen in the next section, P TNA is not restricted to values <1 but instead can reach all values P0. Because the non-absorption probability is multiplied to the particle weight, this may lead to high variances, which must be critically observed at every adjoint Monte Carlo simulation. 3. Adjoint scattering processes Photon interaction processes regarded in AMOS are photoelectric effect, coherent scattering, incoherent scattering and pair production. The only interaction treated as absorption process is the photoelectric effect, which is
329
therefore avoided in the simulation by multiplying the non-absorption probability. Coherent and incoherent scattering and pair production are the possible scattering interactions. In adjoint transport a particle’s energy is changed ‘‘upward’’ instead of ‘‘downward’’. This is the most obvious difference to forward transport. As forward upscattering is disregarded the same must be done with adjoint downscattering. The scattering processes are completely reversed – in energy as well as in direction. But since in general no direction in the material is featured and polarisation is disregarded, the reversal of directions is of no effect. Hence coherent scattering is forward and adjoint exactly the same. As proposed in the EPDL, the coherent cross section is calculated including the coherent form factor, and the real and imaginary anomalous scattering function. For its differential only the coherent form factor is regarded. 3.1. Incoherent scattering According to the Compton formula the forward energy loss in an incoherent scattering is correlated with the polar scattering angle # or its cosine l = cos # through E¼
E0 ; 1 þ E0 =E0 ð1 lÞ
ð26Þ
and because of l 2 [1, 1] the outcome energy is bounded to the interval Emin ¼
E0 6 E 6 E0 1 þ 2 E0 =E0
ð27Þ
depending on the initial energy E 0 . In this notation E0 = 511 keV is the electron rest energy. This gives the energy relation E¼
E0 1 E =E0 ð1 lÞ 0
ð28Þ
for the adjoint process. Note again that the primed variables sign the properties of a Monte Carlo particle before its interaction, so, compared to the forward energy formula (26), the variables E and E 0 are interchanged. In Eq. (28) for energies E 0 P E0/2 large scattering angles with l just above the singularity at lmin lead to very large energies and l < lmin would formally result in negative energies. So now, instead of the scattering energies the polar scattering angle cosine is restricted to lmin < l 6 1 with the energy dependence ( 1 E0 < E20 lmin ¼ ; ð29Þ 1 EE00 E0 P E20 so that for l ! lmin, l > lmin the energy E diverges to infinity for E 0 P E0/2, see Fig. 1. Accordingly, the integration of the adjoint differential cross section
330
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
E0 ¼
Emax ; 1 þ 2 Emax =E0
ð32Þ
and decreases to zero at E = Emax as shown in Fig. 2. According to these equations, the adjoint data for incoherent scattering can be calculated to be used in AMOS the same way as the forward data. In addition, the different formulas for energy or scattering angle have to be implemented. 3.2. Pair production
Fig. 1. Relative energy of incoherently scattered photons in dependence of polar scattering angle cosine for several primary photon energies.
RTi ðr; E0 ;
l
Þ ¼ qa ðrÞ pr2e
E E0 þ 1 þ l2 E0 E
Sðr; E; lÞ; ð30Þ
including the scattering function S(r,E,l) and E = E(E 0 ,l) from (28), would result in an infinite adjoint cross section RTi ðr; E0 Þ. Therefore, a maximum energy level has to be set above which scattering is not possible. This is also suitable considering that in physical problems there is always a highest source energy above which no particles are emitted and hence adjoint scattering into these energies would be irrelevant. Finally and most important, an upper energy boundary is needed to ensure convergence by stopping an otherwise endless adjoint calculation. Thus a problem dependent maximum energy Emax is set, to which there is a minimum scattering angle cosine ( 1 E0 < 1þ2EEmax max =E0 0 lmin ðE ; Emax Þ ¼ : ð31Þ 0 1 þ EEmax EE00 E0 P E20 Now the adjoint cross section may be calculated. Depending on the maximum scattering energy it has a point of discontinuity in the first derivation at
Fig. 2. Comparison of forward and adjoint incoherent scattering cross sections (germanium) with two different maximum energies.
In pure photon transport pair production at the nucleus is treated as a scattering process from energies above 2E0 into E0. The two outcome photons are diametrically emitted. Adjoining this is not possible, because this would mean that two photons, each of the energy E0, meet at a point which they approach from exactly opposite directions. Instead, the forward process can be described as single photon scattering, still into the energy E0 but with no direction bound. For this process the forward differential production cross section is Rpð1Þ ðr; E0 ! E; X0 ! XÞ 1 ð33Þ 4p with the common forward cross section Rp(r, E 0 ) of pair production, the term 1/4p accounting for the isotropic emission and the Dirac d-distribution forcing the outcome energy to equal E0. The increase of the cross section by the factor 2 associates to the two annihilation photons. Their emission is not simulated in correlation to each other. Observing Eq. (33) in the adjoint view leads to the same problem as the adjoint incoherent scattering. A maximum energy has to be set, which cuts the differential cross section ¼ 2 Rp ðr; E0 Þ dðE E0 Þ
RTp ðr; E E0 ; X X0 ; Emax Þ 1 2 Rp ðr; EÞ dðE0 E0 Þ 4p ¼ 0
E < Emax ; E P Emax
ð34Þ
so that the integral is limited. In the simulation what is required is that adjunctions, the adjoint Monte Carlo particles, with energy E0, are available to enter adjoint pair production. This is accomplished by using a so-called energy point detector through which particles are incoherently scattered into E0. In contrast to the usual point detector, the outcome energy instead of a space point is the desired parameter. So, incoherent scattering is compelled by multiplying the conditional probability for adjoint incoherent scattering P Ti from Eq. (24). Then, if possible, the particle is scattered into E = E0. The scattering angle cosine l is already set by Eq. (28) and the pdf fiT at l has to be multiplied to the particle’s weight. The azimuthal scattering angle is chosen out of 1/2p. After the following transport to the next interaction point, the particle may undergo coherent scattering. If this interaction process is drawn a further transport follows and again the possibility of coherent
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
scattering. However, if incoherent scattering is drawn, this history ends. But at each of these interaction points, a subhistory is split off and forced to succumb pair production. Here the conditional probability, but this time for adjoint pair production, P Tp is multiplied to the weight. The outcome energy and direction are chosen via fpT . Such sub-histories with the production of adjoint annihilation photons are split off at every interaction point and run parallel to the original Monte Carlo particle history. As the scattering by adjoint pair production leads to high energies above 2E0 and as incoherent scattering and even more coherent scattering at these energies are very improbable, these sub-histories are destined to be very short. Hence the computational effort, though not minimised, is still kept within limits. This energy point detector is also used for realising discrete energy sources where in the adjoint simulation the same problem arises. Before every collision a sub-history is split off and, if possible, incoherently scattered into a discrete source energy. If its following interaction point is inside the source volume it is counted, or else this pseudo particle is abandoned. This is not a particularly efficient procedure, but regarding the special geometries at which adjoint transport is used and thereby its advantage over the forward transport this drawback is tolerable. 4. Example To prove the applicability of AMOS Pt a simple airover-ground problem is calculated by both AMOS and AMOS Pt and the results are compared. Requested is the attenuation by soil cover of varying thickness d on contaminated ground. The contamination is uranium ore with an activity density of 1 Bq/g of its reference nuclide 226Ra. As all following nuclides of the uranium decay series are treated the highest source energy is 2448 keV emitted from 214 Bi. Hence pair production must be considered. The reference point of the detector is located 1 m above the surface. Only by taking advantage of the homogeneous geometry the detector volume may be expanded significantly, so a
331
Fig. 3. Comparison of forward and adjoint results of the air-over-ground problem.
forward simulation may be performed with tolerable statistical error. The resulting ambient dose equivalent rates H_ ð10Þ are compared in Fig. 3, in which the absolute agreement within the statistical errors can be observed. For example the values without soil cover are forward : H_ ð10Þðd ¼ 0 cmÞ ¼ ð506:2 9:4Þ nSv=h; adjoint : H_ ð10Þðd ¼ 0 cmÞ ¼ ð506:8 7:2Þ nSv=h: Besides the verification of the reliability of AMOS Pt, one advantage of the adjoint simulation can also be demonstrated through this example. Via the adjoint calculation it is possible to generate precise energy spectra. As an example, Fig. 4 shows the spectral ambient dose equivalent rate for again d = 0 cm. Some significant source energies are signed out. A comparison with the analogue forward calculated spectrum would bring no gain because of its large uncertainty for high energies. In contrast, the adjoint simulation yields good statistics for large energies, degrading slightly to lower energies. As can be observed in Fig. 4, it has to be noted that only the radiation of the uranium ore is simulated. Nuclides such as potassium or thorium, which must not be disregarded to simulate reality, are not taken into account here.
Fig. 4. Adjoint simulated spectrum of the ambient dose equivalent rate for d = 0 cm.
332
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
5. Conclusions and discussion In this paper the realisation of adjoint photon transport in AMOS Pt is presented. Its validity is shown through comparison with forward simulated results. The use of the energy point detector firstly leads to a loss of simulation speed. But still the adjoint simulation is clearly preferable at geometries with large physical source and very small physical detector volume and energy range. The energy point detector acts only in one dimension, in which moreover, because of the down- or respectively upscattering condition, one direction is featured. Instead, in a forward simulation a 3D space point detector would have to be used to score in the small detector volume. The loss of speed in the adjoint simulation is therefore not nearly as big as it would be in the forward analogon. As soon as the symmetry is broken forward Monte Carlo transport in air-over-ground problems cannot be carried out as in the here presented example. Acknowledgement We would like to thank David Legrady for his constructive comments on the manuscript. References [1] I. Lux, L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press, Boca Raton, Florida, 1991. [2] R.R. Coveyou, V.R. Cain, K.J. Yost, Adjoint and importance in Monte Carlo application, Nucl. Sci. Eng. 27 (1967) 219. [3] D.C. Irving, The adjoint Boltzmann equation and its simulation by Monte Carlo, Nucl. Eng. Des. 14 (1971) 273.
[4] M.B. Emmett, The MORSE Monte Carlo Radiation Transport Code System, ORNL-4972, UC-32 – Mathematics and Computers, Oak Ridge National Laboratory, Oak Ridge, Tennessee, 1975. [5] T.E. Booth, F.B. Brown, J.S. Bull, L.J. Cox, R.A. Forster, J.T. Goorley, H.G. Hughes, R.D. Mosteller, R.E. Prael, E.C. Selcow, A. Sood, J.E. Sweezy, R.F. Barrett, S.E. Post, T.L. Roberts, MCNP – A General Monte Carlo N-Particle Transport Code, Version 5, Volume I: Overview and Theory, Los Alamos National Laboratory Report LA-UR-03-1987, April 2003. [6] M.H. Kalos, Monte Carlo integration of the adjoint gamma-ray transport equation, Nucl. Sci. Eng. 33 (1968) 284. ¨ zisik, Y. Yener, On the inverse problem for [7] C.E. Siewert, M.N. O multigroup neutron transport theory, Tech. Notes (rev. 1976) 95. [8] A. De Matteis, R. Simonini, A new approach to the adjoint Boltzmann equation, Nucl. Sci. Eng. 65 (1978) 93. [9] A. De Matteis, R. Simonini, A Monte Carlo biasing scheme for adjoint photon transport, Nucl. Sci. Eng. 67 (1978) 309. [10] S.N. Cramer, Adjoint gamma-ray estimation to the surface of a cylinder – Analysis of a remote reprocessing facility, Nucl. Sci. Eng. 79 (1981) 417. [11] R.W. Roussin, F.A. Schmidt, Adjoint calculations of coupled neutron and gamma-ray transport through concrete slabs, Nucl. Eng. Des. 15 (1971) 319. [12] J.E. Hoogenboom, Adjoint Monte Carlo photon transport in continuous energy mode with discrete photons from annihilation, in: PHYSOR, 2000. [13] J.E. Hoogenboom, Methodology of continuous energy adjoint Monte Carlo for neutron, photon, and coupled neutron-photon transport, Nucl. Sci. Eng. 143 (2003) 99. [14] N.M. Borisov, M.P. Panin, Generalized particle concept for adjoint Monte Carlo calculations of coupled gamma-ray-electron-positron transport, Nucl. Sci. Eng. 150 (2005) 284. [15] J. Henniger, Grundlagen der Abschirmberechnungen mit MonteCarlo-Methoden, Wiss. Z. TU Dresden 39 (3) (1990) 9. [16] D.E. Cullen, J.H. Hubbell, L. Kissel, EPDL97; Lawrence Livermore National Laboratory, Livermore, CA, UCRL50400, Vol. 6, Rev. 5, 1997.
Nuclear Instruments and Methods in Physics Research B 251 (2006) 326–332 www.elsevier.com/locate/nimb
AMOS – An effective tool for adjoint Monte Carlo photon transport Dorothea Gabler a, Ju¨rgen Henniger a
a,*
, Uwe Reichelt
b
Radiation Physics Group, Institute of Nuclear and Particle Physics, Technische Universita¨t Dresden, IKTP, 01062 Dresden, Germany b Institute of Nuclear and Hadron Physics, Forschungszentrum Rossendorf, Germany Received 20 March 2006; received in revised form 30 May 2006 Available online 6 September 2006
Abstract In order to expand the photon version of the Monte Carlo transport program AMOS to an adjoint photon version, AMOS Pt, the theory of adjoint radiation transport is reviewed and evaluated in this regard. All relevant photon interactions, photoelectric effect, coherent scattering, incoherent scattering and pair production, are taken into account as proposed in the EPDL 97. In order to simulate pair production and to realise physical source terms with discrete energy levels, an energy point detector is used. To demonstrate the qualification of AMOS Pt a simple air-over-ground problem is simulated by both the forward and the adjoint programs. The results are compared and total agreement is shown. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Monte Carlo; Radiation transport; Photon transport; Transport equation; Adjoint; Air-over-ground
1. Introduction The theory of adjoint transport has been long since developed for photon and neutron transport [1–3]. Adjoint transport is already provided by various Monte Carlo transport codes which are often restricted to adjoint neutron transport and are mostly held in multigroup representation [4–11,13,14]. If the latter is the case adjoint transport is rather easily implemented by simply using the adjoint transition matrices instead of the forward ones. A totally different approach is necessary for continuous energy mode. Only recently a continuous energy adjoint photon transport code was developed using analytical forms for the different interaction cross sections [12]. But with this solution low energy problems cannot be handled satisfactorily, when for example the Klein–Nishina formula for incoherent scattering is not sufficient. AMOS [15] is a validated Monte Carlo tool especially designed for photon and electron transport providing a *
Corresponding author. Tel.: +49 351 463 32479; fax: +49 351 463 37040. E-mail address: [email protected] (J. Henniger). 0168-583X/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nimb.2006.07.005
continuous energy scale for particle energies while using tabulated interaction data grouped into up to 5000 energy intervals of variable width. This large number of energy groups leads to the benefit of speeding up the simulation, because an interpolation is not necessary and the interaction data can be used directly. Hence the here presented algorithm of adjoint photon transport does not originate from neutron or adjoint neutron transport. It is instead optimised for photons. Elementarily all relevant photon interactions, photoelectric effect, coherent and incoherent scattering, and pair production, must be treated upon the same fundamentals as in forward transport. For example, the cross sections, the coherent form factors and incoherent scattering functions are to be regarded as given in the EPDL 97 [16] which provides all photon interaction data used in AMOS. In order to inherit the AMOS Monte Carlo method, adjoint group data have to be generated in the same form as the forward analogon. This leads to the same profit of widely unrestricted geometries and fast simulation. However, to be able to realise pair production and discrete energy source terms an energy point detector similar to a spatial point detector had to be implemented slowing down
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
327
the simulation to a certain degree. Still, adjoint transport is always preferable for geometries with large source and small detector volumes in phase space.
wE,X multiplied with the appropriate response function fw!R Z Z R ¼ . . . d3 r dE d2 X dtf w!R ðr; E; X; tÞ wE;X ðr; E; X; tÞ:
2. Formalism of forward and adjoint Monte Carlo transport
ð6Þ
A detailed description of the forward and adjoint formalism for radiation transport can be found for example in [1]. The shortened version here presented aims for pure photon transport in the non-analogue Monte Carlo simulation program AMOS and it is therefore already an interpretation of any general formalism. The basis of forward Monte Carlo radiation transport is the integral Boltzmann equation in operator notation
In the case of the effect of interest being the fluence U, the response function equals 1=ðV det Rt ðr; EÞÞ r 2 V det ; w!U f ðr; E; V det Þ ¼ ð7Þ 0 r 62 V det
b Tb vE;X ðr0 ; E0 ; X0 ; t0 Þ vE;X ðr; E; X; tÞ ¼ sE;X ðr; E; X; tÞ þ C ð1Þ of the emergent particle density vE,X with the source term sE,X. The subscripts, in this case E and X, denote the varib and Tb ables the radiation quantities are differentiated to. C are integral operators acting to the right on a function b is the collision operator f(r, E, X, t). C b ¼ Cðr; b E0 ! E; X0 ! X; tÞf ðr; E0 ; X0 ; tÞ ¼ Cf
Rðr; E0 ! E; X0 ! XÞ f ðr; E0 ; X0 ; tÞ Rt ðr; E0 Þ
Z Z
2
0
Z
dX 4p
1
dE0
0
ð2Þ
and Tb the transport operator Tb f ¼ Tb ðr0 ! r; E0 ; X0 ; t0 ! tÞf ðr0 ; E0 ; X0 ; t0 Þ Z 1 h i 0 0 dL Rt ðr; E0 Þ ebðr;L;E ;X Þ f ðr0 ; E0 ; X0 ; t0 Þ ¼
ð3Þ
RL with bðr; L; E0 ; X0 Þ ¼ 0 dL0 Rt ðr L0 X0 ; E0 Þ and r 0 = r LX 0 , t 0 = t L/v. L and L 0 are path lengths to integrate over, v is the particle velocity. The differential cross section R(r, E 0 ! E, X 0 ! X) describes the scattering from (E 0 , X 0 ) into (E, X). Rt(r, E 0 ) is the total cross section for all interactions which may take place at r for photons of energy E 0 . As already noted, differentiated radiation field quantities are subscripted with the differentiating parameters. Instead, differentiated cross sections as d3 Rðr; E0 Þ dE d2 X
Rt ðr; E0 Þ ¼ Rs ðr; E0 Þ þ Ra ðr; E0 Þ;
ð4Þ
are marked by arrows. In contrast to this notation the arrows in transport and collision operators are there for reasons of illustration and do not denote differentiation. As resulting quantity an effect of interest R measured by a detector is desired. In general it is calculated via the collision density wE;X ðr; E; X; tÞ ¼ Tb ðr0 ! r; E; X; t0 ! tÞ vE;X ðr0 ; E; X; t0 Þ;
ð8Þ
where Ra is the absorption cross section and Rs is the scattering cross section excluding scattering into E = 0 keV. With this the collision operator can be written as Z 1 Z Z b ¼ Cf d2 X0 dE0 4p
0
Rðr; E0 ! E; X0 ! XÞ ¼
and yields the sensitive detector volume Vdet. Eq. (7) directly shows the advantage of using wE,X instead of vE,X for the integration to R. Because fw!U is nonzero only in the detector region, the integral is limited to the volume of the physical detector. Before further discussion, a few assumptions have to be made. First, absorption processes are taken to be scattering into E = 0 keV. Second, they are the only interactions scattering downward into E = 0 keV. And last it is presumed that only downscattering occurs, possible upscattering processes are not taken into account. Then the total cross section Rt(r,E 0 ) can be split into two terms
0
Rs ðr; E0 ! E; X0 ! XÞ þ Ra ðr; E0 Þ 0 0 f ðr; E ; X ; tÞ Rt ðr; E0 Þ ð9Þ and because particles of energy E = 0 keV do not contribute to any effect of interest the second summand can be disregarded. This results in Z Z Z 1 2 0 b Cf ¼ dX dE0 4p
0
Rs ðr; E0 Þ Rs ðr; E0 ! E; X0 ! XÞ 0 0 f ðr; E ; X ; tÞ Rt ðr; E0 Þ Rs ðr; E0 Þ ð10Þ and the first term in the right hand side integral can be formulated as the so-called non-absorption probability PNA(r,E 0 ) = 1 Ra(r,E 0 )/Rt(r,E 0 ). The second term can be expanded into a sum of products with two terms each " Z Z Z 1 X 2 0 0 b Cf ¼ dX dE P NA ðr; E0 Þ ðP j ðr; E0 Þ 0 8j # 4p fj ðr; E0 ! E; X0 ! XÞÞ f ðr; E0 ; X0 ; tÞ
ð11Þ
ð5Þ which results from applying Tb once more to the emergent particle density. Then, R results from the integration over
constituting all scattering interactions which change energy and direction of a particle from (E 0 , X 0 ) to (E, X) with
328
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
E > 0 keV. For each interaction j fj(r, E 0 E, X 0 X)/ 0 Rj(r, E ) is the probability density function (pdf) for (E, X) with E > 0 keV. Pj(r, E 0 ) = Rj(r, E 0 )/Rs(r, E 0 ) is the conditional probability for the respective scattering interaction under the condition that no absorption takes place. After having settled that upscattering is excluded, particles with energies E = 0 keV can be disregarded. From here on all particle density quantities and cross sections are defined to except energies E = 0 keV and Eq. (1) is applicable for these new quantities. The path of an emergent particle in AMOS is therefore followed by first choosing a free path length out of the transport kernel. Then a scattering process is compelled by multiplying PNA to the particle’s weight, whereby absorption is regarded implicitly. Again it has to be said, that this formulation is only for non-analogue transport, treating the absorption processes as pure leaks. Now, one of the possible scattering processes is chosen via the probabilities Pj. The new energy and direction are sampled from the corresponding fj. In AMOS a particle has energy and direction in a continuous scale (double precision floatingpoint numbers) and the change through scattering is calculated via the exact formulas. The scattering data, which is to say all PNA, Pj and fj, are instead provided for energy groups of E 0 . Between those no interpolation is carried out. After this short description of the forward formulation, the transition to the adjoint can now be carried out. Regarding Eq. (6) the adjoint collision density wþ E;X is fixed by Z Z R ¼ . . . d3 r dE d2 X dt½wþ ð12Þ E;X sE;X and therefore has to obey the adjoint Boltzmann equation w!R wþ þ C þ T þ wþ E;X ¼ f E;X :
ð13Þ
The problem with these adjoint quantities is that the collision and transport operators look totally different compared to the forward operators. Therefore the Monte Carlo simulation may not be derived from Eqs. (12) and (13). Instead, an adjoint integral Boltzmann equation similar to (1) can be set up ~ E;X ðr; E; X; tÞ ¼ f~ w!R ðr; E; X; tÞ þ C ~ E;X ðr0 ; E0 ; X0 ; t0 Þ; ~ T~ w w ð14Þ [4], using a modified adjoint collision density ~ E;X ðr; E; X; tÞ ¼ Rt ðr; EÞ wþ ðr; E; X; tÞ. The accordw E;X ingly modified response function f~ w!R ðr; E; X; tÞ ¼ Rt ðr; EÞ f w!R ðr; E; X; tÞ acts as the adjoint source term. Further, the modified adjoint collision and transport operators are contained, now similar to the forward ones. Here the integration to the effect of interest could be for~ E;X multiplied with mulated through the integration over w an appropriately modified source term of wE,X. But instead, as in the forward case, in practice R is calculated via a different adjoint radiation field quantity, the modified adjoint emergent particle density
~ E;X ðr0 ; E; X; t0 Þ w ~vE;X ðr; E; X; tÞ ¼ T~ ðr0 ! r; E; X; t0 ! tÞ ; Rt ðr0 ; EÞ ð15Þ and the modified physical source term is simply given by ~sE;X ðr; E; X; tÞ ¼ sE;X ðr; E; X; tÞ:
ð16Þ
So again the integration to R Z Z R ¼ . . . d3 r dE d2 X dt ~sE;X ðr; E; X; tÞ ~vE;X ðr; E; X; tÞ ð17Þ is limited to a definite volume by ~sE;X – the volume of the physical source instead of the physical detector. The modified adjoint transport operator has exactly the same form as the forward transport operator T~ ¼ Tb
ð18Þ
and is only labelled differently to keep consistency in the adjoint formulation. Hence the Monte Carlo simulation of the free path between two collision points is forward and adjoint exactly the same. In contrast, the modified adjoint collision operator ~ ¼ Cðr; ~ E Cf E0 ; X X0 tÞf ðr; E0 ; X0 tÞ Z Z Z 1 ¼ d2 X0 dE0 4p
0
Rs ðr; E ! E0 ; X ! X0 Þ 0 0 f ðr; E ; X ; tÞ ; Rt ðr; E0 Þ
ð19Þ
differs substantially from the forward one. In particular, ~ and the different arrow directions in the arguments of C Rs are to be pointed out. At this point it is substantial to discriminate between the different meanings of designations. The physical quantity Rs ðr; E ! E0 ; X ! X0 Þ ¼
d3 Rs ðr; EÞ dE0 d2 X0
ð20Þ
is the differential scattering cross section, in which the arrow directions are set by definition. The arrow directions ~ are put this way for reasons of illustration. Looking in C at the adjoint Monte Carlo particles the primed quantities, for example E 0 , still describe their properties before the interaction process. But in the adjoint everything happens ‘‘backward’’ compared to the physical reality, which is supposed to be pointed out by the backward arrow directions. On the same argumentation as in the forward case, the ~ can be expanded into a sum of products kernel of C ~ E Cðr;
E0 ; X
X0 ; tÞ f ðr; E0 ; X0 ; tÞ ¼
h P TNA ðr; E0 Þ Rj ðP Tj ðr; E0 Þ fjT ðr; E
Z Z
d2 X0
4p
E0 ; X
Z
1
dE0
0
i X0 ÞÞ f ðr; E0 ; X0 ; tÞ ;
ð21Þ
where the sum runs over all upscattering interactions which may change the particle’s energy and direction from (E 0 ,X 0 ) to (E,X) regarding E 0 > 0 keV. In this formalism the fjT are
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
the adjoint pdfs of (E,X) for the single scattering processes j fjT ðr; E
E0 ; X
X0 Þ ¼
Rj ðr; E ! E0 ; X ! X0 Þ RTj ðr; E0 Þ
with the adjoint cross sections Z Z Z 1 d2 X dERj ðr; E ! E0 ; X ! X0 Þ: RTj ðr; E0 Þ ¼ 4p
ð22Þ
ð23Þ
0
As the forward ones, in general they do not depend on the particle’s initial direction X 0 but on its energy E 0 and the scattering material at r. P Tj are the conditional probabilities for the single scattering processes P Tj ðr; E0 Þ ¼
RTj ðr; E0 Þ RTs ðr; E0 Þ
ð24Þ
because RTs is the sum of all adjoint scattering cross sections RTj , again excluding absorption. P TNA is the ‘‘adjoint nonabsorption probability’’ P TNA ðr; E0 Þ ¼
RTs ðr; E0 Þ ; Rt ðr; E0 Þ
ð25Þ
which as a physical quantity is not analogous to the forward one because it is not a probability at all. It must be pointed out that RTs is the adjoint scattering cross section but Rt is the forward total cross section. These two quantities are not directly connected to each other. But still P TNA is part of the adjoint collision operator and accounts for the loss of particles over absorption processes. Since P TNA is not a probability there is no way of playing an analogous adjoint Monte Carlo game. But as AMOS is a non-analogue Monte Carlo simulation program this is of no consequence, and the adjoint radiation transport may be carried out exactly the same way as the forward simulation. Therein P TNA is multiplied to the particle’s weight. Through P Tj the scattering process j and through fjT energy and direction are chosen. Again the scales of energy and direction are held continuous while group data is used for the interactions. The main difference between forward and adjoint transport are the interaction data which must be newly processed for the adjoint transport. One more remark about the adjoint non-absorption probability has to be stated at this point. As to be seen in the next section, P TNA is not restricted to values <1 but instead can reach all values P0. Because the non-absorption probability is multiplied to the particle weight, this may lead to high variances, which must be critically observed at every adjoint Monte Carlo simulation. 3. Adjoint scattering processes Photon interaction processes regarded in AMOS are photoelectric effect, coherent scattering, incoherent scattering and pair production. The only interaction treated as absorption process is the photoelectric effect, which is
329
therefore avoided in the simulation by multiplying the non-absorption probability. Coherent and incoherent scattering and pair production are the possible scattering interactions. In adjoint transport a particle’s energy is changed ‘‘upward’’ instead of ‘‘downward’’. This is the most obvious difference to forward transport. As forward upscattering is disregarded the same must be done with adjoint downscattering. The scattering processes are completely reversed – in energy as well as in direction. But since in general no direction in the material is featured and polarisation is disregarded, the reversal of directions is of no effect. Hence coherent scattering is forward and adjoint exactly the same. As proposed in the EPDL, the coherent cross section is calculated including the coherent form factor, and the real and imaginary anomalous scattering function. For its differential only the coherent form factor is regarded. 3.1. Incoherent scattering According to the Compton formula the forward energy loss in an incoherent scattering is correlated with the polar scattering angle # or its cosine l = cos # through E¼
E0 ; 1 þ E0 =E0 ð1 lÞ
ð26Þ
and because of l 2 [1, 1] the outcome energy is bounded to the interval Emin ¼
E0 6 E 6 E0 1 þ 2 E0 =E0
ð27Þ
depending on the initial energy E 0 . In this notation E0 = 511 keV is the electron rest energy. This gives the energy relation E¼
E0 1 E =E0 ð1 lÞ 0
ð28Þ
for the adjoint process. Note again that the primed variables sign the properties of a Monte Carlo particle before its interaction, so, compared to the forward energy formula (26), the variables E and E 0 are interchanged. In Eq. (28) for energies E 0 P E0/2 large scattering angles with l just above the singularity at lmin lead to very large energies and l < lmin would formally result in negative energies. So now, instead of the scattering energies the polar scattering angle cosine is restricted to lmin < l 6 1 with the energy dependence ( 1 E0 < E20 lmin ¼ ; ð29Þ 1 EE00 E0 P E20 so that for l ! lmin, l > lmin the energy E diverges to infinity for E 0 P E0/2, see Fig. 1. Accordingly, the integration of the adjoint differential cross section
330
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
E0 ¼
Emax ; 1 þ 2 Emax =E0
ð32Þ
and decreases to zero at E = Emax as shown in Fig. 2. According to these equations, the adjoint data for incoherent scattering can be calculated to be used in AMOS the same way as the forward data. In addition, the different formulas for energy or scattering angle have to be implemented. 3.2. Pair production
Fig. 1. Relative energy of incoherently scattered photons in dependence of polar scattering angle cosine for several primary photon energies.
RTi ðr; E0 ;
l
Þ ¼ qa ðrÞ pr2e
E E0 þ 1 þ l2 E0 E
Sðr; E; lÞ; ð30Þ
including the scattering function S(r,E,l) and E = E(E 0 ,l) from (28), would result in an infinite adjoint cross section RTi ðr; E0 Þ. Therefore, a maximum energy level has to be set above which scattering is not possible. This is also suitable considering that in physical problems there is always a highest source energy above which no particles are emitted and hence adjoint scattering into these energies would be irrelevant. Finally and most important, an upper energy boundary is needed to ensure convergence by stopping an otherwise endless adjoint calculation. Thus a problem dependent maximum energy Emax is set, to which there is a minimum scattering angle cosine ( 1 E0 < 1þ2EEmax max =E0 0 lmin ðE ; Emax Þ ¼ : ð31Þ 0 1 þ EEmax EE00 E0 P E20 Now the adjoint cross section may be calculated. Depending on the maximum scattering energy it has a point of discontinuity in the first derivation at
Fig. 2. Comparison of forward and adjoint incoherent scattering cross sections (germanium) with two different maximum energies.
In pure photon transport pair production at the nucleus is treated as a scattering process from energies above 2E0 into E0. The two outcome photons are diametrically emitted. Adjoining this is not possible, because this would mean that two photons, each of the energy E0, meet at a point which they approach from exactly opposite directions. Instead, the forward process can be described as single photon scattering, still into the energy E0 but with no direction bound. For this process the forward differential production cross section is Rpð1Þ ðr; E0 ! E; X0 ! XÞ 1 ð33Þ 4p with the common forward cross section Rp(r, E 0 ) of pair production, the term 1/4p accounting for the isotropic emission and the Dirac d-distribution forcing the outcome energy to equal E0. The increase of the cross section by the factor 2 associates to the two annihilation photons. Their emission is not simulated in correlation to each other. Observing Eq. (33) in the adjoint view leads to the same problem as the adjoint incoherent scattering. A maximum energy has to be set, which cuts the differential cross section ¼ 2 Rp ðr; E0 Þ dðE E0 Þ
RTp ðr; E E0 ; X X0 ; Emax Þ 1 2 Rp ðr; EÞ dðE0 E0 Þ 4p ¼ 0
E < Emax ; E P Emax
ð34Þ
so that the integral is limited. In the simulation what is required is that adjunctions, the adjoint Monte Carlo particles, with energy E0, are available to enter adjoint pair production. This is accomplished by using a so-called energy point detector through which particles are incoherently scattered into E0. In contrast to the usual point detector, the outcome energy instead of a space point is the desired parameter. So, incoherent scattering is compelled by multiplying the conditional probability for adjoint incoherent scattering P Ti from Eq. (24). Then, if possible, the particle is scattered into E = E0. The scattering angle cosine l is already set by Eq. (28) and the pdf fiT at l has to be multiplied to the particle’s weight. The azimuthal scattering angle is chosen out of 1/2p. After the following transport to the next interaction point, the particle may undergo coherent scattering. If this interaction process is drawn a further transport follows and again the possibility of coherent
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
scattering. However, if incoherent scattering is drawn, this history ends. But at each of these interaction points, a subhistory is split off and forced to succumb pair production. Here the conditional probability, but this time for adjoint pair production, P Tp is multiplied to the weight. The outcome energy and direction are chosen via fpT . Such sub-histories with the production of adjoint annihilation photons are split off at every interaction point and run parallel to the original Monte Carlo particle history. As the scattering by adjoint pair production leads to high energies above 2E0 and as incoherent scattering and even more coherent scattering at these energies are very improbable, these sub-histories are destined to be very short. Hence the computational effort, though not minimised, is still kept within limits. This energy point detector is also used for realising discrete energy sources where in the adjoint simulation the same problem arises. Before every collision a sub-history is split off and, if possible, incoherently scattered into a discrete source energy. If its following interaction point is inside the source volume it is counted, or else this pseudo particle is abandoned. This is not a particularly efficient procedure, but regarding the special geometries at which adjoint transport is used and thereby its advantage over the forward transport this drawback is tolerable. 4. Example To prove the applicability of AMOS Pt a simple airover-ground problem is calculated by both AMOS and AMOS Pt and the results are compared. Requested is the attenuation by soil cover of varying thickness d on contaminated ground. The contamination is uranium ore with an activity density of 1 Bq/g of its reference nuclide 226Ra. As all following nuclides of the uranium decay series are treated the highest source energy is 2448 keV emitted from 214 Bi. Hence pair production must be considered. The reference point of the detector is located 1 m above the surface. Only by taking advantage of the homogeneous geometry the detector volume may be expanded significantly, so a
331
Fig. 3. Comparison of forward and adjoint results of the air-over-ground problem.
forward simulation may be performed with tolerable statistical error. The resulting ambient dose equivalent rates H_ ð10Þ are compared in Fig. 3, in which the absolute agreement within the statistical errors can be observed. For example the values without soil cover are forward : H_ ð10Þðd ¼ 0 cmÞ ¼ ð506:2 9:4Þ nSv=h; adjoint : H_ ð10Þðd ¼ 0 cmÞ ¼ ð506:8 7:2Þ nSv=h: Besides the verification of the reliability of AMOS Pt, one advantage of the adjoint simulation can also be demonstrated through this example. Via the adjoint calculation it is possible to generate precise energy spectra. As an example, Fig. 4 shows the spectral ambient dose equivalent rate for again d = 0 cm. Some significant source energies are signed out. A comparison with the analogue forward calculated spectrum would bring no gain because of its large uncertainty for high energies. In contrast, the adjoint simulation yields good statistics for large energies, degrading slightly to lower energies. As can be observed in Fig. 4, it has to be noted that only the radiation of the uranium ore is simulated. Nuclides such as potassium or thorium, which must not be disregarded to simulate reality, are not taken into account here.
Fig. 4. Adjoint simulated spectrum of the ambient dose equivalent rate for d = 0 cm.
332
D. Gabler et al. / Nucl. Instr. and Meth. in Phys. Res. B 251 (2006) 326–332
5. Conclusions and discussion In this paper the realisation of adjoint photon transport in AMOS Pt is presented. Its validity is shown through comparison with forward simulated results. The use of the energy point detector firstly leads to a loss of simulation speed. But still the adjoint simulation is clearly preferable at geometries with large physical source and very small physical detector volume and energy range. The energy point detector acts only in one dimension, in which moreover, because of the down- or respectively upscattering condition, one direction is featured. Instead, in a forward simulation a 3D space point detector would have to be used to score in the small detector volume. The loss of speed in the adjoint simulation is therefore not nearly as big as it would be in the forward analogon. As soon as the symmetry is broken forward Monte Carlo transport in air-over-ground problems cannot be carried out as in the here presented example. Acknowledgement We would like to thank David Legrady for his constructive comments on the manuscript. References [1] I. Lux, L. Koblinger, Monte Carlo Particle Transport Methods: Neutron and Photon Calculations, CRC Press, Boca Raton, Florida, 1991. [2] R.R. Coveyou, V.R. Cain, K.J. Yost, Adjoint and importance in Monte Carlo application, Nucl. Sci. Eng. 27 (1967) 219. [3] D.C. Irving, The adjoint Boltzmann equation and its simulation by Monte Carlo, Nucl. Eng. Des. 14 (1971) 273.
[4] M.B. Emmett, The MORSE Monte Carlo Radiation Transport Code System, ORNL-4972, UC-32 – Mathematics and Computers, Oak Ridge National Laboratory, Oak Ridge, Tennessee, 1975. [5] T.E. Booth, F.B. Brown, J.S. Bull, L.J. Cox, R.A. Forster, J.T. Goorley, H.G. Hughes, R.D. Mosteller, R.E. Prael, E.C. Selcow, A. Sood, J.E. Sweezy, R.F. Barrett, S.E. Post, T.L. Roberts, MCNP – A General Monte Carlo N-Particle Transport Code, Version 5, Volume I: Overview and Theory, Los Alamos National Laboratory Report LA-UR-03-1987, April 2003. [6] M.H. Kalos, Monte Carlo integration of the adjoint gamma-ray transport equation, Nucl. Sci. Eng. 33 (1968) 284. ¨ zisik, Y. Yener, On the inverse problem for [7] C.E. Siewert, M.N. O multigroup neutron transport theory, Tech. Notes (rev. 1976) 95. [8] A. De Matteis, R. Simonini, A new approach to the adjoint Boltzmann equation, Nucl. Sci. Eng. 65 (1978) 93. [9] A. De Matteis, R. Simonini, A Monte Carlo biasing scheme for adjoint photon transport, Nucl. Sci. Eng. 67 (1978) 309. [10] S.N. Cramer, Adjoint gamma-ray estimation to the surface of a cylinder – Analysis of a remote reprocessing facility, Nucl. Sci. Eng. 79 (1981) 417. [11] R.W. Roussin, F.A. Schmidt, Adjoint calculations of coupled neutron and gamma-ray transport through concrete slabs, Nucl. Eng. Des. 15 (1971) 319. [12] J.E. Hoogenboom, Adjoint Monte Carlo photon transport in continuous energy mode with discrete photons from annihilation, in: PHYSOR, 2000. [13] J.E. Hoogenboom, Methodology of continuous energy adjoint Monte Carlo for neutron, photon, and coupled neutron-photon transport, Nucl. Sci. Eng. 143 (2003) 99. [14] N.M. Borisov, M.P. Panin, Generalized particle concept for adjoint Monte Carlo calculations of coupled gamma-ray-electron-positron transport, Nucl. Sci. Eng. 150 (2005) 284. [15] J. Henniger, Grundlagen der Abschirmberechnungen mit MonteCarlo-Methoden, Wiss. Z. TU Dresden 39 (3) (1990) 9. [16] D.E. Cullen, J.H. Hubbell, L. Kissel, EPDL97; Lawrence Livermore National Laboratory, Livermore, CA, UCRL50400, Vol. 6, Rev. 5, 1997.