Health Phzsics Pergamon Press 1976. Vol. 30 (February), pp. 167-171. Printed in Northern Ireland
MONTE CARL0 SIMULATED IRRADIATION IN A SPATIALLY DEPENDENT FIELD* T. D. JONES and L. T. DILLMAN? Health Physics Division, Oak Ridge National Laboratory, Oak Ridge, T N 37830
(Received 2 May 1975; accepted 29 JUIJ1975)
Abstract-For Monte Carlo simulated irradiation of any detector having a non-vanishing volume, immersed in a spatially dependent field of radiation, points of incidence of irradiating particles may be determined and then be either accepted or rejected according to the evaluation of a probabilistic model. The model could be based on either the differential arc method or the projected area method discussed in this paper. These two sampling algorithms were developed in order to precisely approximate experimental irradiation conditions for purposes of computer simulated irradiation. I n addition, direction cosines describing the paths of travel of particles incident from any spatially dependent source field may be assigned according to a sampling algorithm applicable to any analytically defined isotimic surface where the partial derivatives exist and are continuous. INTRODUCTION
FOR MONTECARLOsimulated irradiation of a detector1 immersed in a spatially dependent field of radiation, the points of entrance of the irradiating particles may be selected uniformly and randomly according to differential units of surface area. Each position of entrance should then be either accepted or rejected (i.e. the “fate” of a particle) according to the comparison of the magnitude of a random number, chosen between zero and unity, with an evaluation of a probabilistic model, derived from the angular fluence intensities of the irradiating source field about the point of interest. A method of “particle current to fluence” conversion for instruments of detection or for Monte Carlo simulated irradiation was published previously (JONES,1975) and is applicable to radiation source field analogues discussed in this paper. I n fact, this paper is a straightforward extension of the “current to fluence” treatment. A probabilistic model of the angular fluence intensities, P{O, d}, should be
* Research sponsored by the Energy Research and Development Administration under contract with Union Carbide Corporation. t Consultant, Ohio Wesleyan University. $ I n this paper, detector indicates any irradiated surface and can mean a dosimeter or a bioorganism but not a point of detection as is often used in neutron physics.
selected so that the relative responses for different directions are linear in magnitude and, for economy, a value of unity indicates the direction of the maximum irradiating fluence. For a n analytically defined isotimic surface, i.e. a surface defined byf ( X , Y, Z ) equals a constant, it is often convenient to select points of entrance of the irradiating particles randomly but nonuniformly with either surface area or angular fluence intensity and to additionally modify the probabilistic model in order to take this biasing condition into account. This modification of P{O, $1 by M{B, yields a new density function
+>
w, 41 = M@, dl w, 41
(1)
and describes the expectation probability of a particle’s acceptance. Two algorithms which may be used to obtain M{B, 4) distributions will be presented later in this paper. I n addition, a general algorithm for the determination of direction cosines describing the paths of travel of incident particles will be discussed. Although the methods of JONES (1975) and those presented in this paper are completely general, boundaries defined by planes or quadratic surfaces (where cylindrical or spherical symmetry is present) usually simplify the mathematical and computer application of these techniques, because many terms of the equations become zeros. 167
168
IRRADATION IN A SPATIALLY DEPENDENT FIELD
sothat
SAMPLING ALGORITHM BASED ON THE DIFFERENTIAL ARC METHOD
- P d A 4 B 4 + ( B 2 - A 2 ) 2 P 4 cos2 Asin2 I A density function M{O, t$}, which helps to A2B2 describe the “fate” [as discussed in the Intro- dA duction) of a point, may be determined from a (5) comparison of differential arc Iengths in exposure or situations where irradiated areas are determined P3 dS/dl = d B 4 cos2 I A4 sin2 A . (6) by the magnitudes of the arc lengths. Such is the A2B2 case for most surfaces used to approximate a bioorganism for computer simulations and with I n order to obtain the density function for the almost all sensitive volumes of instruments of probabilistic schema, it is necessary to compute detection. dS(4 For a n arbitrarily shaped curve described in dA polar coordinates by S(p, A), the differential arc M{R} = ( 7) length, dS, relative to a differential change in the max angle of rotation, dA, may be expressed as for each and every point of entrance. (2) From equation ( 6 ) dA AB(A2 - B 2 ) d 2S as shown in Fig. 1. = d B 4 cos2 1 A4 sin2 3, I\P (sin 1 cos /I X ( A 2sin2 i. + B 2 cos2 x {(A2 B2)( A 2sin23, B 2 cos2 A) - 3(B4 cos2A + A4 sin2 A)}
+
ps)]
($T
dS++
+
I
c
+
I
+
(8)
and zeros of equation (8) occur when
FIG. 1. Differential arc length of a plane curve.
sin 3, = 0 cos 3, = 0
which implies 1, which implies A 2
= =
0, T, (9) &7~/2,
(10) The normalization of the M(3,) distribution requires the determination of all roots of d2S/d12 and an investigation as to which root yields the maximum value of dS/dA. This rnax [dS(l)/dA] then becomes the normalizing factor of the probabilistic density function denoted by M{A}. I n order to obtain d2S/dA2and then evaluate the max [dS(A)d/A], it is necessary to express p as a function of R and to form dp/dA. This technique is illustrated most profitably for the practical case of a n elliptical cylinder where A2B2 (3) pa = ~2 sin2 1 ~2 cOs2 I
+
and dp/dA
=
p 3 ( B 2- A2) sin A cos 3, A2B2
and (A2 + B2) (A2 sin2 A + B2 C O S ~3,) - 3(B4 cos2 1 A 4 sin2 3,) = 0. If
+
f= B/A
(11)
(12)
then roots of equation (1 1) are found when
A,
= arc sin
If/
1 -2f2 2(1 - f 4 )
and
)
+
dS(1,) Afdsin2 A, f cos2 A, dA (sina I , f cos2 A,)3/2
+
(13),
(14)
which is then compared with A and B from roots (9) and (10) in order to select the rnax [dS(A)/ dA]. Equation (13) breaks down for a narrow range of loas illustrated in Fig. 2 and in this (4) region the normalizing factor becomes either A
169
T. D. JONES and L. T. DILLMAN I .4
1.2
.C
0.a - N
L
O :
0.6
c .-
n
0.4
0.2
C
I
2
3
4
5 f =B/A
6
7
8
9
10-
FIG. 2. Differential elliptical arc lengths for a probabilistic model of detector incidence parameters.
or B, depending on which is larger. This method of obtaining the normalizing factor is similar conceptually to that used by DILLMAN (1974) in the infinite isotropic cloud calculations although the results and application are different. SAMPLING ALGORITHM BASED ON THE PROJECTED AREA METHOD
The preceding probability density function modification is usually applied to cases where the differential units of surface are proportional to differential arc lengths and, in these cases, the points of entry are usually predicted easily from a systematic determination of the upper limits of integration of
*-A
s””
X,.P,.Z, p‘’zu
da(X, Y, 2) =
area of the detector, and da is a differential unit of surface area (JONES et al., 1975). Unfortunately, equation (15) is sometimes difficult to solve for the upper limits of integration as in situations where the irradiatedsurface is not “well-defined” by differential arc lengths; hence, one must resort to other means. One such alternate technique is to uniformly select a random distribution of points on another surface, project these points onto the surface exposed to the irradiating particle so that there is an isomorphism of points between the two surfaces, and then play the Monte Carlo game with a probability density function. An example of the application of this principle is the determination of points of entrance uniformly and randomly on the lateral surface area of a frustrum of a right elliptical cone defined by
net interval zero to unity, A is the total surface Equation (16) forA
= 20, B =
10, C = IOOand
170
IRRADATION IN A SPATIALLY DEPENDENT FIELD
-80 5 2 0 is used to simulate the combined two legs of a person for the reference man phantom (MIRD, 1969) which is used for computer experiments involving spatially dependent radiation fields. To sample uniformly with surface area, one is tempted to choose Z from a density function obtained from the ratios of perimeters in a plane parallel to 2 = 0, i.e. 7
4
I n order to apply this technique, it is necessary to determine the minimum value of cos A so that
and the density function for the probabilistic model is then obtained from (cos A) min
MONTE CARL0 SIMULATED IRRADIATION IN A SPATIALLY DEPENDENT FIELD* T. D. JONES and L. T. DILLMAN? Health Physics Division, Oak Ridge National Laboratory, Oak Ridge, T N 37830
(Received 2 May 1975; accepted 29 JUIJ1975)
Abstract-For Monte Carlo simulated irradiation of any detector having a non-vanishing volume, immersed in a spatially dependent field of radiation, points of incidence of irradiating particles may be determined and then be either accepted or rejected according to the evaluation of a probabilistic model. The model could be based on either the differential arc method or the projected area method discussed in this paper. These two sampling algorithms were developed in order to precisely approximate experimental irradiation conditions for purposes of computer simulated irradiation. I n addition, direction cosines describing the paths of travel of particles incident from any spatially dependent source field may be assigned according to a sampling algorithm applicable to any analytically defined isotimic surface where the partial derivatives exist and are continuous. INTRODUCTION
FOR MONTECARLOsimulated irradiation of a detector1 immersed in a spatially dependent field of radiation, the points of entrance of the irradiating particles may be selected uniformly and randomly according to differential units of surface area. Each position of entrance should then be either accepted or rejected (i.e. the “fate” of a particle) according to the comparison of the magnitude of a random number, chosen between zero and unity, with an evaluation of a probabilistic model, derived from the angular fluence intensities of the irradiating source field about the point of interest. A method of “particle current to fluence” conversion for instruments of detection or for Monte Carlo simulated irradiation was published previously (JONES,1975) and is applicable to radiation source field analogues discussed in this paper. I n fact, this paper is a straightforward extension of the “current to fluence” treatment. A probabilistic model of the angular fluence intensities, P{O, d}, should be
* Research sponsored by the Energy Research and Development Administration under contract with Union Carbide Corporation. t Consultant, Ohio Wesleyan University. $ I n this paper, detector indicates any irradiated surface and can mean a dosimeter or a bioorganism but not a point of detection as is often used in neutron physics.
selected so that the relative responses for different directions are linear in magnitude and, for economy, a value of unity indicates the direction of the maximum irradiating fluence. For a n analytically defined isotimic surface, i.e. a surface defined byf ( X , Y, Z ) equals a constant, it is often convenient to select points of entrance of the irradiating particles randomly but nonuniformly with either surface area or angular fluence intensity and to additionally modify the probabilistic model in order to take this biasing condition into account. This modification of P{O, $1 by M{B, yields a new density function
+>
w, 41 = M@, dl w, 41
(1)
and describes the expectation probability of a particle’s acceptance. Two algorithms which may be used to obtain M{B, 4) distributions will be presented later in this paper. I n addition, a general algorithm for the determination of direction cosines describing the paths of travel of incident particles will be discussed. Although the methods of JONES (1975) and those presented in this paper are completely general, boundaries defined by planes or quadratic surfaces (where cylindrical or spherical symmetry is present) usually simplify the mathematical and computer application of these techniques, because many terms of the equations become zeros. 167
168
IRRADATION IN A SPATIALLY DEPENDENT FIELD
sothat
SAMPLING ALGORITHM BASED ON THE DIFFERENTIAL ARC METHOD
- P d A 4 B 4 + ( B 2 - A 2 ) 2 P 4 cos2 Asin2 I A density function M{O, t$}, which helps to A2B2 describe the “fate” [as discussed in the Intro- dA duction) of a point, may be determined from a (5) comparison of differential arc Iengths in exposure or situations where irradiated areas are determined P3 dS/dl = d B 4 cos2 I A4 sin2 A . (6) by the magnitudes of the arc lengths. Such is the A2B2 case for most surfaces used to approximate a bioorganism for computer simulations and with I n order to obtain the density function for the almost all sensitive volumes of instruments of probabilistic schema, it is necessary to compute detection. dS(4 For a n arbitrarily shaped curve described in dA polar coordinates by S(p, A), the differential arc M{R} = ( 7) length, dS, relative to a differential change in the max angle of rotation, dA, may be expressed as for each and every point of entrance. (2) From equation ( 6 ) dA AB(A2 - B 2 ) d 2S as shown in Fig. 1. = d B 4 cos2 1 A4 sin2 3, I\P (sin 1 cos /I X ( A 2sin2 i. + B 2 cos2 x {(A2 B2)( A 2sin23, B 2 cos2 A) - 3(B4 cos2A + A4 sin2 A)}
+
ps)]
($T
dS++
+
I
c
+
I
+
(8)
and zeros of equation (8) occur when
FIG. 1. Differential arc length of a plane curve.
sin 3, = 0 cos 3, = 0
which implies 1, which implies A 2
= =
0, T, (9) &7~/2,
(10) The normalization of the M(3,) distribution requires the determination of all roots of d2S/d12 and an investigation as to which root yields the maximum value of dS/dA. This rnax [dS(l)/dA] then becomes the normalizing factor of the probabilistic density function denoted by M{A}. I n order to obtain d2S/dA2and then evaluate the max [dS(A)d/A], it is necessary to express p as a function of R and to form dp/dA. This technique is illustrated most profitably for the practical case of a n elliptical cylinder where A2B2 (3) pa = ~2 sin2 1 ~2 cOs2 I
+
and dp/dA
=
p 3 ( B 2- A2) sin A cos 3, A2B2
and (A2 + B2) (A2 sin2 A + B2 C O S ~3,) - 3(B4 cos2 1 A 4 sin2 3,) = 0. If
+
f= B/A
(11)
(12)
then roots of equation (1 1) are found when
A,
= arc sin
If/
1 -2f2 2(1 - f 4 )
and
)
+
dS(1,) Afdsin2 A, f cos2 A, dA (sina I , f cos2 A,)3/2
+
(13),
(14)
which is then compared with A and B from roots (9) and (10) in order to select the rnax [dS(A)/ dA]. Equation (13) breaks down for a narrow range of loas illustrated in Fig. 2 and in this (4) region the normalizing factor becomes either A
169
T. D. JONES and L. T. DILLMAN I .4
1.2
.C
0.a - N
L
O :
0.6
c .-
n
0.4
0.2
C
I
2
3
4
5 f =B/A
6
7
8
9
10-
FIG. 2. Differential elliptical arc lengths for a probabilistic model of detector incidence parameters.
or B, depending on which is larger. This method of obtaining the normalizing factor is similar conceptually to that used by DILLMAN (1974) in the infinite isotropic cloud calculations although the results and application are different. SAMPLING ALGORITHM BASED ON THE PROJECTED AREA METHOD
The preceding probability density function modification is usually applied to cases where the differential units of surface are proportional to differential arc lengths and, in these cases, the points of entry are usually predicted easily from a systematic determination of the upper limits of integration of
*-A
s””
X,.P,.Z, p‘’zu
da(X, Y, 2) =
area of the detector, and da is a differential unit of surface area (JONES et al., 1975). Unfortunately, equation (15) is sometimes difficult to solve for the upper limits of integration as in situations where the irradiatedsurface is not “well-defined” by differential arc lengths; hence, one must resort to other means. One such alternate technique is to uniformly select a random distribution of points on another surface, project these points onto the surface exposed to the irradiating particle so that there is an isomorphism of points between the two surfaces, and then play the Monte Carlo game with a probability density function. An example of the application of this principle is the determination of points of entrance uniformly and randomly on the lateral surface area of a frustrum of a right elliptical cone defined by
net interval zero to unity, A is the total surface Equation (16) forA
= 20, B =
10, C = IOOand
170
IRRADATION IN A SPATIALLY DEPENDENT FIELD
-80 5 2 0 is used to simulate the combined two legs of a person for the reference man phantom (MIRD, 1969) which is used for computer experiments involving spatially dependent radiation fields. To sample uniformly with surface area, one is tempted to choose Z from a density function obtained from the ratios of perimeters in a plane parallel to 2 = 0, i.e. 7
4
I n order to apply this technique, it is necessary to determine the minimum value of cos A so that
and the density function for the probabilistic model is then obtained from (cos A) min
