- Email: [email protected]
Monte Carlo calculation of the average solid angle subtended by a parallelepiped detector from a distributed source
436
Nuclear Instruments and Methods in Physics Research 226 (1984) 436-448 North-Holland, Amsterdam
M O N T E C A R L O C A L C U L A T I O N O F T H E A V E R A G E S O L I D A N G L E S U B T E N D E D BY A PARALLELEPIPED DETECTOR FROM A DISTRIBUTED SOURCE * Lucian W I E L O P O L S K I
Brookhaven National Laboratory, Medical Research Center, Upton, New York 11973, USA Received 1 February 1984
A Monte Carlo method utilizing total variance reduction is applied to the calculation of the solid angle subtended by a parallelepiped detector from an arbitrarily positioned point. Although there is an analytical expression for the solid angle subtended by a rectangle from an arbitrary point, the method should prove valuable when a mean solid angle from a distributed source is required. In addition, the algorithm provides an efficient way to sample the directions of gamma rays which intercept a parallelepiped detector, required in the Monte Carlo calculations of the efficiency of square detectors.
1. Introduction
Knowledge of the solid angle subtended by a detector is required in a variety of problems involving the measurement of nuclear radiation. The general definition of a solid angle subtended by an object at a point P whose position vector is ra is: ~2p = f n ._(_r-- r~) a s ,
6
(1)
I t - r~l 3
where r is the variable position vector of the surface element ds "visible" at P, n is the unit vector normal to ds pointing away from P (the angle between n and r - rp is always less than or equal to ~r/2) and s is the surface over which the integration is to be carried out. Eq. (1) cannot be solved analytically but for the simple cases. W h e n a parallelepiped is viewed from point P the solid angle can be calculated analytically by summing up the partial solid angles subtended by the surfaces of the parallelepiped as they are seen from the point P. Solid angle subtended by a rectangular shape is given by an inverse trigonometric function [1],
~p=arctan((x2--Xp)(y2--yp)/Zp((X2--Xp)2+(y2--yp)
2 +z2)
--arctan{(xl--xp)(Y2--YP)/zP((xl-XP)2+(y2-yP)2+
_arctan{(x2_Xp)(y 1 _ y p ) / / g p ( ( X 2 _ X p ) 2 + ( y l + a r c t a n { ( x I - X p ) ( y 1 -yp)/zp((x
1-Xp)2
_yp)2
]/2} 2 \1/2)
zP)
)
..]_ z 2 ) l / 2 }
+ (yl-yp)2+z2)'/2),
(2)
where the rectangle is enclosed b y the four straight lines x = x I , x = x 2, y = yl, and y = y2(x~ < x 2, yx < Y2). The coordinates of the observation point P are (xp, yp, zp). However, when a distributed source is considered eq. (1) has to be integrated according to the source distribution. The additional integration, in most cases, precludes the possibility of deriving an analytical * Research carried out under the auspices of the Department of energy under contract no. DE-AC02-76CH00016.
0168-9002/84/$03.00 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
L. Wielopolski/ Monte Carlo calculation of average solid angle
437
solution and one has to use a proper quadrature scheme [2-4]. An efficient Monte Carlo algorithm, based on the one developed for a right circular cylinder [5], is applied here for a parallelepiped detector. In this algorithm each history (run) is a success (hit) with a proper weight and the estimate of the solid angle is expressed as the mean weight of all the histories. In addition, the algorithm provides the directional information of the g a m m a radiation required to evaluate the efficiency of parallelepiped detectors.
2. Solid angle subtended by a parailelepiped at a point P It has been shown in ref. [5] that isotropic emission from an arbitrary point source, P, into a unit sphere can be sampled from two marginal density distribution functions for 0 and a O
p(O)dO=l/2sinOdO,
(3)
and p ( a ) d a = 1 / 2 7rda,
0 _< a _< 2~r,
(4)
where 0 and a are the azimuthal and the horizontal angles respectively. To consider a parallelepiped detector, set the origin of a Cartesian coordinate system at one of the corners of the base of the parallelepiped and let the x, y, and z axes coincide with the segments a, b, and c of the parallelepiped (fig. 1). The points P1, P2, and P3 outside the parallelepiped in fig. 1, represent the cases for which only one, two, or three faces of the parallelpiped can be seen, respectively. Thus, these three points represent any possible configuration of an external point relative to the parallelepiped. In each case a different range of azimuthal (0) and horizontal ( a ) angles needs to be sampled according to eqs. (3) and (4). Each of the points is dealt with separately in the following sections. There are 6 subcases from which only one face of the parallelepiped can be seen, 12 subcases from which only two faces can be seen and 8 subcases from which all three faces can be seen. All these subcases are described in the appendix. The mathematical treatment for all the subcases in each case remains the same, only the coordinate system is rotated according to the subcase under consideration. All the possible configuration are treated in the computer program given in table 3. 2.1. Point P~
The point Pj (Xp, yp, Zp) under consideration satisfies the conditions 0 _< xp _< a, 0 _ c and it is depicted in fig. 2 with the pertinent notation. P3
6j
/
/
D
Y i C
I
b .i
.i_.f
1 ~"
B
/
(
Fig. 1. Parallelepiped seen by three different points P1, P2, P3-
L. Wielopolski / Monte Carlo calculation of average solid angle
438
The normal projection of the point P1 intercepts the top surface F G D E of the parallelepiped at point O'. The segment O ' O " is perpendicular to the x axis and is a reference line from which the horizontal angle a is sampled. After integrating eq. (4) the angle a is given by a = ami n + Ran(o/ma
x
--
Otmin) ,
(5)
where Otmin and ama x are the range over which the horizontal angle a is sampled. Ran is a r a n d o m n u m b e r uniformly distributed between 0 and 1 (0 < Ran < 1) [5]. The weight associated with this particular selection of the angle a, We, is obtained .after integrating eq. (4) [5] Wet = (o/max -- o/rain)/2qr"
(6)
The angles o/max and timi n are adjusted to assure that the selected direction of the ray always intercepts the detector. In addition, four angles o/l, a 2, a3, a 4 (defined between the segment O ' O ' and the segments O ' G , O ' D , O'E, O ' F respectively, see fig. 2) define four regions in which the azimuthal angle 0 is calculated. The angles o/min, Ot. . . . O/l, 0(2, O/3, and 0/4 depend on the position of the point P i ( x p , yp) and are given as follows: (xp,yp)=
(0, 0),
o/min=q'/",
o~1=q7 ,
a 3 = ~r - arctan ( a / b ) , 0£4 =
(Xp, y p ) = ( O , b ) ,
3~r.2, 2rr,
o/min = O/ma x =
(7a)
3¢r/2;
a I = 31r/2, o/2 = 3~r/2, (7b)
a 3 = 31r/2, a 4 = 2~r - a r c t a n ( a / b ) ; (Xp,yp)=(a,O),
O:mi n =
~r/2,
a I = 7r/2, 0/2 = 7/" - -
O~ma x = q'g~
arctan(a/b),
(7c)
O/3 = qT"~
O/4 = qT;
//~moPIx(xp'yp'zp) / / / / /
I
x
A
\\e
//i
i
I
1
I
F i g . 2. T h e c o n f i g u r a t i o n
,
o f PI f r o m w h i c h o n l y o n e f a c e o f t h e p a r a l l e l e p i p e d
c a n b e seen.
439
L. Wielopolski / Monte Carlo calculation of average ~olid angle
(xp, yp) = (a, b),
ami. - O,
a, = a r c t a n ( a / b ) ,
am~ ~ = ~r/2,
a z = ¢r/2,
(7d)
a 3 = ~r/2, a 4 = ~r/2; (xp,yp)=(O,O
0/mi. = ~r, o/max =
a, = rr,
2~r,
0/2 = ~r,
(7e)
0/3 = ~r + a r c t a n [ a / ( b " y p ) ] , 0/4 =
(xp,yp)=(a,O
2.r - a r c t a n ( a / y p ) ;
0/,~i. --- 0,
cq=arctan(a/yp),
0/max=~r,
0/2=Ir-arctan
[/( a b . y)] p ,
(7f)
0/3 ~ '//' ~ Ol 4 =
( x p , y . ) = (0 < a , 0),
q'?;
O~min = ~r/2,
0/, = ~r/2,
a m . x = 3~r/2,
0/2 = ~ r - a r c t a n ( x p / b ) ,
(7g)
0/3 = qr + a r c t a n [ ( a - x e ) / b ] , 0/4 = 3 ~ r / 2 ; ( y p , y p ) = (0 < a , b ) ,
ami. = - ~ r . 2 ,
0/1 "-- a r c t a n ( x e / b ) ,
0~max = ~r/ 2 ,
0/z = ¢r/2,
(7h)
0/3 = 3~r/2, 0/. = 2.r - a r c t a n [ ( a - X p ) / b ] ; ( y p , y p ) = (0 < a , 0 < b ) ,
0/~i. - 0,
0/1 = a r c t a n ( x p . y p ) ,
0/m.~ = 2 ~ r ,
0/2 = ~ r -
arctan[xo/(b-yp)
0/3 = ~r + a r c t a n [ ( a - X p ) / ( b 0/4 = 2~r + a r c t a n [ ( a -
]. -yp)],
(7i)
xo)/yp].
T h e a z i m u t h a l a n g l e 0 is s a m p l e d b e t w e e n 0mi n a n d 0ma~, Omi n = O,
Omax = a r c t a n ( p / ( z p -
c),
(S) (9)
w h e r e p is given b y p = y e / c o s 0/,
0 _< 0 / < 0/1,
p = X p / s i n et,
0/1 <-- a < a2 ,
p ~- ( y v
--
b ) / c o s 0/,
p = (x e - a)/sin p
a,
- - y p / / C O S 0/,
0/2 -< 0/-< 0/3, 0/3 ~ 0/ < 0/4' 0/4 ~ 0/ < 2'/7".
OOa) OOb) (10c) (lOd) (lOe)
I n t e g r a t i o n o f eq. (3) y i e l d s for the a n g l e 0 a n d the w e i g h t W e a s s o c i a t e d w i t h the a n g l e 0 [5] 0 = arcos[cos
0mi n - -
R a n ( c o s 0 ~ . - cos 0max) ] ,
(n)
and ~ V0 = (COS Omi n - - COS Oma x ) / 2 .
(12)
440
L. Wielopolski
/
Monte Carlo calculation of average solid angle
It should be noted that there is no need to sample angle 0, eq. (11), for solid angle calculations. It is provided for completeness when the interception point of the g a m m a ray with the detector is required. The combined weight for this history, W,, is W~ = W, W0,
(13)
and the solid angle I2p is given by
1 N
(14)
where N is the total n u m b e r of histories. The standard deviation in $2p is given by
o~, =
[,
N(N-
1)
(15)
y" W/2 - NQ2 i=1
2.2. Point P2
The point P2 under consideration satisfies the conditions Xe > a, 0 < y e < b and Zp > c (see fig. 3). The horizontal angle a is allowed in the range of Otmax only, where arnax is defined by A O ' C and is given by area X= a 1 + t~=,
(16)
cq = arctan[ y p / ( X p - a ) ] ,
(17)
where
P2 (Xp ,yp, Zp) Omaix/z'~zJlomin //
II
G z /4
/ I
);
/
i1'
" \\\\/ /
iii/
/','l
/I i
,
F~I'" " / i~
l /
A t
" ".P~
'i
It/ /\ / l I,,\
\ ,, "I I x If ,~J
J
///~E / ' //I
I °
Y
!
"<< IC
Fig. 3. The configuration of P2 from which two faces of the parallelepiped can be seen.
L. Wielopolski / Monte Carlo calculation of average solid angle
441
and a2 = a r c t a n [ ( b - y p ) / ( x p
- a)].
(18)
N o w the horizontal angle a can be sampled from a = R a n . a ....
(19)
and its weight is W a = O~max//2qT.
(20)
F o r a point source W~ remains constant for all the following histories. The range of 0 varies for each selection of the angle a. The calculation of the range of 0 is determined by the region into which the selected angle a falls. In fig. 3 the three regions are defined by the angles all and o~22. The radius #2 defined by O T is given by 02 = y p / s i n ( c q
- a),
P2 = Xp/COS(Otl
-- a),
02 = ( b - y p ) / s i n ( a
- a~),
0 < a < a,1,
(21a)
all
(21b)
_~< a < Ot22 ,
a22 < a < area x,
(21C)
where a n = a I - arctan(yp/xp),
(22)
0t22 = OtI d- a r c t a n [ ( b - y p ) / X p ] ,
(23)
and the angle Omax is given by 0max = arctan O 2 / ( Z a - c).
(24)
Angle 0mi. is independent of the selected region and is always given by Omi n =
arctan p l / z p ,
(25)
where Oa = ( x p - a ) / c o s ( aa - a ) .
(26)
The angle 0 can now be sampled according to eq. (11), although there is no need for it when calculating a solid angle. The weight ~ is given by eq. (12). The solid angle I2p and its standard deviation is given by eqs. (14) and (15) respectively. 2.3. Point P3
The point P3 under consideration here satisfies the conditions a < xp, b < yp, and c < z e (see fig. 4). The angle ctmax over which the angle et is allowed to vary is given b y Otmax = ~r//2 -- Otl -- otz ,
(27)
a I = arctan(xp - a)/yp,
(28)
a 2 = arctan(yp - b)xp.
(29)
where
and
442
L. Wielopolski / Monte Carlo calculation of average solid angle
T w o critical angles a]] and 0t22 defined by F O ' G and F O ' E respectively, are given by a,] = a r c t a n ( x p / y p )
-
(30)
a I,
and a22 = arctan[(xp -
a)/(yp
b)] - a ] .
-
(31)
The angle a is sampled according to eq. (19) and the weight factor I,V~ is given by eq. (10). For a given angle a the angles 6max and 0 m i n c a n be calculated as follows: Ümax =
arctan[ yp/COS(a
0ma x =
arctan[xp/COS(ama x - a + a 2 ) ( z P - c ) ] ,
+ al)(Z P - c)],
0 ~ a < all ,
(32a)
or
a,1 < a < a . . . .
(32b)
and 0mi. = a r c t a n [ ( x p -
a)/sin(a
+
a,)zp],
0 < a < a22,
(33a)
a22 < a < ama x .
(33b)
or
0min = a r c t a n [ ( y o - b ) / c o s ( a
+ cq)Zp],
AS before, there is no need to sample angle 0. The weight factor I,Vo is given by eq. (12) and the solid angle 12p and oap are given by eqs. (14) and (15), respectively.
3. Surface and volumetric sources In the previous section an algorithm for the calculation of a solid angle subtended by a parellelepiped at a point was presented. The calculation required the three dimensions of the parallelepiped and the three
e.,o~.~ s
I
~s
,Zp)
iI
_
I
/ s
[
T'
I
I ~
\
\
i
, / i
'/
Fig. 4. The configuration of P3 from which three faces of the parallelepiped can be seen.
L. Wielopolski / Monte Carlo calculation of average solid angle
443
coordinates of a point. For a point source these parameters remain constant. However, when the average solid angle from a surface or a volumetric source is considered, a new point sampled properly from the source distribution must be determined separately for each history. In other words, x e, yp, and z 0 will accept new values for each history. At this stage it is also possible to introduce special weight which will account for nonisotropic emission from the source, or self-attenuation in a volumetric source.
4. Results and discussion The results of. the calculations are summarized in table 1. Results 1 - 5 represent the P1 case. The P2 and P3 cases are represented by results 6, 7, 8, and 9, respectively. Results 10-13 represent a distributed square source concentric and parallel to one of the faces of the parallelepiped. Table 1 includes the point coordinates (in cm), number N of histories used in each run, the solid angle 12p calculated according to eq. (14), and the standard deviation of the solid angle oap calculated according to eq. (15). The percent relative standard deviation and the central processing unit (cpu) time for each run. Since the calculations were performed on two computer systems, namely VAX-11/780 (DEC) and the 6600 (CDC), there are two columns for the solid angle and the cpu time. The last column is the analytical solution according to eq. (2). In the case of P2 the analytical solution involved summation of the solid angle due to two faces and in the case of P3, summation over 3 faces of the parallelepiped. The analytical solution in case 10 is from ref. [2], and for cases 11-13 from ref. [3]. It is apparent from table 1 that the Monte Carlo solutions are in good agreement with the exact analytical solutions, although only 1000 histories were used. Solutions 1 - 4 involve an identical case except that the number of histories increase successively by a factor 10 (table 1). The convergence of the solution is rapid and the standard deviation decreases successively by a factor ~ . The execution time on the C D C computer is faster than on the VAX system. The time increase for larger number of histories increases faster on the C D C system than on the VAX system, until in both systems the time increases proportionally to the increased computational task. The differences are due to different overhead time schemes in the computer systems. In table 1, point pairs (2,5), (6,7) and (8,9) are symmetrical. The standard deviations of the points 6, 7, 8, 9, and 10 are lower than those of points 2 and 5, although they run for the same number of histories. This is attributed to the range of values which the variables accept during the calculations. For example in the case of P1 a point with the coordinateds (xp, yp, zp) -- (10, 5, 10) yields 0.8% relative standard deviation. Point 2 in table 1 which belongs to the same has coordinates (1, 1, 10) and yields 3.77% relative standard deviation. Rowlands [2] reported one value for the solid angle calculation from a distributed source, entry 10 in table 1. The Monte Carlo result was divided by 4 to fit the analytical solution. There may be some ambiguity in ref. [2] with regard to the dimensions of the variables. The numerical approximation of the solid angle from a square distributed source suggested by Gillespie [3], entries 11-13 in table 1, appears to underestimate the solid angle at short distances between the source and the detector. At large distances the solutions match. In the Monte Carlo calculations, as the source approached the detector, the solid angle converged to 0.5 fractional solid angle.
Appendix A few of the problems associated with Monte Carlo calculations are: the large number of histories the programs have to run, maintaining tract of the particles and their tallies. Another problem arises when there is a need to transfer to calculations to a designated place in the program or a subroutine according to the values of the selected variables under considerations. This last problem has been addressed in the present work. It is based on the idea that by assigning a proper numerical value (NV) to each individual case it might be possible to identify a pattern according to
25
25
8
9
15
15
15
15
- 5
10
10
10
10 10
10
10
10
z 1,
100 1000
1000
1000
1000
1000
100000 1000
10000
N
histories
Number of
3.84 × 1 0 -
3 . 8 4 × 10 - 2
8.28×10 -2
8.29×10 -2
1.11×10 -1 1.03×10 -1
1.12×10 -1
1 . 1 6 × 10 -1
1 . 1 1 × 1 0 -1
Solid a n g l e ~2p
Separation
1.25
10.000
11
12
13
1000
1000
1000
1000
7.94×10 -4
3.95 × 1 0 -
2 . 2 2 × 10 - l
2.36×10 -2
4
4
2.64×10 -5
1.17 × 1 0 - 3
4.04×10 -3
3.07 × 10
3.09 × 1 0 - 4
3 . 0 9 × 10 - 4
5.47×10 -4
5.42 × 10
4.34×10 -4 4.18×10 -3
3.32
2.95
1.82
1.30
0.80
0.80
0.66
0.65
0.39 4.07
1.22
3.77
12.20
Relative S.D.
-
-
-
-
2.78
2.75
2.65
2.80
86.84 2.51
9.99
2.66
1.70
CPU time (s)
-
-
-
2.33×10 -2
3.74 × 1 0 - 2
3 . 7 4 x 10 - 2
8.15×10 -2
8.16×10 -2
1.12×10 -1 1.13×10 -1
1 . 1 2 × 1 0 -1
1 . 1 2 × 10 -1
1.01 × 10 -1
Solid a n g l e ~V
M . C . C D C 6600
-
-
-
0.51
0.48
0.48
0.48
0.47
41.34 0.46
4.18
0.46
0.07
CPU time (S)
1
2
7.91×10 -4
3.66 × 1 0 - 2
1.58×10 -2
2.36×10 -2
3.79 × 10
3.79×10 -2
8.24x10 -2
8.24×10 -2
1 . 1 2 × 1 0 -1 1 . 1 2 × 1 0 -1
1.12×10 -1
1.12×10 -I
1.12×10
solution
Analytical
c m , a n d c = 0 cm.
a) T h e d i m e n s i o n s o f the p a r a l l e l e p i p e d f o r the cases 1 - 9 a r e a = 20 c m , b = 10 c m , a n d c = 5 c m . F o r c a s e n u m b e r 10 a = b = 2 c m , a n d c = 0. F o r c a s e n u m b e r s 1 1 - 1 3 a = b = 1
1.5
0.25
10
distance
3
1 . 3 7 x 10 - 3
4.37 x 10
1.36×10 -2
Standard deviation
M.C. V A X - 1 1 / 7 8 0
Distributed square sources I cm by 1 cm parallel a n d concentric to the detector
5
15
7
1 19
6
1
1
3
4 5
1 9
1
1
1
1
yp
2
x r,
Point coordinates (cm)
1 a)
Case
Solid angle b y M o n t e C a r l o a n d a n a l y t i c a l s o l u t i o n
Table 1
t,,
¢
t~
L. Wielopolski / Monte Carlo calculation of average solid angle
445
which a proper transfer can be performed. Thus, saving a considerable a m o u n t of computational effort (multiple I F statements) in the preliminary assignment of the problem. In general if the case depends on n variables and each variable m a y assume m values from 0 to m - 1 then N V is defined as
~a,b'-',
NV=
(34)
i~l
where b is the base (b = m) and the coefficient a i is an integer 0 < ai < m - 1. Eq. (34) is nothing else but a numerical system with base b. In the present work there are three variables x e, yp, and ze. Each variable assumes three values; 0 on the left side of the parallelepiped, 1 facing the parallelepiped, and 2 on the right side of the parallelepiped. Thus, the base equals 3 and N V can be constructed as follows: NV=a
I +3a 2+9a
(35)
3.
With the coordinate system used in the text
a i=O,
xp,yp, zp
(36)
a,=l,
O<_xp,yp, Zp<_a,b,c,
(37)
ai=2,
x p , y p , gp, > a , b , c .
(38)
It is possible now to construct table 2.
Table 2 Numerical values of the P1, P2, and P3 cases zp (a3)
yp (a2)
xp (al)
NV
Number of faces seen
0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2
0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2
0 1 2 3 4 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
3 2 3 2 1 2 3 2 3 2 1 2 1 inside 1 2 1 2 3 2 3 2 1 2 3 2 3
P1 1 face
P2 2 faces
P2 3 faces
4
7 9 10 1i 12 14 !5 16 17 18 19 20 21 22 23 24 25 26
!OO CONTINUE IF(A!.EQ.I.
.AND.
S I N G L E FACE
A2.~Q.I.)T~EN
IF(XP. LT.O.0)THEN A1=0.0 ELSE IF(XP. GT.A)THEN A1=2.0 ELSE A1=1.0 END I F IF(YP. LT.O. 0)THEN A2=0.0 ELSE I F ( Y P . GT.B)THEN A2-2.0 ELSE A2=1.0 E NDI F IF(ZP. L T . O . O ) T H E N A3=0.0 ELSE I F ( Z P . GT.C)THEN A3=2.0 ELSE A3=1.0 END I F ORDER=AI+3.*A2+9.*A3 IINT=2*INT(ORDER/2.) I F ( I I N T . LT. ORDER)THEN GO TO 2 0 0 ELSE IF (A1.EQ.$. . O R . A2. EQ.I. * .OR. A3. EQ. 1.)THEN GO TO 1 0 0 ELSE GO TO 3 0 0 END I F
CASE SELECTION
ZRAN=RAN ( I I ) XP=ZRAN ZRAN=RAN( I I ) YP=Z~AN
POINT SELECTION FROM A IxlCM SOURCE
BUR=@. • SUMSQ=O. 0 TYPE 1000 1 0 0 0 FORMAT (IH$~ 2X, ' INSERT Aw B~ C, XP, YP~ ZP, N = ' ) ACCEPT *~A~B~C, XP, YP~ ZP, N DO 5 0 0 I = I , N
PROGRAM SOLANG PI=3. 141592654 11=987654321 STDERR-O. • SOLAN=O. •
Table 3 Solid angle program in Fort:an V language
AP=XP BP=ZP BB=C AA=A I F ( A 2 . EQoO. ) CP--YP I F ( A 2 . EQ. 2. ) CP=YP-B
) CP=--XP ) CP=XP-A
) C.P=-ZP ) CP-ZP-C .AND. A3. EQ. 1. )THEN
END I F IF(AP. EQ.O.)THEN IF(BP. EQ.O.)THEN ALFMIN=Pl ALFMAX=3.*Pl/2. ALFAI=Pl RLFA2=PI ALFA3=PI+ATAN(AA/BB) RLFA4=3.*Pl/2. ELSE IF(BP. EQ. BB)THEN ALFMIN-3.*PI/2. ALFMAX=2.*Pl ALFAI=3.*PI/2. ALFA2=3.*PI/2. ALFA3=3.*PI/2 ALFA4=2.*PI-ATAN(AA/BB) ELSE ALFMIN=Pl ALFMAX=2.*Pl ALFAI=PI ALFA2=PI ALFA3=PI+ATAN(AA/(BB-BP)) ALFA4=2.*PI-ATAN(AA/BP) END I F E L S E I F ( A P . EQ. AA)THEN IF(BP. EQ.O.)THEN ALFMIN=PII2. ALFMAX=Pl ~LFAI=Pl/2. RLFA2=PI-ATAN(AP/BB) ALFA3=Pl ALFA4=Pl ELSE IF(BP. EQ. BB)THEN ALFMIN=O. ALFMAX=PI/2. QLFAI=ATAN(AP/BB)
ELSE
AP-XP BP=YP BB-B RA=A I F ( A 3 . EQ.O. I F ( A 3 . EQ. 2. ELSE I F ( A 2 . EO. 1. AP=Zp BP=YP BB=B AA-C I F ( A t . EQ. 0. IF(AI.EQ.2.
•
I F ( B P . EO.O.)THEN ALFMIN,=PI/2. ALFMAX=3.*PI/2.
END I F
RLFNINaO. ALFMRX=Pl ALFRI-ATRN(Rp/Bp) AL#I~mPI-ATAN(API(BB-BP)) ALFA3-PI RLFA4=PI
ALFA=ALFMIN+ZRAN*(ALFMAX-ALFMIN) IF(ALFA. L T . O . ) ALFA=2.*PI+ALFA CALFA=COS(ALFA) SALFA=SIN(ALFA) IF(ALFA. GE.Q. .AND. ALFA. LT. ALFAI)THEN RO=BP/CALFA ELSE IF(ALFA. GE. ALFA1 .AND. ALFA.LT. * ALFA2)THEN RO=AP/SALFA ELSE IF(ALFA. GE. ALFA2 .AND. ALFA.LT. * ALFA3)THEN RO=(BP-BB)/CALFA ELSE IF(ALFR. GE. ALFA3 .AND. ALFA.LT. * ALFA4)THEN RO=(AP-AA)/SALFA ELSE RO=BP/CALFA END I F
ZRAN=RAN(II)
ALFAI=Pl/2. ALFR2=PI-ATAN(AP/BB) ALFA3=PI+ATAN((I~-AP)/BB) PI.FA4=3.*Pl/2. ELSE IF(BP, EQ.B~)THEN ALFMIN=-PI/2. ALFMAX=PI/2. ALFAI=ATAN(AP/BB) ALFA2=PI/2. ALFA3-3.*PI/2. ALFA4=2.*PI-ATAN((AA-AP)/BB) ELSE ALFMIN=O. ALFMAX=2. * P I ALFAI=ATAN(AP/Bp) ALFR2=PI-ATAN(AP/
ELSE
ELSE
RLFR~PII2. RLFR3-PII2. RLFA4~PII2.
C C C 'C
DOUBLE FACE
IF(AS. EQ.O.) THEN CP=C-ZP ELSE CD=ZP END IF ELSE P~=Ym IF~A3. E Q . ~ . ) THEN CP=C-ZP ELSE CP=ZD END IF END IF ELSE I F ( A ~ . E Q . I . ) THEN AA=A BB=B CC=C BP=YP IF(A1.EQ.O.) THEN AP=A-XP IF(A3. EO.O.) THEN CP=C-ZP ELSE CP=ZP END I F ELSE AD=XP IF(A3. EQ.O.) THEN CP=C-ZD ELSE Cp=zP END I F END I F ELSE AA=A
THEN AA=B BB=A CC=C BP=XP IF(A2. EQ.O.) T~EN
200 CONTINUE lC(A!.EO.I.)
I!
WTETA=(I.-COS(TETMAX))/2. W=WTETA*WALFA SUM=SUM+W SUMSQ=SUMSQ+W*W CONTINUE GO TO 4 0 0
ZRAN=RAN(II) TETA=ACOS(I.-ZRAN*(!.-COS(TETMAX>))
TETMAX=ATAN(RO/CP)
C C C
C
C E C C C
C
3 0 0 CONTINUE AA-A BB=B CC=C
T R I P L E FACE
WTETA=(EOS(TETMIN)-COS(TETNAX))/2. W=WTETA*WALFA SUM=SUM+W SUMSQ=SUMSQ+W*W 21 CONTINUE GO TO 4 0 0
ZRAN=RAN(II) TETA=ACOS(COS(TETMIN)-ZRAN*(COS(TETHIN) COS(TETMAX)))
BB=C CC=B BP=C-ZP !F(A!,EQ.O.) THEN AP=A-XO IF(A2. EQ.O.) THEN CO=B-yP ELSE CD=YP END I F ELSE AP=XP ~ A 2 . E~.O.) THEN CP=B-YD ELSE CP=YP END IF END IF END IF ALFAI=ATAN(BP/(AP-AA)) ALFA2=ATAN((BB-BP)/(AP-AA)) RLFMAX=ALFAI+ALFA2 WALFA=ALFMAX/2./Pl DO 21 I = I , N ZRAN=RAN(II) ALFA=IRAN*ALFMAX TETMIN=ATAN((AP-AA)/COS(ALFAI-ALFA)/CP) ALFAII=ALFAI-ATAN(BP/AO) ALFA22=ALFAI+ATAN((BB-BP)/AP) IF(ALFA. LT. ALFA11)THEN RO=BP/SIN(ALFA1-ALFA) ELSE IF(ALFA. SE. ALFAll .AND. ALFA.LT° * ALFA22)THEN RO=AP/COS(ALFA1-ALFA) ELSE RO=(BB-BP)/SIN(ALFA-ALFA1) END IF TETMAX-ATAN(RO/(CP-CC))
~ ( B p - ~ ) J~p)
THEN
(CP-CC))
SERR=SQRT(S/N) T Y P E *~DLAN~SERR END
SLIMS@=SUMSQ+W*W C 31 CONTINUE 4 0 0 CONTINUE 5 0 0 CONTI N U E SDLAN-BU~/N S= (SUI~Q-N*SOLAN*SOLAN) / ( N - I )
ELSE TETMIN=ATAN((BP-BB)/COS(ALFA+ALFA1)/CP) END I F WTETA=(COS(TETMIN)-COS(TETMAX))/2. W=WTETA*WALFA SUM=SUM+W
TETNIN=ATRN((AP-AA)ISIN(ALFA+ALFAI)/CP)
*END IF I F ( A L F A . LE. ALFA22)
TETMAX=ATAN(AP/COSCALFMAX-ALFA+ALFA2)/
WALFA=ALFMAXI2./Pl, DO 31 I = I , N : ZRAN=RAN(II) " ALFA-ZRAN*ALFMAX IF(ALFA. LE. ALFA11) THEN TETMAX-ATAN(BP/COSIALFA+ALFA1)/IOP-CCII ELSE
ALFA22-ATAN(~AP-AA')/~BP-BB))-ALFAI
ALFMAX=PII2.@ ~-~FAI-ALFA2 ALFA11=ATAN(AP/BP)~ALFA1
~LFP~=AT~
THE~ CP=C-ZP IF(A2. EO.O.~ THEN BP=B-yP IF(~I.EQ.O.) AP=~q-XP IF(AI.EQ. 2 . ) AP=XP ELSE BP=Yp IF(AI.EQ.@o) AP=A-XP IF(A1.EQ. 2 . ) AP=XP END IF ELSE CP-ZP IF(A2. EQ.O.) THEN BP=B-YP IF(A1.EQ.O.) AP-A-XP IF(AI.EQ. 2 . ) AP=XP ELSE Bp=yp IF(A1.EQ.O.) AP=R-XP IF(A1.EQ. 2 . ) AP=XP END IF END IF ALFAI=ATAN((AA-AA)/BP)
IF(A3.EO.O.)
~.
~
~"
~.
~
448
L. Wielopolski / Monte Carlo calculation of auerage solid angle
It is a p p a r e n t from the table that all the o d d values of N V can be a d d r e s s e d to P2 calculations, except p o i n t 13 which represents a p o i n t inside the p a r a l l e l e p i p e d and is not being c o n s i d e r e d here. It also can be o b s e r v e d that in o r d e r to a d d r e s s a p o i n t for P1 calculations it suffices to assure that any two of the coefficients equal to one, a n d b y default all the o t h e r cases can be a d d r e s s e d to P3, calculations. This p r o c e d u r e is i m p l e m e n t e d in the p r o g r a m in table 2. T h e p r o g r a m in table 3 is written in F O R T R A N V language a n d it i m p l e m e n t s all the features discussed in the article. The p r o g r a m is n o t written as a general c o d e thus, different cases require small changes in the p r o g r a m itself. The p r o g r a m in table 3 represent the case for s a m p l i n g u n i f o r m l y a 1 c m x 1 cm square source. F o r a p o i n t source the four statements following P O I N T S E L E C T I O N F R O M 1 x 1 C M S Q U A R E have to be removed. F o r a specific case of a p o i n t source the calculations can be p e r f o r m e d faster b y r e m o v i n g D O - L O O P 500 a n d i m p l e m e n t i n g D O - L O O P s 11, 21, a n d 31 which are neutralized at present b y letter C in the first column.
References [1] [2] [3] [4] [5]
H. Gotoh and H. Yagi, Nucl. Instr. and Meth. 96 (1971) 485. G. Rowlands, Int. J. Appl. Radiat. Isot. 10 (1061) 86. C.R. Gillespie, Rev. Sci. Instr. 41 (1970) 42. J. Cook, Nucl. Instr. and Meth. 178 (1980) 561. L. Wielopolski, Nucl. Instr. and Meth. 143 (1977) 577.
Nuclear Instruments and Methods in Physics Research 226 (1984) 436-448 North-Holland, Amsterdam
M O N T E C A R L O C A L C U L A T I O N O F T H E A V E R A G E S O L I D A N G L E S U B T E N D E D BY A PARALLELEPIPED DETECTOR FROM A DISTRIBUTED SOURCE * Lucian W I E L O P O L S K I
Brookhaven National Laboratory, Medical Research Center, Upton, New York 11973, USA Received 1 February 1984
A Monte Carlo method utilizing total variance reduction is applied to the calculation of the solid angle subtended by a parallelepiped detector from an arbitrarily positioned point. Although there is an analytical expression for the solid angle subtended by a rectangle from an arbitrary point, the method should prove valuable when a mean solid angle from a distributed source is required. In addition, the algorithm provides an efficient way to sample the directions of gamma rays which intercept a parallelepiped detector, required in the Monte Carlo calculations of the efficiency of square detectors.
1. Introduction
Knowledge of the solid angle subtended by a detector is required in a variety of problems involving the measurement of nuclear radiation. The general definition of a solid angle subtended by an object at a point P whose position vector is ra is: ~2p = f n ._(_r-- r~) a s ,
6
(1)
I t - r~l 3
where r is the variable position vector of the surface element ds "visible" at P, n is the unit vector normal to ds pointing away from P (the angle between n and r - rp is always less than or equal to ~r/2) and s is the surface over which the integration is to be carried out. Eq. (1) cannot be solved analytically but for the simple cases. W h e n a parallelepiped is viewed from point P the solid angle can be calculated analytically by summing up the partial solid angles subtended by the surfaces of the parallelepiped as they are seen from the point P. Solid angle subtended by a rectangular shape is given by an inverse trigonometric function [1],
~p=arctan((x2--Xp)(y2--yp)/Zp((X2--Xp)2+(y2--yp)
2 +z2)
--arctan{(xl--xp)(Y2--YP)/zP((xl-XP)2+(y2-yP)2+
_arctan{(x2_Xp)(y 1 _ y p ) / / g p ( ( X 2 _ X p ) 2 + ( y l + a r c t a n { ( x I - X p ) ( y 1 -yp)/zp((x
1-Xp)2
_yp)2
]/2} 2 \1/2)
zP)
)
..]_ z 2 ) l / 2 }
+ (yl-yp)2+z2)'/2),
(2)
where the rectangle is enclosed b y the four straight lines x = x I , x = x 2, y = yl, and y = y2(x~ < x 2, yx < Y2). The coordinates of the observation point P are (xp, yp, zp). However, when a distributed source is considered eq. (1) has to be integrated according to the source distribution. The additional integration, in most cases, precludes the possibility of deriving an analytical * Research carried out under the auspices of the Department of energy under contract no. DE-AC02-76CH00016.
0168-9002/84/$03.00 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
L. Wielopolski/ Monte Carlo calculation of average solid angle
437
solution and one has to use a proper quadrature scheme [2-4]. An efficient Monte Carlo algorithm, based on the one developed for a right circular cylinder [5], is applied here for a parallelepiped detector. In this algorithm each history (run) is a success (hit) with a proper weight and the estimate of the solid angle is expressed as the mean weight of all the histories. In addition, the algorithm provides the directional information of the g a m m a radiation required to evaluate the efficiency of parallelepiped detectors.
2. Solid angle subtended by a parailelepiped at a point P It has been shown in ref. [5] that isotropic emission from an arbitrary point source, P, into a unit sphere can be sampled from two marginal density distribution functions for 0 and a O
p(O)dO=l/2sinOdO,
(3)
and p ( a ) d a = 1 / 2 7rda,
0 _< a _< 2~r,
(4)
where 0 and a are the azimuthal and the horizontal angles respectively. To consider a parallelepiped detector, set the origin of a Cartesian coordinate system at one of the corners of the base of the parallelepiped and let the x, y, and z axes coincide with the segments a, b, and c of the parallelepiped (fig. 1). The points P1, P2, and P3 outside the parallelepiped in fig. 1, represent the cases for which only one, two, or three faces of the parallelpiped can be seen, respectively. Thus, these three points represent any possible configuration of an external point relative to the parallelepiped. In each case a different range of azimuthal (0) and horizontal ( a ) angles needs to be sampled according to eqs. (3) and (4). Each of the points is dealt with separately in the following sections. There are 6 subcases from which only one face of the parallelepiped can be seen, 12 subcases from which only two faces can be seen and 8 subcases from which all three faces can be seen. All these subcases are described in the appendix. The mathematical treatment for all the subcases in each case remains the same, only the coordinate system is rotated according to the subcase under consideration. All the possible configuration are treated in the computer program given in table 3. 2.1. Point P~
The point Pj (Xp, yp, Zp) under consideration satisfies the conditions 0 _< xp _< a, 0 _
6j
/
/
D
Y i C
I
b .i
.i_.f
1 ~"
B
/
(
Fig. 1. Parallelepiped seen by three different points P1, P2, P3-
L. Wielopolski / Monte Carlo calculation of average solid angle
438
The normal projection of the point P1 intercepts the top surface F G D E of the parallelepiped at point O'. The segment O ' O " is perpendicular to the x axis and is a reference line from which the horizontal angle a is sampled. After integrating eq. (4) the angle a is given by a = ami n + Ran(o/ma
x
--
Otmin) ,
(5)
where Otmin and ama x are the range over which the horizontal angle a is sampled. Ran is a r a n d o m n u m b e r uniformly distributed between 0 and 1 (0 < Ran < 1) [5]. The weight associated with this particular selection of the angle a, We, is obtained .after integrating eq. (4) [5] Wet = (o/max -- o/rain)/2qr"
(6)
The angles o/max and timi n are adjusted to assure that the selected direction of the ray always intercepts the detector. In addition, four angles o/l, a 2, a3, a 4 (defined between the segment O ' O ' and the segments O ' G , O ' D , O'E, O ' F respectively, see fig. 2) define four regions in which the azimuthal angle 0 is calculated. The angles o/min, Ot. . . . O/l, 0(2, O/3, and 0/4 depend on the position of the point P i ( x p , yp) and are given as follows: (xp,yp)=
(0, 0),
o/min=q'/",
o~1=q7 ,
a 3 = ~r - arctan ( a / b ) , 0£4 =
(Xp, y p ) = ( O , b ) ,
3~r.2, 2rr,
o/min = O/ma x =
(7a)
3¢r/2;
a I = 31r/2, o/2 = 3~r/2, (7b)
a 3 = 31r/2, a 4 = 2~r - a r c t a n ( a / b ) ; (Xp,yp)=(a,O),
O:mi n =
~r/2,
a I = 7r/2, 0/2 = 7/" - -
O~ma x = q'g~
arctan(a/b),
(7c)
O/3 = qT"~
O/4 = qT;
//~moPIx(xp'yp'zp) / / / / /
I
x
A
\\e
//i
i
I
1
I
F i g . 2. T h e c o n f i g u r a t i o n
,
o f PI f r o m w h i c h o n l y o n e f a c e o f t h e p a r a l l e l e p i p e d
c a n b e seen.
439
L. Wielopolski / Monte Carlo calculation of average ~olid angle
(xp, yp) = (a, b),
ami. - O,
a, = a r c t a n ( a / b ) ,
am~ ~ = ~r/2,
a z = ¢r/2,
(7d)
a 3 = ~r/2, a 4 = ~r/2; (xp,yp)=(O,O
0/mi. = ~r, o/max =
a, = rr,
2~r,
0/2 = ~r,
(7e)
0/3 = ~r + a r c t a n [ a / ( b " y p ) ] , 0/4 =
(xp,yp)=(a,O
2.r - a r c t a n ( a / y p ) ;
0/,~i. --- 0,
cq=arctan(a/yp),
0/max=~r,
0/2=Ir-arctan
[/( a b . y)] p ,
(7f)
0/3 ~ '//' ~ Ol 4 =
( x p , y . ) = (0 < a , 0),
q'?;
O~min = ~r/2,
0/, = ~r/2,
a m . x = 3~r/2,
0/2 = ~ r - a r c t a n ( x p / b ) ,
(7g)
0/3 = qr + a r c t a n [ ( a - x e ) / b ] , 0/4 = 3 ~ r / 2 ; ( y p , y p ) = (0 < a , b ) ,
ami. = - ~ r . 2 ,
0/1 "-- a r c t a n ( x e / b ) ,
0~max = ~r/ 2 ,
0/z = ¢r/2,
(7h)
0/3 = 3~r/2, 0/. = 2.r - a r c t a n [ ( a - X p ) / b ] ; ( y p , y p ) = (0 < a , 0 < b ) ,
0/~i. - 0,
0/1 = a r c t a n ( x p . y p ) ,
0/m.~ = 2 ~ r ,
0/2 = ~ r -
arctan[xo/(b-yp)
0/3 = ~r + a r c t a n [ ( a - X p ) / ( b 0/4 = 2~r + a r c t a n [ ( a -
]. -yp)],
(7i)
xo)/yp].
T h e a z i m u t h a l a n g l e 0 is s a m p l e d b e t w e e n 0mi n a n d 0ma~, Omi n = O,
Omax = a r c t a n ( p / ( z p -
c),
(S) (9)
w h e r e p is given b y p = y e / c o s 0/,
0 _< 0 / < 0/1,
p = X p / s i n et,
0/1 <-- a < a2 ,
p ~- ( y v
--
b ) / c o s 0/,
p = (x e - a)/sin p
a,
- - y p / / C O S 0/,
0/2 -< 0/-< 0/3, 0/3 ~ 0/ < 0/4' 0/4 ~ 0/ < 2'/7".
OOa) OOb) (10c) (lOd) (lOe)
I n t e g r a t i o n o f eq. (3) y i e l d s for the a n g l e 0 a n d the w e i g h t W e a s s o c i a t e d w i t h the a n g l e 0 [5] 0 = arcos[cos
0mi n - -
R a n ( c o s 0 ~ . - cos 0max) ] ,
(n)
and ~ V0 = (COS Omi n - - COS Oma x ) / 2 .
(12)
440
L. Wielopolski
/
Monte Carlo calculation of average solid angle
It should be noted that there is no need to sample angle 0, eq. (11), for solid angle calculations. It is provided for completeness when the interception point of the g a m m a ray with the detector is required. The combined weight for this history, W,, is W~ = W, W0,
(13)
and the solid angle I2p is given by
1 N
(14)
where N is the total n u m b e r of histories. The standard deviation in $2p is given by
o~, =
[,
N(N-
1)
(15)
y" W/2 - NQ2 i=1
2.2. Point P2
The point P2 under consideration satisfies the conditions Xe > a, 0 < y e < b and Zp > c (see fig. 3). The horizontal angle a is allowed in the range of Otmax only, where arnax is defined by A O ' C and is given by area X= a 1 + t~=,
(16)
cq = arctan[ y p / ( X p - a ) ] ,
(17)
where
P2 (Xp ,yp, Zp) Omaix/z'~zJlomin //
II
G z /4
/ I
);
/
i1'
" \\\\/ /
iii/
/','l
/I i
,
F~I'" " / i~
l /
A t
" ".P~
'i
It/ /\ / l I,,\
\ ,, "I I x If ,~J
J
///~E / ' //I
I °
Y
!
"<< IC
Fig. 3. The configuration of P2 from which two faces of the parallelepiped can be seen.
L. Wielopolski / Monte Carlo calculation of average solid angle
441
and a2 = a r c t a n [ ( b - y p ) / ( x p
- a)].
(18)
N o w the horizontal angle a can be sampled from a = R a n . a ....
(19)
and its weight is W a = O~max//2qT.
(20)
F o r a point source W~ remains constant for all the following histories. The range of 0 varies for each selection of the angle a. The calculation of the range of 0 is determined by the region into which the selected angle a falls. In fig. 3 the three regions are defined by the angles all and o~22. The radius #2 defined by O T is given by 02 = y p / s i n ( c q
- a),
P2 = Xp/COS(Otl
-- a),
02 = ( b - y p ) / s i n ( a
- a~),
0 < a < a,1,
(21a)
all
(21b)
_~< a < Ot22 ,
a22 < a < area x,
(21C)
where a n = a I - arctan(yp/xp),
(22)
0t22 = OtI d- a r c t a n [ ( b - y p ) / X p ] ,
(23)
and the angle Omax is given by 0max = arctan O 2 / ( Z a - c).
(24)
Angle 0mi. is independent of the selected region and is always given by Omi n =
arctan p l / z p ,
(25)
where Oa = ( x p - a ) / c o s ( aa - a ) .
(26)
The angle 0 can now be sampled according to eq. (11), although there is no need for it when calculating a solid angle. The weight ~ is given by eq. (12). The solid angle I2p and its standard deviation is given by eqs. (14) and (15) respectively. 2.3. Point P3
The point P3 under consideration here satisfies the conditions a < xp, b < yp, and c < z e (see fig. 4). The angle ctmax over which the angle et is allowed to vary is given b y Otmax = ~r//2 -- Otl -- otz ,
(27)
a I = arctan(xp - a)/yp,
(28)
a 2 = arctan(yp - b)xp.
(29)
where
and
442
L. Wielopolski / Monte Carlo calculation of average solid angle
T w o critical angles a]] and 0t22 defined by F O ' G and F O ' E respectively, are given by a,] = a r c t a n ( x p / y p )
-
(30)
a I,
and a22 = arctan[(xp -
a)/(yp
b)] - a ] .
-
(31)
The angle a is sampled according to eq. (19) and the weight factor I,V~ is given by eq. (10). For a given angle a the angles 6max and 0 m i n c a n be calculated as follows: Ümax =
arctan[ yp/COS(a
0ma x =
arctan[xp/COS(ama x - a + a 2 ) ( z P - c ) ] ,
+ al)(Z P - c)],
0 ~ a < all ,
(32a)
or
a,1 < a < a . . . .
(32b)
and 0mi. = a r c t a n [ ( x p -
a)/sin(a
+
a,)zp],
0 < a < a22,
(33a)
a22 < a < ama x .
(33b)
or
0min = a r c t a n [ ( y o - b ) / c o s ( a
+ cq)Zp],
AS before, there is no need to sample angle 0. The weight factor I,Vo is given by eq. (12) and the solid angle 12p and oap are given by eqs. (14) and (15), respectively.
3. Surface and volumetric sources In the previous section an algorithm for the calculation of a solid angle subtended by a parellelepiped at a point was presented. The calculation required the three dimensions of the parallelepiped and the three
e.,o~.~ s
I
~s
,Zp)
iI
_
I
/ s
[
T'
I
I ~
\
\
i
, / i
'/
Fig. 4. The configuration of P3 from which three faces of the parallelepiped can be seen.
L. Wielopolski / Monte Carlo calculation of average solid angle
443
coordinates of a point. For a point source these parameters remain constant. However, when the average solid angle from a surface or a volumetric source is considered, a new point sampled properly from the source distribution must be determined separately for each history. In other words, x e, yp, and z 0 will accept new values for each history. At this stage it is also possible to introduce special weight which will account for nonisotropic emission from the source, or self-attenuation in a volumetric source.
4. Results and discussion The results of. the calculations are summarized in table 1. Results 1 - 5 represent the P1 case. The P2 and P3 cases are represented by results 6, 7, 8, and 9, respectively. Results 10-13 represent a distributed square source concentric and parallel to one of the faces of the parallelepiped. Table 1 includes the point coordinates (in cm), number N of histories used in each run, the solid angle 12p calculated according to eq. (14), and the standard deviation of the solid angle oap calculated according to eq. (15). The percent relative standard deviation and the central processing unit (cpu) time for each run. Since the calculations were performed on two computer systems, namely VAX-11/780 (DEC) and the 6600 (CDC), there are two columns for the solid angle and the cpu time. The last column is the analytical solution according to eq. (2). In the case of P2 the analytical solution involved summation of the solid angle due to two faces and in the case of P3, summation over 3 faces of the parallelepiped. The analytical solution in case 10 is from ref. [2], and for cases 11-13 from ref. [3]. It is apparent from table 1 that the Monte Carlo solutions are in good agreement with the exact analytical solutions, although only 1000 histories were used. Solutions 1 - 4 involve an identical case except that the number of histories increase successively by a factor 10 (table 1). The convergence of the solution is rapid and the standard deviation decreases successively by a factor ~ . The execution time on the C D C computer is faster than on the VAX system. The time increase for larger number of histories increases faster on the C D C system than on the VAX system, until in both systems the time increases proportionally to the increased computational task. The differences are due to different overhead time schemes in the computer systems. In table 1, point pairs (2,5), (6,7) and (8,9) are symmetrical. The standard deviations of the points 6, 7, 8, 9, and 10 are lower than those of points 2 and 5, although they run for the same number of histories. This is attributed to the range of values which the variables accept during the calculations. For example in the case of P1 a point with the coordinateds (xp, yp, zp) -- (10, 5, 10) yields 0.8% relative standard deviation. Point 2 in table 1 which belongs to the same has coordinates (1, 1, 10) and yields 3.77% relative standard deviation. Rowlands [2] reported one value for the solid angle calculation from a distributed source, entry 10 in table 1. The Monte Carlo result was divided by 4 to fit the analytical solution. There may be some ambiguity in ref. [2] with regard to the dimensions of the variables. The numerical approximation of the solid angle from a square distributed source suggested by Gillespie [3], entries 11-13 in table 1, appears to underestimate the solid angle at short distances between the source and the detector. At large distances the solutions match. In the Monte Carlo calculations, as the source approached the detector, the solid angle converged to 0.5 fractional solid angle.
Appendix A few of the problems associated with Monte Carlo calculations are: the large number of histories the programs have to run, maintaining tract of the particles and their tallies. Another problem arises when there is a need to transfer to calculations to a designated place in the program or a subroutine according to the values of the selected variables under considerations. This last problem has been addressed in the present work. It is based on the idea that by assigning a proper numerical value (NV) to each individual case it might be possible to identify a pattern according to
25
25
8
9
15
15
15
15
- 5
10
10
10
10 10
10
10
10
z 1,
100 1000
1000
1000
1000
1000
100000 1000
10000
N
histories
Number of
3.84 × 1 0 -
3 . 8 4 × 10 - 2
8.28×10 -2
8.29×10 -2
1.11×10 -1 1.03×10 -1
1.12×10 -1
1 . 1 6 × 10 -1
1 . 1 1 × 1 0 -1
Solid a n g l e ~2p
Separation
1.25
10.000
11
12
13
1000
1000
1000
1000
7.94×10 -4
3.95 × 1 0 -
2 . 2 2 × 10 - l
2.36×10 -2
4
4
2.64×10 -5
1.17 × 1 0 - 3
4.04×10 -3
3.07 × 10
3.09 × 1 0 - 4
3 . 0 9 × 10 - 4
5.47×10 -4
5.42 × 10
4.34×10 -4 4.18×10 -3
3.32
2.95
1.82
1.30
0.80
0.80
0.66
0.65
0.39 4.07
1.22
3.77
12.20
Relative S.D.
-
-
-
-
2.78
2.75
2.65
2.80
86.84 2.51
9.99
2.66
1.70
CPU time (s)
-
-
-
2.33×10 -2
3.74 × 1 0 - 2
3 . 7 4 x 10 - 2
8.15×10 -2
8.16×10 -2
1.12×10 -1 1.13×10 -1
1 . 1 2 × 1 0 -1
1 . 1 2 × 10 -1
1.01 × 10 -1
Solid a n g l e ~V
M . C . C D C 6600
-
-
-
0.51
0.48
0.48
0.48
0.47
41.34 0.46
4.18
0.46
0.07
CPU time (S)
1
2
7.91×10 -4
3.66 × 1 0 - 2
1.58×10 -2
2.36×10 -2
3.79 × 10
3.79×10 -2
8.24x10 -2
8.24×10 -2
1 . 1 2 × 1 0 -1 1 . 1 2 × 1 0 -1
1.12×10 -1
1.12×10 -I
1.12×10
solution
Analytical
c m , a n d c = 0 cm.
a) T h e d i m e n s i o n s o f the p a r a l l e l e p i p e d f o r the cases 1 - 9 a r e a = 20 c m , b = 10 c m , a n d c = 5 c m . F o r c a s e n u m b e r 10 a = b = 2 c m , a n d c = 0. F o r c a s e n u m b e r s 1 1 - 1 3 a = b = 1
1.5
0.25
10
distance
3
1 . 3 7 x 10 - 3
4.37 x 10
1.36×10 -2
Standard deviation
M.C. V A X - 1 1 / 7 8 0
Distributed square sources I cm by 1 cm parallel a n d concentric to the detector
5
15
7
1 19
6
1
1
3
4 5
1 9
1
1
1
1
yp
2
x r,
Point coordinates (cm)
1 a)
Case
Solid angle b y M o n t e C a r l o a n d a n a l y t i c a l s o l u t i o n
Table 1
t,,
¢
t~
L. Wielopolski / Monte Carlo calculation of average solid angle
445
which a proper transfer can be performed. Thus, saving a considerable a m o u n t of computational effort (multiple I F statements) in the preliminary assignment of the problem. In general if the case depends on n variables and each variable m a y assume m values from 0 to m - 1 then N V is defined as
~a,b'-',
NV=
(34)
i~l
where b is the base (b = m) and the coefficient a i is an integer 0 < ai < m - 1. Eq. (34) is nothing else but a numerical system with base b. In the present work there are three variables x e, yp, and ze. Each variable assumes three values; 0 on the left side of the parallelepiped, 1 facing the parallelepiped, and 2 on the right side of the parallelepiped. Thus, the base equals 3 and N V can be constructed as follows: NV=a
I +3a 2+9a
(35)
3.
With the coordinate system used in the text
a i=O,
xp,yp, zp
(36)
a,=l,
O<_xp,yp, Zp<_a,b,c,
(37)
ai=2,
x p , y p , gp, > a , b , c .
(38)
It is possible now to construct table 2.
Table 2 Numerical values of the P1, P2, and P3 cases zp (a3)
yp (a2)
xp (al)
NV
Number of faces seen
0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2
0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2
0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2
0 1 2 3 4 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
3 2 3 2 1 2 3 2 3 2 1 2 1 inside 1 2 1 2 3 2 3 2 1 2 3 2 3
P1 1 face
P2 2 faces
P2 3 faces
4
7 9 10 1i 12 14 !5 16 17 18 19 20 21 22 23 24 25 26
!OO CONTINUE IF(A!.EQ.I.
.AND.
S I N G L E FACE
A2.~Q.I.)T~EN
IF(XP. LT.O.0)THEN A1=0.0 ELSE IF(XP. GT.A)THEN A1=2.0 ELSE A1=1.0 END I F IF(YP. LT.O. 0)THEN A2=0.0 ELSE I F ( Y P . GT.B)THEN A2-2.0 ELSE A2=1.0 E NDI F IF(ZP. L T . O . O ) T H E N A3=0.0 ELSE I F ( Z P . GT.C)THEN A3=2.0 ELSE A3=1.0 END I F ORDER=AI+3.*A2+9.*A3 IINT=2*INT(ORDER/2.) I F ( I I N T . LT. ORDER)THEN GO TO 2 0 0 ELSE IF (A1.EQ.$. . O R . A2. EQ.I. * .OR. A3. EQ. 1.)THEN GO TO 1 0 0 ELSE GO TO 3 0 0 END I F
CASE SELECTION
ZRAN=RAN ( I I ) XP=ZRAN ZRAN=RAN( I I ) YP=Z~AN
POINT SELECTION FROM A IxlCM SOURCE
BUR=@. • SUMSQ=O. 0 TYPE 1000 1 0 0 0 FORMAT (IH$~ 2X, ' INSERT Aw B~ C, XP, YP~ ZP, N = ' ) ACCEPT *~A~B~C, XP, YP~ ZP, N DO 5 0 0 I = I , N
PROGRAM SOLANG PI=3. 141592654 11=987654321 STDERR-O. • SOLAN=O. •
Table 3 Solid angle program in Fort:an V language
AP=XP BP=ZP BB=C AA=A I F ( A 2 . EQoO. ) CP--YP I F ( A 2 . EQ. 2. ) CP=YP-B
) CP=--XP ) CP=XP-A
) C.P=-ZP ) CP-ZP-C .AND. A3. EQ. 1. )THEN
END I F IF(AP. EQ.O.)THEN IF(BP. EQ.O.)THEN ALFMIN=Pl ALFMAX=3.*Pl/2. ALFAI=Pl RLFA2=PI ALFA3=PI+ATAN(AA/BB) RLFA4=3.*Pl/2. ELSE IF(BP. EQ. BB)THEN ALFMIN-3.*PI/2. ALFMAX=2.*Pl ALFAI=3.*PI/2. ALFA2=3.*PI/2. ALFA3=3.*PI/2 ALFA4=2.*PI-ATAN(AA/BB) ELSE ALFMIN=Pl ALFMAX=2.*Pl ALFAI=PI ALFA2=PI ALFA3=PI+ATAN(AA/(BB-BP)) ALFA4=2.*PI-ATAN(AA/BP) END I F E L S E I F ( A P . EQ. AA)THEN IF(BP. EQ.O.)THEN ALFMIN=PII2. ALFMAX=Pl ~LFAI=Pl/2. RLFA2=PI-ATAN(AP/BB) ALFA3=Pl ALFA4=Pl ELSE IF(BP. EQ. BB)THEN ALFMIN=O. ALFMAX=PI/2. QLFAI=ATAN(AP/BB)
ELSE
AP-XP BP=YP BB-B RA=A I F ( A 3 . EQ.O. I F ( A 3 . EQ. 2. ELSE I F ( A 2 . EO. 1. AP=Zp BP=YP BB=B AA-C I F ( A t . EQ. 0. IF(AI.EQ.2.
•
I F ( B P . EO.O.)THEN ALFMIN,=PI/2. ALFMAX=3.*PI/2.
END I F
RLFNINaO. ALFMRX=Pl ALFRI-ATRN(Rp/Bp) AL#I~mPI-ATAN(API(BB-BP)) ALFA3-PI RLFA4=PI
ALFA=ALFMIN+ZRAN*(ALFMAX-ALFMIN) IF(ALFA. L T . O . ) ALFA=2.*PI+ALFA CALFA=COS(ALFA) SALFA=SIN(ALFA) IF(ALFA. GE.Q. .AND. ALFA. LT. ALFAI)THEN RO=BP/CALFA ELSE IF(ALFA. GE. ALFA1 .AND. ALFA.LT. * ALFA2)THEN RO=AP/SALFA ELSE IF(ALFA. GE. ALFA2 .AND. ALFA.LT. * ALFA3)THEN RO=(BP-BB)/CALFA ELSE IF(ALFR. GE. ALFA3 .AND. ALFA.LT. * ALFA4)THEN RO=(AP-AA)/SALFA ELSE RO=BP/CALFA END I F
ZRAN=RAN(II)
ALFAI=Pl/2. ALFR2=PI-ATAN(AP/BB) ALFA3=PI+ATAN((I~-AP)/BB) PI.FA4=3.*Pl/2. ELSE IF(BP, EQ.B~)THEN ALFMIN=-PI/2. ALFMAX=PI/2. ALFAI=ATAN(AP/BB) ALFA2=PI/2. ALFA3-3.*PI/2. ALFA4=2.*PI-ATAN((AA-AP)/BB) ELSE ALFMIN=O. ALFMAX=2. * P I ALFAI=ATAN(AP/Bp) ALFR2=PI-ATAN(AP/
ELSE
ELSE
RLFR~PII2. RLFR3-PII2. RLFA4~PII2.
C C C 'C
DOUBLE FACE
IF(AS. EQ.O.) THEN CP=C-ZP ELSE CD=ZP END IF ELSE P~=Ym IF~A3. E Q . ~ . ) THEN CP=C-ZP ELSE CP=ZD END IF END IF ELSE I F ( A ~ . E Q . I . ) THEN AA=A BB=B CC=C BP=YP IF(A1.EQ.O.) THEN AP=A-XP IF(A3. EO.O.) THEN CP=C-ZP ELSE CP=ZP END I F ELSE AD=XP IF(A3. EQ.O.) THEN CP=C-ZD ELSE Cp=zP END I F END I F ELSE AA=A
THEN AA=B BB=A CC=C BP=XP IF(A2. EQ.O.) T~EN
200 CONTINUE lC(A!.EO.I.)
I!
WTETA=(I.-COS(TETMAX))/2. W=WTETA*WALFA SUM=SUM+W SUMSQ=SUMSQ+W*W CONTINUE GO TO 4 0 0
ZRAN=RAN(II) TETA=ACOS(I.-ZRAN*(!.-COS(TETMAX>))
TETMAX=ATAN(RO/CP)
C C C
C
C E C C C
C
3 0 0 CONTINUE AA-A BB=B CC=C
T R I P L E FACE
WTETA=(EOS(TETMIN)-COS(TETNAX))/2. W=WTETA*WALFA SUM=SUM+W SUMSQ=SUMSQ+W*W 21 CONTINUE GO TO 4 0 0
ZRAN=RAN(II) TETA=ACOS(COS(TETMIN)-ZRAN*(COS(TETHIN) COS(TETMAX)))
BB=C CC=B BP=C-ZP !F(A!,EQ.O.) THEN AP=A-XO IF(A2. EQ.O.) THEN CO=B-yP ELSE CD=YP END I F ELSE AP=XP ~ A 2 . E~.O.) THEN CP=B-YD ELSE CP=YP END IF END IF END IF ALFAI=ATAN(BP/(AP-AA)) ALFA2=ATAN((BB-BP)/(AP-AA)) RLFMAX=ALFAI+ALFA2 WALFA=ALFMAX/2./Pl DO 21 I = I , N ZRAN=RAN(II) ALFA=IRAN*ALFMAX TETMIN=ATAN((AP-AA)/COS(ALFAI-ALFA)/CP) ALFAII=ALFAI-ATAN(BP/AO) ALFA22=ALFAI+ATAN((BB-BP)/AP) IF(ALFA. LT. ALFA11)THEN RO=BP/SIN(ALFA1-ALFA) ELSE IF(ALFA. SE. ALFAll .AND. ALFA.LT° * ALFA22)THEN RO=AP/COS(ALFA1-ALFA) ELSE RO=(BB-BP)/SIN(ALFA-ALFA1) END IF TETMAX-ATAN(RO/(CP-CC))
~ ( B p - ~ ) J~p)
THEN
(CP-CC))
SERR=SQRT(S/N) T Y P E *~DLAN~SERR END
SLIMS@=SUMSQ+W*W C 31 CONTINUE 4 0 0 CONTINUE 5 0 0 CONTI N U E SDLAN-BU~/N S= (SUI~Q-N*SOLAN*SOLAN) / ( N - I )
ELSE TETMIN=ATAN((BP-BB)/COS(ALFA+ALFA1)/CP) END I F WTETA=(COS(TETMIN)-COS(TETMAX))/2. W=WTETA*WALFA SUM=SUM+W
TETNIN=ATRN((AP-AA)ISIN(ALFA+ALFAI)/CP)
*END IF I F ( A L F A . LE. ALFA22)
TETMAX=ATAN(AP/COSCALFMAX-ALFA+ALFA2)/
WALFA=ALFMAXI2./Pl, DO 31 I = I , N : ZRAN=RAN(II) " ALFA-ZRAN*ALFMAX IF(ALFA. LE. ALFA11) THEN TETMAX-ATAN(BP/COSIALFA+ALFA1)/IOP-CCII ELSE
ALFA22-ATAN(~AP-AA')/~BP-BB))-ALFAI
ALFMAX=PII2.@ ~-~FAI-ALFA2 ALFA11=ATAN(AP/BP)~ALFA1
~LFP~=AT~
THE~ CP=C-ZP IF(A2. EO.O.~ THEN BP=B-yP IF(~I.EQ.O.) AP=~q-XP IF(AI.EQ. 2 . ) AP=XP ELSE BP=Yp IF(AI.EQ.@o) AP=A-XP IF(A1.EQ. 2 . ) AP=XP END IF ELSE CP-ZP IF(A2. EQ.O.) THEN BP=B-YP IF(A1.EQ.O.) AP-A-XP IF(AI.EQ. 2 . ) AP=XP ELSE Bp=yp IF(A1.EQ.O.) AP=R-XP IF(A1.EQ. 2 . ) AP=XP END IF END IF ALFAI=ATAN((AA-AA)/BP)
IF(A3.EO.O.)
~.
~
~"
~.
~
448
L. Wielopolski / Monte Carlo calculation of auerage solid angle
It is a p p a r e n t from the table that all the o d d values of N V can be a d d r e s s e d to P2 calculations, except p o i n t 13 which represents a p o i n t inside the p a r a l l e l e p i p e d and is not being c o n s i d e r e d here. It also can be o b s e r v e d that in o r d e r to a d d r e s s a p o i n t for P1 calculations it suffices to assure that any two of the coefficients equal to one, a n d b y default all the o t h e r cases can be a d d r e s s e d to P3, calculations. This p r o c e d u r e is i m p l e m e n t e d in the p r o g r a m in table 2. T h e p r o g r a m in table 3 is written in F O R T R A N V language a n d it i m p l e m e n t s all the features discussed in the article. The p r o g r a m is n o t written as a general c o d e thus, different cases require small changes in the p r o g r a m itself. The p r o g r a m in table 3 represent the case for s a m p l i n g u n i f o r m l y a 1 c m x 1 cm square source. F o r a p o i n t source the four statements following P O I N T S E L E C T I O N F R O M 1 x 1 C M S Q U A R E have to be removed. F o r a specific case of a p o i n t source the calculations can be p e r f o r m e d faster b y r e m o v i n g D O - L O O P 500 a n d i m p l e m e n t i n g D O - L O O P s 11, 21, a n d 31 which are neutralized at present b y letter C in the first column.
References [1] [2] [3] [4] [5]
H. Gotoh and H. Yagi, Nucl. Instr. and Meth. 96 (1971) 485. G. Rowlands, Int. J. Appl. Radiat. Isot. 10 (1061) 86. C.R. Gillespie, Rev. Sci. Instr. 41 (1970) 42. J. Cook, Nucl. Instr. and Meth. 178 (1980) 561. L. Wielopolski, Nucl. Instr. and Meth. 143 (1977) 577.