KEK Internal 2021-***** 2021
R/D

Lecture Notes ofCalculation of air absorbed dosewith ionization chamber(cg Version)

H. Hirayama and Y. Namito

High Energy Accelerator Reserach Organization2©High Energy Accelerator Research Organization (KEK), 2021
KEK Reports are available fromHigh Energy Accelerator Research Organization (KEK)
1-1 Oho, Tsukuba-shi
Ibaraki-ken, 305-0801
JAPAN

Phone:
Fax:
E-mail:Internet:+81-29-864-5137+81-29-864-4604
irdpub@mail.kek.jp
https://www.kek.jp/en/Lecture Notes ofCalculation of air absorbed dosewith ionization chamber(cg Version)

H. Hirayama and Y. Namito

High Energy Accelerator Reserach Organization

Contents
Japanese Parts

1

1 Combinatrial Geometry (CG)
1.1 Body の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 リージョンの定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 リージョン定義の例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2
2
2
3

2 サンプルプログラム ucionch cgv.f の概要
2.1 CG 入力データ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

3 ユーザーコードの内容
3.1 メインプログラム: Step 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 メインプログラム: Step 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 メインプログラム: Step 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 メインプログラム: Step 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 メインプログラム: Step 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 メインプログラム: Step 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.7 メインプログラム: Step 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8 メインプログラム: Step 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.8.1 統計誤差 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.9 メインプログラム: Step 9 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.10 Subroutine ausgab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.11 Subroutine howfar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 function encoea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6
6
7
8
9
10
10
10
11
12
13
15
16
17

4 計算結果

17

5 実習課題5.1 実習課題１：線源光子のエネルギーを変えた計算 . . . . . . . . . . . . . . . . . .
5.2 実習課題 2：0.662MeV 線源について、一方向 (Z-方向) のみに放出している線源光子を、等方線源に変更した計算 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18
18

6 実習課題の解答例
6.1 実習課題１ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 実習課題 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19
19
20

Appendix: Full listings of ucionch cgv.f

24

i

18

EGS5 サンプルプログラム (ucionch cgv.f )電離箱による空気吸収線量の計算(Japanese Parts)1

1

Combinatrial Geometry (CG)

1.1

Body の定義

EGS 用 CG [1] では、以下のような立体（Body）を使用する事ができる。
1. 直方体 (RPP)x-, y- と z-方向の最小値及び最大値で定義する。各面はいずれかの軸と平行である。
2. 球 (SPH)
球の中心を示すベクトル V と半径で定義する。
3. 円筒 (RCC)
円筒の底面の中心を示すベクトル V と、中心からの高さベクトル H 及び円筒の半径で定義する。
4. 円錐台 (TRC)
円錐の底面の中心を示すベクトル V 、底面中心からの上面中心への高さベクトル H、 及び底面と上面のそれぞれの半径 R1 及び R2 で定義する。
5. トーラス (TOR)
いずれかの軸に平行なトーラスの中心を示すベクトル V、 トーラス中心から、チューブの中心までの距離 R1 、チューブの半径 R2 及びトーラスの方向を示す番号、 (n: x/y/z =1/2/3) で定義する。更に、トーラスの始まりの角度 θ1 と終わりの角度 θ2 を指定する。トーラス全体を使用する場合には、θ1=0, 及び θ2=2π とする。

Table 1: 各形状の立体とその記述のためのデータ
形状 通番
各形状の立体を定義するデータ
RPP #Xmin Xmax Ymin Ymax Zmin Zmax
SPH #
Vx
Vy
Vz
R
RCC #
Vx
Vy
Vz
Hx
Hy
Hz
R
TRC #
Vx
Vy
Vz
Hx
Hy
Hz
R1
R2
TOR #
Vx
Vy
Vz
R1
R2
θ1
θ2
n

1.2

リージョンの定義

各リージョンは、body の組み合わせにより定義する。組み合わせには、特別な記号、+, − 及びOR が使われる。
+ 記号の後に body 番号が書かれた場合には、body の内側の領域がリージョンとなる。一方、− 記号の後に body 番号が書かれた場合には、body の外側の領域がリージョンとなる。body 番
号の後に+又は − 記号と body 番号が続く場合には、間に AND 記号があるのと同じである。従って、+1 +2 は、body 1 の内側でなおかつ body 2 の内側を意味するので、body 1 と body 2 の重なった領域となる。一方、+1 −2 は、body 1 の内側でなおかつ body 2 の外側を意味するので、body 1 の領域中で body 2 と重なっていない領域を意味することになる。Body 番号が OR 記号の後に書かれた場合は、OR 記号は結合記号として使用される。リージョンが、OR 記号で結合したサブリージョンの組み合わせで定義される場合もある。２つ以上の OR 記号が使われる場合、OR の機能は、OR 記号の間及び OR 記号からリージョン定義の行の最後までに含まれる全てのbody 番号に、+ や − 記号に関係なく適用される。

2

1.3

リージョン定義の例

(a)

(b)

(c)

2

I

3

(d)
(e)
J

2

K

3

Figure 1: Combinatorial Geometry の例
第 1 図に示すような、球 (body 2) に円筒 (body 3) が挿入している様な体系を考える。もし、球と円筒の物質が同じであれば、リージョン I (図 1c) の様に一つのリージョンとする事ができる。
リージョン I は、
I = +2OR + 3と記述する。これは、リージョン I が、body 2 か body 3 のどちらかに属する領域である事を意味している。
球と円筒が異なった物質の場合、円筒部を除外した球には、円筒部のリージョン番号 (K) と異なったリージョン番号を付ける (例えば J)。
リージョン J (図 1d) は、
J = +2 − 3と記述する。これは、body 2 に属するが、body 3 に属さない領域を意味する。
リージョン K (図 2e) は、単に
K = +3と記述する。これは、body 3 の属する領域を意味する。

3

2 つ以上の body を組み合わせる場合には、+ 、− や OR 記号を含む長い記述となる。しかしながら、形状中の全ての点は、どれか一つのリージョンとして定義される様にしなければなら
ない。

B
2

A

3
1

Figure 2: Use of OR operator.
OR 記号を使ったもっと複雑な例として、第２図の斜線部の領域 A と斜線を引いていない領
域 B を考える。これらのリージョンは、２つの直方体 (body 1 と 3) と、一つの円筒 (body 2) で記述される。それぞれのリージョンは、
A = +1 + 2
そして

B = +3 − 1OR + 3 − 2と記述する。OR 記号は、次に OR 記号が現れるまで、それに続く全ての body 番号に適用される事に注意する必要がある。

4

サンプルプログラム ucionch cgv.f の概要

2

ucionch cgv.f は、CG を使って形状を記述するユーザコードである。CG 入力データは、ユニット
4 のデータファイルに記述する。

2.1

CG 入力データ

図 3 に示すように円筒の組み合わせで体系を定義している。

R

10.0
20 degree 1 atm air (region3)
5.2
5.0

Polyethylene case (region2)0.662 MeV
photon
0.0
-12.0
-0.2

20 degree 10 atm air(region1)
10.2
10.0
20.0

Z

Outer vacuum (region4)

Figure 3: ucionch cgv.f のジオメトリー
この形状の入力データは、以下のように記述する。

RCC

1

RCC

2

RCC

3

RCC

4

0.00
10.0
0.00
10.4
0.00
30.0
0.00
34.0

0.0
5.00
0.0
5.20
0.0
10.0
0.0
12.0

0.0

0.00

0.0

-0.2

0.00

0.0

-10.0

0.00

0.0

-12.0

0.00

0.0

END
5

Z1
Z2
Z3
Z4
END
1

+1
+2
+3
+4
2

1

-1
-2
-3
0

最後の行は、対応するリージョンの物質番号である。

1. 線源条件
 入射粒子は、エネルギー 0.662MeV の光子 入射粒子の角度は、Z-軸に沿って (0.0, 0.0, -5.0) から垂直入射

2. 得られる情報
 CGview 用の飛跡データ (egs5job.pic) 計算結果 (egs5job.out)– 使用する物質に関するデータ– 各リージョンに関する情報– 電離箱内の 1 気圧空気の吸収エネルギーから求めた空気吸収線量* LARCH 機能を使って、電離箱内の空気にエネルギーを付与する電子の生成場
所（電離箱内の空気中、前面のケース、後面のケース及び側面のケース）毎の吸収エネルギーと合計のエネルギーから求めた空気吸収線量* 空気のエネルギー質量エネルギー吸収係数と光子の飛程長から求めた空気吸収
線量* 光子の飛程長から求めた電離箱内での光子の平均スペクトル3

ユーザーコードの内容

3.1

メインプログラム: Step 1

egs5 は、Fortran で書かれているので、egs5 やジオメトリーや、ユーザーコードで使われている変数の配列の大きさは、別のファイルに parameter 文で指定し、include 機能によりユーザーコー
ドに取り入れている。common についても、同じく include 機能を用いている。
egs5 に直接関係する include 関係のファイルは、include/ディレクトリ（egs に関係するもの）、
pegscommons/（pegs に関係するもの）および auxcommons/（egs5 の著者から提供しているジオ
メトリー関係のサブルーティン等ユーザーコードにのみ関係するもの）とリンクすることにより、使用できるようにしている。1
この点が、Mortran のマクロ機能により、ユーザーコードで再設定できた EGS4 の場合と最も異なることである。従って、egs5 に直接関係する配列の大きさを変更場合は、include/egs5 h.f
内の、その他の場合は、auxcommons/aux h.f の当該 parameter 文の値を変更することになる。
最初の設定は、egs に直接関連する include 文である。
implicit none!
!
!

-----------EGS5 COMMONs-----------include ’include/egs5_h.f’! Main EGS "header" fileinclude ’include/egs5_bounds.f’include ’include/egs5_brempr.f’include ’include/egs5_edge.f’include ’include/egs5_media.f’
1

これらの設定は、egs5run スクリプト又は egs5run.bat で設定される。

6

include ’include/egs5_misc.f’include ’include/egs5_stack.f’
include ’include/egs5_thresh.f’include ’include/egs5_uphiot.f’include ’include/egs5_useful.f’
include ’include/egs5_usersc.f’include ’include/egs5_userxt.f’include ’include/randomm.f’include ’include/egs5 h.f’ は必ず必要であるが、それ以外の common に関連する include
文は、メインプログラムで使用する可能性があるものだけで良い。2次の設定は、ジオメトリー関係のサブルーティン等ユーザーコードに関連する include 文で
ある。
!
!
!

---------------------Auxiliary-code COMMONs---------------------include ’auxcommons/aux_h.f’! Auxiliary-code "header" fileinclude ’auxcommons/edata.f’
include ’auxcommons/etaly1.f’include ’auxcommons/instuf.f’
include ’auxcommons/lines.f’include ’auxcommons/watch.f’!
!
!

-----------------cg related COMMONs-----------------include ’auxcommons/geom_common.f’ ! geom-common file
integer irinn
最後の include 文が、CG に関連したもので、CG を使用する場合には常にこの表現とする。
ユーザーコード内で使用する common を次に定義する。

common/totals/* dose(5),dosep,deltae,spec(50),maxpict
real*8 dose,dosep,deltae,spec
integer maxpict! Variables to score

メインプログラムの先頭で、implicit none 宣言をしているので、メインプログラムで使用している全ての変数の型式宣言をする必要がある。
実行文の先頭で、使用するユニットを open する。egs5 では、pegs5 をプログラムの一部として含む構造を標準としている。pegs の実行に伴い、ユニット 7-26 は、close されることから、メインプログラムで open していても、pegs 実行後に、再度 open することが必要となる。そのため、ユニット 7-26 の使用を避ける方が良い。

!
!
!

---------Open files---------open(6,FILE=’egs5job.out’,STATUS=’unknown’)open(4,FILE=’egs5job.inp’,STATUS=’old’)open(39,FILE=’egs5job.pic’,STATUS=’unknown’)
ユニット 39 は、飛跡情報の出力ファイルである。
その後、各カウンターをリセットするサブルーティン counters out(0) を call する。

3.2

メインプログラム: Step 2

ユーザーコードで使用する物質数を nmed で定義する。
物質データ及び各物質の Characteristic Dimension を設定した後で、pegs5 を call する。medarrのデータは必ず 24 文字分を指定する必要がある。物質名が 24 文字未満の場合には空白を補って、
2

EGS4 の COMIN マクロに対応する扱いである。

7

合計 24 文字とする。Chracteristic Dimension は、当該物質で構成されるリージョンの最も小さい
サイズ (1 cm × 1 cm × 1 cm の立方体であれば 1cm) に設定する。

nmed=2if(nmed.gt.MXMED) thenwrite(6,’(A,I4,A,I4,A/A)’)*’ nmed (’,nmed,’) larger than MXMED (’,MXMED,’)’,*’ MXMED in iclude/egs5_h.f must be increased.’stop
end if!
!

==============
call block_set==============! Initialize some general variables

medarr(1)=’20 degree 1 atm air
medarr(2)=’POLYETHYLENE’’

do j=1,nmed
do i=1,24
media(i,j)=medarr(j)(i:i)
end doend dochard(1) = 5.0d0
chard(2) = 0.2d0!

automatic step-size control

write(6,fmt="(’chard =’,5e12.5)") (chard(j),j=1,nmed)
!
!
!
100
!

----------------------------------Run KEK PEGS5 before calling HATCH----------------------------------write(6,100)
FORMAT(’PEGS5-call comes next’/)!

==========
call pegs5==========

3.3

メインプログラム: Step 3

飛跡データファイルのフォーマットを指定する npreci を設定する。このユーザーコードでは、フリーフォーマットの 3 を指定する。計算結果の出力ファイルに、CG データの開始を示す CG dataを書き込み、その後 CG の入力データを読み込み、CG データを指定したファイル (この場合は、
6) に出力する処理を行うサブルーティン geomgt を call する。その後、CG データの終了を意味
する End of CG data を出力する。次に、飛跡データファイルに必要な情報を出力する。出力ユニットである ifto は、39 に設定している。PICT のデータモードを示す文字列 (CSTA-FREE 又は、CSTA) を出力し、再度 subroutine geomgt により CG データを飛跡データファイルに出力する。
最後に CG データの終了を意味する CEND を出力する。これらの処理後、cg データから、リージョン総数である nreg を引き出す。
CG を使用する場合には、この部分が必ず必要であり、変更する必要はない。

write(6,*) ’Read cg-related data’!----------------------------------------------!
initilize cg related parameter!----------------------------------------------npreci=3! PICT data mode for CGView in free format
ifti = 4
ifto = 39! Input unit number for cg-data! Output unit number for PICT

write(6,fmt="(’ CG data’)")call geomgt(ifti,6) ! Read in CG data
write(6,fmt="(’ End of CG data’,/)")

8

if(npreci.eq.3) write(ifto,fmt="(’CSTA-FREE-TIME’)")if(npreci.eq.2) write(ifto,fmt="(’CSTA-TIME’)")rewind ifti
call geomgt(ifti,ifto)! Dummy call to write geom info for ifto
write(ifto,’(A)’) ’CEND’!-------------------------------!
Get nreg from cg input data!-------------------------------nreg=izonin
物質番号の設定を、CG データの最後で定義した情報から読み込んだ後、各リージョンの物質
番号、egs カットオフエネルギー、オプションの設定（このユーザーコードでは、光電子の角度分
布、特性 X 線の発生）を行う。
anlux 乱数のシード inseed の値を設定し、初期化する。

!

Read material for each refion from egs5job.data
read(4,*) (med(i),i=1,nreg)!

Set option except vacuum region
do i=1,nreg
if(med(i).ne.0) theniphter(i) = 1! Switches for PE-angle sampling
iedgfl(i) = 1! K & L-edge fluorescenceiauger(i) = 0! K & L-Auger
iraylr(i) = 0! Rayleigh scattering
lpolar(i) = 0! Linearly-polarized photon scattering
incohr(i) = 0! S/Z rejection
iprofr(i) = 0! Doppler broadening
impacr(i) = 0! Electron impact ionization
end ifend do!
!
!
!

-------------------------------------------------------Random number seeds. Must be defined before call hatch
or defaults will be used. inseed (1- 2^31)-------------------------------------------------------luxlev = 1
inseed=1
write(6,120) inseed
FORMAT(/,’ inseed=’,I12,5X,*’ (seed for generating unique sequences of Ranlux)’)

120
!
!

=============
call rluxinit=============

3.4

メインプログラム: Step 4! Initialize the Ranlux random-number generator

入射粒子のパラメータを設定する。この例では、単一エネルギーの光子 (1.253MeV) が、垂直に検出器の中心に入射するとしている。なお CG を用いる場合には、irin=0 と指定することにより
irin を自動的に設定することができる。このユーザーコードでは線源位置は固定なので、irin=0の場合は、直後に irin を CG により求めるルーチンが書かれている。

! Define initial variables for incident particle normally incident! on the slab
iqin=0! Incident particle charge - photons
ekein=0.662! Incident particle kinetic energy
xin=0.0! Source position
yin=0.0
zin=-5.0
uin=0.0! Moving along z axis
vin=0.0
win=1.0
irin=0! Starting region (0: Automatic search in CG)
wtin=1.0! Weight = 1 since no variance reduction used9!

pdf data for many source
deltae=0.05! Energy bin of response!----------------------------------------!
Get source region from cg input data!----------------------------------------!
if(irin.le.0.or.irin.gt.nreg) thencall srzone(xin,yin,zin,iqin+2,0,irin)
if(irin.le.0.or.irin.ge.nreg) thenwrite(6,fmt="(’ Stopped in MAIN. irin = ’,i5)")irin
stop
end ifcall rstnxt(iqin+2,0,irin)
end if3.5

メインプログラム: Step 5

pegs で作成した物質データの最高エネルギー（全エネルギー）の最小値 emaxe を hatch で求める
ために、emaxe を 0 に設定後に、subroutine hatch を呼ぶ。emaxe は、入射粒子のエネルギーが、物質データの最高エネルギーを超えていないことをチェックするために用いる。hatch で読み込まれた物質データや、リージョンに設定した情報を確認のために出力するようにしている。
emaxe = 0.D0 ! dummy value to extract min(UE,UP+RM).
write(6,’(/A)’) ’ Call hatch to get cross-section data’!
!
!

-----------------------------Open files (before HATCH call)-----------------------------open(UNIT=KMPI,FILE=’pgs5job.pegs5dat’,STATUS=’old’)open(UNIT=KMPO,FILE=’egs5job.dummy’,STATUS=’unknown’)!
!

==========
call hatch==========

3.6

メインプログラム: Step 6

普通のユーザーコードでは、このステップで形状に関する情報（平板、円筒、球等）を記述する
が、本ユーザーコードでは cg で形状を指定しているので、このステップで記述する事項はない。

3.7

メインプログラム: Step 7

ausgab に必要な設定を行う。エネルギービンの幅を、入射エネルギーとビン数 (50) から計算する。Ncases は、ヒストリー数で、maxpict は、飛跡情報を記録するヒストリー数である。
!

Energy bin widthdeltae=ekein / 50!

Zero the variables
do i=1,5dose(i)=0.D0
doses(i)=0.D0
dose2s(i)=0.D0
end dodosep=0.D0
doseps=0.D0
dosep2s=0.D0
do j=1,50
spec(j)=0.D0
end do!

Set histories
ncases=10000000
10!

Set maximum number for pictmaxpict=50

3.8

メインプログラム: Step 8

飛跡表示ファイルにバッチ番号を出力した後で、ncases ヒストリーだけ subroutine shower を
call する。入射粒子の運動エネルギーに対応する電子の全エネルギーが物質データで設定した最高エネルギー emaxe を超えていないか調べる。各ヒストリー毎に、NaI 領域での吸収エネルギーの有無を調べ、吸収エネルギーがある場合には、全検出効率の数に wtin を加え、そのエネルギーが入射粒子の 99.9%以上の時は、ピーク検出効率の数に wtin を加える。また、吸収エネルギーの値により、波高分布の対応するチャンネルに wtin を加える。統計誤差評価のために、粒子スペク
トルを含め上記の変数のヒストリー毎の結果の自乗値の和を求めておく。

!
!
!

Write batch number
write(39,fmt="(’01’)")========================
if(iwatch.gt.0) call swatch(-99,iwatch)========================! ------------------------! Start of shower call-loop! -------------------------

do i=1,ncases
latchi=0!
!
!

---------------------Select incident energy---------------------wtin = 1.0
wtsum = wtsum + wtin
etot = ekein + iabs(iqin)*RM
if(iqin.eq.1) thenavailke = ekein + 2.0*RM
else
availke = ekein
end if! Keep running sum of weights! Incident total energy (MeV)! Available K.E. (MeV) in system! for positron! Available K.E. (MeV) in system! for photon and electron

totke = totke + availke! Keep running sum of KE!
!
!

---------------------Select incident angle----------------------!
!
!

--------------------------------------------------Print first NWRITE or NLINES, whichever comes first--------------------------------------------------if (ncount .le. nwrite .and. ilines .le. nlines) thenilines = ilines + 1
write(6,280) etot,xin,yin,zin,uin,vin,win,iqin,irin,idin
FORMAT(7G15.7,3I5)
end if

280
!
!
!
1

!
!
!

----------------------------------------------------------Compare maximum energy of material data and incident energy----------------------------------------------------------if(etot+(1-iabs(iqin))*RM.gt.emaxe) thenwrite(6,fmt="(’ Stopped in MAIN.’,’ (Incident kinetic energy + RM) > min(UE,UP+RM).’)")
stop
end if---------------------------------------------------Verify the normalization of source direction cosines---------------------------------------------------if(abs(uin*uin+vin*vin+win*win-1.0).gt.1.e-6) thenwrite(6,fmt="(’ Following source direction cosines are not’,
11

1

!

’ normalized.’,3e12.5)")uin,vin,win
stop
end if!

==========================================================
call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin)==========================================================!
!
!

---------------------------Sum variable and its square.
---------------------------do ie=1,5
doses(ie)=doses(ie)+dose(ie)dose2s(ie)=dose2s(ie)+dose(ie)*dose(ie)dose(ie)=0.d0
end dodoseps=doseps+dosep
dosep2s=dosep2s+dosep*dosep
dosep=0.d0
do ie=1,50
specs(ie)=specs(ie)+spec(ie)spec2s(ie)=spec2s(ie)+spec(ie)*spec(ie)spec(ie)=0.D0
end doncount = ncount + 1!

! Count total number of actual cases========================
if(iwatch.gt.0) call swatch(-1,iwatch)========================!

! ----------------------! End of CALL SHOWER loop
! -----------------------

end do!

========================
if(iwatch.gt.0) call swatch(-88,iwatch)========================!

call plotxyz(99,0,0,0.D0,0.D0,0.D0,0.D0,0,0.D0,0.D0)
write(39,fmt="(’9’)")! Set end of batch for CG View

300

tt=etime(tarray)tt1=tarray(1)cputime=tt1-tt0
write(6,300) cputime
format(’ Elapsed Time (sec)=’,G15.5)

3.8.1

統計誤差

x をモンテカルロ計算で計算したい量（スコアーする量）とする。モンテカルロ計算の結果には、
その統計誤差が必要である。ucnaicgv.f では、次のような MCNP で使用している方法を採用している。
 ヒストリー数を N とする。
 xi を i 番目のヒストリーの結果とする。
 x の平均値を計算する:

N
1 ∑
xi
N i=1

(1)

N
N
1 ∑
1 ∑
(xi − x)2 ≃ x2 − x2 (x2 =x2 ).
N − 1 i=1
N i=1 i(2)x= xi の分散値を以下の式から求める。:

s2 =

12

 x の分散値は、

s2x =

1 2
1
s ≃ [x2 − x2 ]
N
N

(3)

1 2
(x − x2 )]1/2
N

(4)

となる。

 統計誤差として、

sx ≃ [を用いる。

先の計算した結果とその自乗の和は、上記の処理に用いる。

3.9

メインプログラム: Step 9

得られた結果を処理して打ち出す処理を行う。電離箱中での空気 9 の吸収エネルギーから求めた発生光子当たりの「空気の吸収線量」、空気の質量エネルギー吸収係数と電離箱中での空気の飛程
長から求めた発生光子当たりの「空気の吸収線量」及び電離箱中での空気の飛程長から求めた電
離箱中での光子の平均スペクトルについて、ヒストリー毎の結果の和と結果の自乗和から、平均
値とその統計誤差を求めて出力する。
光子の飛程長を使った計算では、各粒子の検出器中での飛程長 (ustep) に光子のエネルギーと空気の質量エネルギー吸収係数を掛けたの合計及び、光子のエネルギー区分毎の飛程長の合計を
電離箱の体積で割ることにより、電離箱中での吸収エネルギー及び平均の光子スペクトルを求めることができる。

330

tdet=10.0
rdet=5.0
tcase=0.2
rcase=0.2
write(6,330) tdet,rdet,tcase,rcase
FORMAT(/’ Detector length=’,G15.5,’ cm’/*
’ Detector radius=’,G15.5,’ cm’/*
’ Polyetylene case thickness=’,G10.2,’ cm’/*
’ Polyetylene case side thickness=’,G10.2,’ cm’/)340

write(6,340) ekein
FORMAT(’ Results for ’,G15.5,’MeV photon’/)350

write(6,350) (media(i,2),i=1,24)
FORMAT(’ Cover material is ’,24A1/)

360

write(6,360) (media(i,1),i=1,24)
FORMAT(’ Inside chamber is ’,24A1/)

365

write(6,365) zin
FORMAT(’ Cs-137 exists at

!
!
!
!
!
!
!

370
!
!

z=’,F8.2,’ cm, x=y=0.0’/)----------------------------------Calculate average and its deviation----------------------------------Absorbed energy (MeV) devided by weight of air inside chamber (dweight)
MeV/g
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
tdet=10.0
rdet=5.0
vol=pi*rdet*rdet*tdet! volume of air cm3
weight=vol*rhom(med(1))! weight of air inside chamber g! rhom(1) is a density of air (medium 1)
write(6,370) vol,weight
FORMAT(’ Volume of chamber : ’,G15.5,’ cm3’/*’ weight of chamber : ’,G15.5,’ g’)--------------Absorbed dose
13

!
380

--------------write(6,380)
FORMAT(/’ Result from absorbed energy in air inside chamber’)write(6,’(A)’) ’ Air absorbed dose from ’do ie=1,5
avdose = doses(ie)/ncountsigdose=dose2s(ie)/ncountsigdose=dsqrt((sigdose-avdose*avdose)/ncount)
avdose = avdose/weight*1.602E-10
sigdose = sigdose/weight*1.602E-10
if(ie.eq.1) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron produced in air =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.2) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from front wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.3) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from bottom wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.4) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from side wall=’,*
avdose,’+-’,sigdose,’ Gy/incident photon’
else
write(6,’(A,G15.5,A,G15.4,A)’)*’ Total=’,*avdose,’+-’,sigdose,’ Gy/incident photon’
end ifend do

400

!
!
!
!
!
!
!

410
!
!
!
420

write(6,400)
FORMAT(/’ Result from absorbed energy in air’/*’using track length of photon inside chamber and’,*’ mass energy absorption coefficient of air’)Track lenght (cm) gives photon fluence deviding by volume
Conversion from air absorbed dose in MeV cm2/g to that in Gy
Unit of mass energy absorption coefficient mu_en is cm2/g
Ausgab scores energy (MeV) times mu_en in unit of MeV cm2/g.
Main routine converts this into Gy (J/kg).
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
avdose = doseps/ncountdosep2s=dosep2s/ncountsigdose=dsqrt((dosep2s-avdose*avdose)/ncount)
avdose = avdose/vol*1.602E-10
sigdose = sigdose/vol*1.602E-10
write(6,410) avdose,sigdose
FORMAT(’ Air absorbed dose =’,G15.4,’+-’,G15.4,*’ Gy/incident photon’)---------------------------------------Average photon spectrum inside chamber.
---------------------------------------write(6,420)
FORMAT(/’ Average photon spectrum inside chamber’/)do ie=1,50
elow=deltae*(ie-1)
eup=deltae*ie
avspe = specs(ie)/ncountspec2s(ie)=spec2s(ie)/ncount14

sigspe=dsqrt((spec2s(ie)-avspe*avspe)/ncount)avspe=avspe/vol/deltae
sigspe= sigspe/vol/deltae
430

3.10

write(6,430) eup,avspe,sigspe
FORMAT(’ Upper energy’, G15.5,’ MeV--’,G12.5,’+-’,G12.5,*’ photons/cm2/MeV’)end doSubroutine ausgab

AUSGAB は、ユーザが求める情報をスコアするサブルーチンである。最初に、メインプログラムと同様に、include 文及びローカル変数の型式宣言を行う。
iarg < 5 の場合には、それぞれのリージョンでの吸収エネルギーを計算する。
電子が周囲の壁から空気中に入ってきた（現在のリージョン (irl) が前のリージョンと異なり、前のリージョン番号が周囲のケース 2）場合には、粒子の z(np) により、前面 (1)、後面 (2) あるいは側面 (3) かを調べ、当該番号を latch(np) とする。電離箱内の空気中で発生した電子の場合は、
latch(np) は 0 のままである。
if(iql.eq.-1.and.irl.ne.irold.and.irold.eq.2) thenif(z(np).le.0.0) thenlatch(np)=1
elseif(z(np).ge.10.0) thenlatch(np)=2
else
latch(np)=3
end ifend if
リージョン番号が 1 (電離箱内の空気）の時は、ステップ中での吸収エネルギーを、latch(np)を使って、電子の生成場所毎に加算する。dose(5) は、合計なので、latch(np) に関係なく常に吸収エネルギーを加算する。

if (irl .eq. 1) then ! inside ion chamberdose(5) = dose(5) + edepwt
if(latch(np).ge.0.and.latch(np).le.3) thendose(latch(np)+1) = dose(latch(np)+1) + edepwt
end if
リージョン 1 内で光子が移動する場合には、当該光子のエネルギーに対応する空気の質量エネ
ルギー吸収係数 (dcon) を求め、エネルギーと光子の飛程長 (ustep) を使って、空気の吸収エネル
ギーを求め、dosep に加算する。また、電離箱内の空気中での平均光子スペクトルを得るために、光子のエネルギー区分毎の飛程長を加算する。

if (iql .eq. 0) then! photon
dcon=encoea(e(np))
dosep=dosep+wt(np)*e(np)*dcon*ustep
ie = e(np)/deltae +1
if(ie .gt. 50) ie = 50
spec(ie) = spec(ie) + wt(np)*e(np)*dcon*ustep
end if
ヒストリー数が maxpict 以下の時は、飛積情報を出力するために plotxyz を call する。

!
!
!

-----------------------Set some local variables-----------------------irl = ir(np)iql = iq(np)edepwt = edep*wt(np)if(iql.eq.-1.and.irl.ne.irold.and.irold.eq.2) then15

if(z(np).le.0.0) thenlatch(np)=1
elseif(z(np).ge.10.0) thenlatch(np)=2
else
latch(np)=3
end ifend if!
!
!

----------------------------------------------------------Keep track of energy deposition (for conservation purposes)----------------------------------------------------------if (iarg .lt. 5) thenesum(iql+2,irl,iarg+1) = esum(iql+2,irl,iarg+1) + edepwt
nsum(iql+2,irl,iarg+1) = nsum(iql+2,irl,iarg+1) + 1
end if!
!
!
!

----------------------------------------------------------------Print out particle transport information (if switch is turned on)----------------------------------------------------------------========================
if (iwatch .gt. 0) call swatch(iarg,iwatch)========================
if(iarg .ge. 5) return!

if (irl .eq. 1) then ! inside ion chamberdose(5) = dose(5) + edepwt
if(latch(np).ge.0.and.latch(np).le.3) thendose(latch(np)+1) = dose(latch(np)+1) + edepwt
end if!
!
!

-------------------------------------------------------Score particle information if it proceeds inside chamber-------------------------------------------------------if (iql .eq. 0) then! photon
dcon=encoea(e(np))
dosep=dosep+wt(np)*e(np)*dcon*ustep
ie = e(np)/deltae +1
if(ie .gt. 50) ie = 50
spec(ie) = spec(ie) + wt(np)*e(np)*dcon*ustep
end ifend if!
!
!

---------------------------------------------------------------Print out stack information (for limited number cases and lines)---------------------------------------------------------------if (ncount .le. nwrite .and. ilines .le. nlines) thenilines = ilines + 1
write(6,100) e(np),x(np),y(np),z(np),u(np),v(np),w(np),*iql,irl,iarg
100
FORMAT(7G15.7,3I5)
end if!
!
!

-----------------------------------Output particle information for plot-----------------------------------if (ncount.le.maxpict) thencall plotxyz(iarg,np,iq(np),x(np),y(np),z(np),e(np),ir(np),*wt(np),time(np))
end ifreturnend

3.11

Subroutine howfar

CG を利用するかぎりユーザーが howfar を変更する必要は一切ない。
以下、参考のため howfar の機能を述べる。howfar は、粒子の進行方向でのリージョン境界ま16

での距離を計算し、反応点までの距離との比較をし、境界までの距離の方が短い場合には粒子の移動距離を境界までの距離に置き換え、リージョンが変わるという処理を行う。
その他に、howfar では、ユーザが粒子の追跡を止める設定を行う (idisc=1) 。通常は、粒子が検討している領域の外に出て追跡を終了する場合にこの設定を行う。

3.12

function encoea

光子のエネルギーに対応する空気の質量エネルギー吸収係数を log-log 内挿で計算する function である。S. M. Seltzer と J. H. Hubbell のデータ [3] を用いている。このデータは、NIST のホームページ
http://www.physics.nist.gov/PhysRefData/Xcom/html/xcom1-t.html
で公開されているデータと同じものである。

計算結果

4

 ucionch cgv.f を egs5run で実行する。

– Linux 又は Cygwin の場合
ユーザーコード名として、ucionch cgv を, ユニット 4 及びユニット 25 のファイル名には，何も入力しないでリターンする。
”Does this user code read from the terminal?”に対して 1 を入力する。
– DOS の場合
egs5run ucionch cgv– ucionch cgv 等が、egs5run.bat を実行しているディレクトリーと別なディレクトリーにある場合は、ディレクトリー名からを記載する。DOS の場合、ディレクトリーの識
別子は、/ ではなく =Yであるので、間違わないように注意する。
 計算が終了したら、egs5job.out を調べる。電離箱中の空気の吸収エネルギーから求めた空気の吸収線量と空気の質量エネルギー吸収係数と電離箱中での光子の飛程長から求めた空気の吸収線量を下の表に示す。
計算結果
空気吸収線量 (Gy per source photon)
空気吸収エネルギー使用
電離箱内の空気中で生成した電子
前面のケースで生成した電子
後面のケースで生成した電子
側面のケースで生成した電子
合計
光子飛程長使用

0.290E-14 ± 0.397E-16
0.757E-14 ± 0.686E-16
0.899E-15 ± 0.271E-16
0.894E-14 ± 0.373E-16
0.203E-13 ± 0.120E-15
0.392E-13 ± 0.143E-17

空気の吸収エネルギーから求めた結果が小さいことがわかる。
なお、計算機、コンパイラーの違いにより、計算の状況がかわり、数値は統計誤差程度の変動を示すことがありうる。

17

5

実習課題

5.1

実習課題１：線源光子のエネルギーを変えた計算

次のように変更して、2 つの空気吸収線量の関係を調べよ。

1. 線源を、0.2MeV) に変える。
2. 線源を、1.0 MeV に変える。

5.2

実習課題 2：0.662MeV 線源について、一方向 (Z-方向) のみに放出している線源
光子を、等方線源に変更した計算

線源位置を次のように変更して 1 Bq 当たりの 2 つの空気吸収線量率の関係を調べよ。

1. -2.0 cm
2. -5.0 cm
3. -10.0 cm。

18

実習課題の解答例

6

比較のために、ucionch cgv.f の計算モードでの結果 (egs5job.out) を別な名称のファイル名 (例えば、chamber 662.out) で保存しておく。

6.1

実習課題１

1. ファイルのの変更
 cp ucionch cgv.f uciobch cgv1.f
これは、UNIX 又は Cygwin の場合である。DOS の場合は、copy ucionch cgv.f uciobch cgv1.f
又は、Windows 上でファイルのコピーを行う。以下の操作でも同様。
 cp ucionch cgv.data uciobch cgv1.data
 cp uciobch cgv.inp uciobch cgv1.inp
 uciobch cgv1.f を以下のように変更する。

– uciobch cgv1.f 中のekein=0.662
を
ekein=0.2
又は、
ekein=1.0に変更する。

! Incident particle kinetic energy! Incident particle kinetic energy! Incident particle kinetic energy

2. uciobch cgv1.f を egs5run で実行する。
 Linux 又は Cygwin の場合
ユーザーコード名として、uciobch cgv1 を, ユニット 4 及びユニット 25 のファイル名には，何も入力しないでリターンする。
”Does this user code read from the terminal?”に対して 1 を入力する。
 DOS の場合
egs5run uciobch cgv1
 uciobch cgv1 等が、egs5run.bat を実行しているディレクトリーと別なディレクトリーにある場合は、ディレクトリー名から記載する。DOS の場合、ディレクトリーの識別
子は、/ ではなく =Yであるので、間違わないように注意する。

3. 計算が終了したら、egs5job.out を調べ、0.662 MeV の結果と比較する。
空気吸収エネルギーの比較

0.662MeV
空気吸収エネルギー使用
電離箱内の空気中で生成した電子
前面のケースで生成した電子
後面のケースで生成した電子
側面のケースで生成した電子
合計
光子飛程長使用

空気吸収線量 (Gy per source photon)
0.200 MeV
1.00 MeV

0.290E-14 ± 0.397E-16
0.757E-14 ± 0.686E-16
0.899E-15 ± 0.271E-16
0.894E-14 ± 0.373E-16
0.203E-13 ± 0.120E-15
0.392E-13 ± 0.541E-17

190.775E-14 ± 0.801E-16
0.110E-14 ± 0.318E-16
0.125E-15 ± 0.834E-17
0.117E-14 ± 0.328E-16
0.101E-13 ± 0.397E-16
0.107E-13 ± 0.124E-17

0.194E-14 ± 0.303E-16
0.107E-13 ± 0.727E-16
0.114E-14 ± 0.291E-16
0.128E-13 ± 0.872E-16
0.266E-13 ± 0.127E-15
0.565E-13 ± 0.181E-17

6.2

実習課題 2

1. ファイルの変更
 cp uciobch cgv.f uciobch cgv2.f
 cp uciobch cgv.data uciobch cgv2.data
 cp uciobch cgv.inp uciobch cgv2.inp
 uciobch cgv2.f を以下の様に変更する。

– 崩壊当たりの γ 線放出率 erate を変数に追加する。
real*8! Local variables* availke,doses,dose2s,doseps,dosep2s,avdose,sigdose,avspe,sigspe,* rr0,wtin,wtsum, xi0,yi0,zi0,vol,weightを
real*8! Local variables* availke,doses,dose2s,doseps,dosep2s,avdose,sigdose,avspe,sigspe,* rr0,wtin,wtsum, xi0,yi0,zi0,vol,weight,erateに変更する。
– erate の値を指定する data 文の挿入
character*24 medarr(MXMED)
を
character*24 medarr(MXMED)
data erate/0.851/に変更する。
– 線源位置の変更-5.0 cm 以外の場合は、
zin=-5.0を
zin=-2.0
又は、
zin=-10.0に変更する。
– 4π での等方角度分布をサンプリングする文を追加する。
!
---------------------!
Select incident angle!
---------------------を!
---------------------!
Select incident angle!
---------------------275
call randomset(rnnow)
zi0=2.0*rnnow-1.0
call randomset(rnnow)
xi0=2.0*rnnow-1.0
call randomset(rnnow)
yi0=2.0*rnnow-1.0
rr0=dsqrt(xi0*xi0+yi0*yi0+zi0*zi0)
if(rr0.gt.1.0) go to 275
win = zi0/rr0
uin = xi0/rr0
vin = yi0/rr0に変更。
– 出力部の変更
20

do ie=1,5
avdose = doses(ie)/ncountsigdose=dose2s(ie)/ncountsigdose=dsqrt((sigdose-avdose*avdose)/ncount)
avdose = avdose/weight*1.602E-10
sigdose = sigdose/weight*1.602E-10
if(ie.eq.1) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron produced in air =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.2) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from front wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.3) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from bottom wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.4) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from side wall=’,*
avdose,’+-’,sigdose,’ Gy/incident photon’
else
write(6,’(A,G15.5,A,G15.4,A)’)*’ Total=’,*avdose,’+-’,sigdose,’ Gy/incident photon’
end ifend do
400

!
!
!
!
!
!
!

410

write(6,400)
FORMAT(/’ Result from absorbed energy in air’/*’using track length of photon inside chamber and’,*’ mass energy absorption coefficient of air’)Track lenght (cm) gives photon fluence deviding by volume
Conversion from air absorbed dose in MeV cm2/g to that in Gy
Unit of mass energy absorption coefficient mu_en is cm2/g
Ausgab scores energy (MeV) times mu_en in unit of MeV cm2/g.
Main routine converts this into Gy (J/kg).
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
avdose = doseps/ncountdosep2s=dosep2s/ncountsigdose=dsqrt((dosep2s-avdose*avdose)/ncount)
avdose = avdose/vol*1.602E-10
sigdose = sigdose/vol*1.602E-10
write(6,410) avdose,sigdose
FORMAT(’ Air absorbed dose =’,G15.4,’+-’,G15.4,*’ Gy/source photon’)を

do ie=1,5
avdose = doses(ie)/ncountsigdose=dose2s(ie)/ncountsigdose=dsqrt((sigdose-avdose*avdose)/ncount)
avdose = avdose*erate/weight*1.602E-10*3600.0
sigdose = sigdose*erate/weight*1.602E-10*3600.0
if(ie.eq.1) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron produced in air =’,*avdose,’+-’,sigdose,’ Gy/h per Bq’
elseif(ie.eq.2) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from front wall =’,*avdose,’+-’,sigdose,’ Gy/h per Bq’
elseif(ie.eq.3) thenwrite(6,’(A,G15.5,A,G15.4,A)’)
21

*
*

’ Electron from bottom wall =’,avdose,’+-’,sigdose,’ Gy/h per Bq’
elseif(ie.eq.4) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from side wall=’,*
avdose,’+-’,sigdose,’ Gy/h per Bq’
else
write(6,’(A,G15.5,A,G15.4,A)’)*’ Total=’,*avdose,’+-’,sigdose,’ Gy/h per Bq’
end ifend do

400

!
!
!
!
!
!
!

410
!
!
!
420

write(6,400)
FORMAT(/’ Result from absorbed energy in air’/*’using track length of photon inside chamber and’,*’ mass energy absorption coefficient of air’)Track lenght (cm) gives photon fluence deviding by volume
Conversion from air absorbed dose in MeV cm2/g to that in Gy
Unit of mass energy absorption coefficient mu_en is cm2/g
Ausgab scores energy (MeV) times mu_en in unit of MeV cm2/g.
Main routine converts this into Gy (J/kg).
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
avdose = doseps/ncountsigdose=dosep2s/ncountsigdose=dsqrt((sigdose-avdose*avdose)/ncount)
avdose = avdose*erate/vol*1.602E-10*3600.0
sigdose = sigdose*erate/vol*1.602E-10*3600.0
write(6,410) avdose,sigdose
FORMAT(’ Air absorbed dose rate =’,G15.4,’+-’,G15.4,*’ Gy/h per Bq’)---------------------------------------Average photon spectrum inside chamber.
---------------------------------------write(6,420)
FORMAT(/’ Average photon spectrum inside chamber’/)do ie=1,50
elow=deltae*(ie-1)
eup=deltae*ie

avspe = specs(ie)/ncountspec2s(ie)=spec2s(ie)/ncountsigspe=dsqrt((spec2s(ie)-avspe*avspe)/ncount)avspe=avspe*erate/vol/deltae
sigspe= sigspe*erate/vol/deltae
write(6,430) eup,avspe,sigspe
430
FORMAT(’ Upper energy’, G15.5,’ MeV--’,G12.5,’+-’,G12.5,*’ photons/cm2/MeV/s per Bq’)end doに変更する。
2. uciobch cgv2.f を egs5run で実行する。
3. 計算が終了したら、egs5job.out を調べ、各宣言位置で較する。
空気吸収エネルギーの比較

22空気吸収線量率 (Gy/ per Bq)- 5 cm-2 cm
空気吸収エネルギー使用
電離箱内の空気中で生成した電子
前面のケースで生成した電子
後面のケースで生成した電子
側面のケースで生成した電子
合計
光子飛程長使用

0.100E-11 ± 0.400E-13
0.612E-11 ± 0.108E-12
0.352E-12 ± 0.295E-13
0.101E-10 ± 0.151E-12
0.175E-10 ± 0.202E-120.178E-10 ± 0.111E-13

0.554E-12 ± 0.309E-13
0.258E-11 ± 0.700E-13
0.167E-12 ± 0.191E-13
0.487E-11 ± 0.106E-12
0.817E-11 ± 0.137E-12
0.831E-11 ± 0.832E-14-10 cm0.209E-12 ± 0.184E-13
0.923E-12 ± 0.402E-13
0.942E-13 ± 0.141E-13
0.228E-11 ± 0.740E-13
0.351E-11 ± 0.911E-13
0.352E-11 ± 0.586E-14

References[1] T. Torii and T. Sugita, “Development of PRESTA-CG Incorprating Combinatorial Geometry in EGS4/PRESTA”, JNC TN1410 2002-201, Japan Nuclear Cycle Development
Institute (2002).
[2] T. Sugita, T. Torii, A. Takamura, “Incorporating Combinatorial Geometry to the EGS5
Code and Its Speed-Up”, Twelfth EGS User’s Meeting in Japan, KEK Proc. 2005-10, 7-21,(KEK, Tsukuba, 9 - 11 Aug. 2005).
[3] S. M. Seltzer and J. H. Hubbell, “Tables and Graphs of photon mass attenuation coeﬃcientsand energy-absorption coeﬃcients for photon energies 1 keV to 20 MeV for elements Z=1
to 92 and some dosimetric materials”, 1995 Japanese Society of Radiological Technology.

23

Appendix: Full listings of ucionch cgv.f!***********************************************************************!***************************** KEK, High Energy Accelerator Research *!***************************** Organization*!* u c i o n c h _ c g v ******!*****************************
EGS5.0 USER CODE - 30 May 2021/1430 *!
correct score condition*!***********************************************************************!* This is a general User Code based on the cg geometry scheme.
*!***********************************************************************!
*! PROGRAMMERS: H. Hirayama*!
Applied Research Laboratory*!
KEK, High Energy Accelerator Research Organization*!
1-1, Oho, Tsukuba, Ibaraki, 305-0801
*
!
Japan
*
!
*
!
E-mail:hideo.hirayama@kek.jp*!
Telephone:+81-29-864-5451*!
Fax:+81-29-864-4051*!
*!***********************************************************************!***********************************************************************! The ucionch_cgv.f User Code requires a cg-input file only*
! (e.g., ucionch_cgv.data).
*! The following shows the geometry for ucnaicg.data.
*! Input data for CG geometry must be written at the top of data-input *! file together with material assignment to each region. Cg-data can *! be checked by CGview.
*! Use Ranlux random number generator.
*!***********************************************************************!
*!
------------------------*!
cg Geometry (ucionch_cgv)*!
------------------------*
!
*
!
*
!
R
*
!
^
*
!
|
*
!
+----+----+----+--------+------+--*
!
|
*
!
|
Outer vacuum region*!
++----+----+--------+------+-----+ R=10.0
*
!
|
|
Air
|
*
!
|
|
+--------------------++ R=5.2
*
!
|
|
|
Case
|
|
*
!
|
+
+
+-----------+++ R=5.0*!
|||| 1 atm
|
|
|
*
!
|
|
|
| Air
|
|
|
*
!
|
|
|
|
|
|
|
*
!
|
|
|
|
|
|
|
*
! 0.662 MeV |
|
|
|
|
|
|
*
!
==============>+----+-----------+---+-----+--------> Z*!
photons -10.0 -0.2 -0.0 0.0
10.0 10.2 15.0
*
!
*
!
*
!
*
!
*
!***********************************************************************!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12!----------------------------------------------------------------------!------------------------------- main code ----------------------------!----------------------------------------------------------------------!----------------------------------------------------------------------! Step 1: Initialization!----------------------------------------------------------------------implicit none!
!
!

-----------EGS5 COMMONs-----------include ’include/egs5_h.f’! Main EGS "header" fileinclude ’include/egs5_bounds.f’24include ’include/egs5_brempr.f’include ’include/egs5_edge.f’include ’include/egs5_media.f’
include ’include/egs5_misc.f’include ’include/egs5_stack.f’
include ’include/egs5_thresh.f’include ’include/egs5_uphiot.f’include ’include/egs5_useful.f’
include ’include/egs5_usersc.f’include ’include/egs5_userxt.f’include ’include/randomm.f’!
!
!

---------------------Auxiliary-code COMMONs---------------------include ’auxcommons/aux_h.f’! Auxiliary-code "header" fileinclude ’auxcommons/edata.f’
include ’auxcommons/etaly1.f’include ’auxcommons/instuf.f’
include ’auxcommons/lines.f’include ’auxcommons/nfac.f’include ’auxcommons/watch.f’!
!
!

-----------------cg related COMMONs-----------------include ’auxcommons/geom_common.f’ ! geom-common file
integer irinn
common/totals/* dose(5),dosep,deltae,spec(50),maxpict
real*8 dose,dosep,deltae,spec
integer maxpict! Variables to score!**** real*8
real*8 totke
real*8 rnnow,etot
real*8 esumt! Arguments

real*8! Local variables* availke,doseps,dosep2s,avdose,sigdose,avspe,sigspe,* rr0,wtin,wtsum, xi0,yi0,zi0,vol,weight
real*8* specs(50),spec2s(50),doses(5),dose2s(5)real* elow,eup,rdet,rcase,tcase,tdet! Local variables

real* tarray(2),tt,tt0,tt1,cputime,etime
integer* i,icases,idin,ie,ifti,ifto,ii,iiz,imed,ireg,isam,* izn,nlist,j,k,n,ner
character*24 medarr(MXMED)!
---------!
Open files!
---------!---------------------------------------------------------------!
Units 7-26 are used in pegs and closed. It is better not!
to use as output file. If they are used, they must be opened!
after getcg etc. Unit for pict must be 39.
!---------------------------------------------------------------open(6,FILE=’egs5job.out’,STATUS=’unknown’)open(4,FILE=’egs5job.inp’,STATUS=’old’)open(39,FILE=’egs5job.pic’,STATUS=’unknown’)!
!

====================
call counters_out(0)====================!----------------------------------------------------------------------! Step 2: pegs5-call!----------------------------------------------------------------------!

---------------------------------

25

!
!

Define media before calling PEGS5--------------------------------nmed=2
if(nmed.gt.MXMED) thenwrite(6,’(A,I4,A,I4,A/A)’)*’ nmed (’,nmed,’) larger than MXMED (’,MXMED,’)’,*’ MXMED in iclude/egs5_h.f must be increased.’stop
end if!
!

==============
call block_set==============! Initialize some general variables

medarr(1)=’20 degree 1 atm air
medarr(2)=’POLYETHYLENE’’

do j=1,nmed
do i=1,24
media(i,j)=medarr(j)(i:i)
end doend dochard(1) = 5.0d0
chard(2) = 0.2d0!

automatic step-size control

write(6,fmt="(’chard =’,5e12.5)") (chard(j),j=1,nmed)
!
!
!
100
!
!

----------------------------------Run KEK PEGS5 before calling HATCH----------------------------------write(6,100)
FORMAT(’PEGS5-call comes next’/)==========
call pegs5==========!----------------------------------------------------------------------! Step 3: Pre-hatch-call-initialization!----------------------------------------------------------------------write(6,*) ’Read cg-related data’!----------------------------------------------!
Initialize CG related parameters!----------------------------------------------npreci=3! PICT data mode for CGView in free format
ifti = 4
ifto = 39! Input unit number for cg-data! Output unit number for PICT

write(6,fmt="(’ CG data’)")call geomgt(ifti,6) ! Read in CG data
write(6,fmt="(’ End of CG data’,/)")
if(npreci.eq.3) write(ifto,fmt="(’CSTA-FREE-TIME’)")if(npreci.eq.2) write(ifto,fmt="(’CSTA-TIME’)")110

rewind ifti
call geomgt(ifti,ifto)! Dummy call to write geom info for ifto
write(ifto,110)
FORMAT(’CEND’)!-------------------------------!
Get nreg from cg input data!-------------------------------nreg=izonin!

Read material for each refion from egs5job.data
read(4,*) (med(i),i=1,nreg)!

Set option except vacuum region
do i=1,nreg
if(med(i).ne.0) theniphter(i) = 1! Switches for PE-angle sampling
iedgfl(i) = 1! K & L-edge fluorescenceiauger(i) = 0! K & L-Auger
iraylr(i) = 0! Rayleigh scattering
lpolar(i) = 0! Linearly-polarized photon scattering

26

incohr(i) = 0
iprofr(i) = 0
impacr(i) = 0
end ifend do
!
!
!
!

120
!
!

! S/Z rejection! Doppler broadening! Electron impact ionization-------------------------------------------------------Random number seeds. Must be defined before call hatch
or defaults will be used. inseed (1- 2^31)-------------------------------------------------------luxlev = 1
inseed=1
write(6,120) inseed
FORMAT(/,’ inseed=’,I12,5X,*’ (seed for generating unique sequences of Ranlux)’)=============
call rluxinit=============! Initialize the Ranlux random-number generator!----------------------------------------------------------------------! Step 4: Determination-of-incident-particle-parameters!----------------------------------------------------------------------! Define initial variables for incident particle normally incident! on the slab
iqin=0! Incident particle charge - photons
ekein=0.662! Incident particle kinetic energy
xin=0.0! Source position
yin=0.0
zin=-5.0
uin=0.0! Moving along z axis
vin=0.0
win=1.0
irin=0! Starting region (0: Automatic search in CG)
wtin=1.0! Weight = 1 since no variance reduction used!

pdf data for many source
deltae=0.05! Energy bin of response!----------------------------------------!
Get source region from cg input data!----------------------------------------!
if(irin.le.0.or.irin.gt.nreg) thencall srzone(xin,yin,zin,iqin+2,0,irin)
if(irin.le.0.or.irin.ge.nreg) thenwrite(6,fmt="(’ Stopped in MAIN. irin = ’,i5)")irin
stop
end ifcall rstnxt(iqin+2,0,irin)
end if!----------------------------------------------------------------------! Step 5:
hatch-call!----------------------------------------------------------------------emaxe = 0.D0 ! dummy value to extract min(UE,UP+RM).
130
!
!
!

140
!
!
!
!
!

write(6,130)
format(/’ Call hatch to get cross-section data’)-----------------------------Open files (before HATCH call)-----------------------------open(UNIT=KMPI,FILE=’pgs5job.pegs5dat’,STATUS=’old’)open(UNIT=KMPO,FILE=’egs5job.dummy’,STATUS=’unknown’)
write(6,140)
FORMAT(/,’ HATCH-call comes next’,/)==========
call hatch==========-----------------------------Close files (after HATCH call)-----------------------------close(UNIT=KMPI)
close(UNIT=KMPO)! ---------------------------------------------------------! Print various data associated with each media (not region)! ----------------------------------------------------------

27

150
160
170
180
190

write(6,150)
FORMAT(/,’ Quantities associated with each MEDIA:’)do j=1,nmed
write(6,160) (media(i,j),i=1,24)
FORMAT(/,1X,24A1)
write(6,170) rhom(j),rlcm(j)FORMAT(5X,’ rho=’,G15.7,’ g/cu.cm
rlc=’,G15.7,’ cm’)write(6,180) ae(j),ue(j)FORMAT(5X,’ ae=’,G15.7,’ MeV
ue=’,G15.7,’ MeV’)
write(6,190) ap(j),up(j)FORMAT(5X,’ ap=’,G15.7,’ MeV
up=’,G15.7,’ MeV’,/)end do! ------------------------------------------------------! Print media and cutoff energies assigned to each region! ------------------------------------------------------do i=1,nreg
if (med(i) .eq. 0) thenwrite(6,200) i
200
FORMAT(’ medium(’,I3,’)=vacuum’)
else
write(6,210) i,(media(ii,med(i)),ii=1,24),ecut(i),pcut(i)210FORMAT(’ medium(’,I3,’)=’,24A1,*’ecut=’,G12.5,’ MeV, pcut=’,G12.5,’ MeV’)!
----------------------------------------------!
Print out energy information of K- and L-X-rays!
----------------------------------------------if (iedgfl(i) .ne. 0) then
! Output X-ray energy
ner = nne(med(i))
do iiz=1,ner
izn = zelem(med(i),iiz) ! Atomic number of this element
write(6,220) izn
220
FORMAT(’
X-ray information for Z=’,I3)
write(6,230) (ekx(ii,izn),ii=1,10)
230
FORMAT(’
K-X-ray energy in keV’,/,*
4G15.5,/,4G15.5,/,2G15.5)
write(6,240) (elx1(ii,izn),ii=1,8)
240
FORMAT(’
L-1 X-ray in keV’,/,4G15.5,/,4G15.5)
write(6,250) (elx2(ii,izn),ii=1,5)
250
FORMAT(’
L-2 X-ray in keV’,/,5G15.5)
write(6,260) (elx3(ii,izn),ii=1,7)260
FORMAT(’
L-3 X-ray in keV’,/,4G15.5,/,3G15.5)
end doend ifend ifend dowrite(39,fmt="(’MSTA’)")write(39,fmt="(i4)") nreg
write(39,fmt="(15i4)") (med(i),i=1,nreg)
write(39,fmt="(’MEND’)")!----------------------------------------------------------------------! Step 6: Initialization-for-howfar!----------------------------------------------------------------------!----------------------------------------------------------------------! Step 7: Initialization-for-ausgab!----------------------------------------------------------------------ncount = 0
ilines = 0
nwrite = 10
nlines = 10
idin = -1
totke = 0.
wtsum = 0.
iwatch=0
!
!
270

!

=========================
call ecnsv1(0,nreg,totke)
call ntally(0,nreg)=========================
write(6,270)
FORMAT(//,’ Energy/Coordinates/Direction cosines/etc.’,/,*
6X,’e’,14X,’x’,14X,’y’,14X,’z’,*
14X,’u’,14X,’v’,14X,’w’,11X,’iq’,3X,’ir’,1X,’iarg’,/)
Energy bin width

28

deltae=ekein / 50!

Zero the variables
do i=1,5dose(i)=0.D0
doses(i)=0.D0
dose2s(i)=0.D0
end dodosep=0.D0
doseps=0.D0
dosep2s=0.D0
do j=1,50
spec(j)=0.D0
end do!

Set histories
ncases=10000000
Set maximum number for pictmaxpict=50!

tt=etime(tarray)tt0=tarray(1)!----------------------------------------------------------------------! Step 8: Shower-call!----------------------------------------------------------------------!
Write batch number
write(39,fmt="(’0
1’)")
!
!

========================
if(iwatch.gt.0) call swatch(-99,iwatch)========================! ------------------------! Start of shower call-loop! -------------------------

do i=1,ncases
latchi=0!
!
!

---------------------Select incident energy---------------------wtin = 1.0
wtsum = wtsum + wtin
etot = ekein + iabs(iqin)*RM
if(iqin.eq.1) thenavailke = ekein + 2.0*RM
else
availke = ekein
end if! Keep running sum of weights! Incident total energy (MeV)! Available K.E. (MeV) in system! for positron! Available K.E. (MeV) in system! for photon and electron

totke = totke + availke! Keep running sum of KE!
!
!

---------------------Select incident angle----------------------!
!
!

--------------------------------------------------Print first NWRITE or NLINES, whichever comes first--------------------------------------------------if (ncount .le. nwrite .and. ilines .le. nlines) thenilines = ilines + 1
write(6,280) etot,xin,yin,zin,uin,vin,win,iqin,irin,idin
FORMAT(7G15.7,3I5)
end if

280
!
!
!
1

!
!
!

----------------------------------------------------------Compare maximum energy of material data and incident energy----------------------------------------------------------if(etot+(1-iabs(iqin))*RM.gt.emaxe) thenwrite(6,fmt="(’ Stopped in MAIN.’,’ (Incident kinetic energy + RM) > min(UE,UP+RM).’)")
stop
end if---------------------------------------------------Verify the normalization of source direction cosines---------------------------------------------------if(abs(uin*uin+vin*vin+win*win-1.0).gt.1.e-6) thenwrite(6,fmt="(’ Following source direction cosines are not’,

29

1

!

’ normalized.’,3e12.5)")uin,vin,win
stop
end if!

==========================================================
call shower (iqin,etot,xin,yin,zin,uin,vin,win,irin,wtin)==========================================================!
!
!

---------------------------Sum variable and its square.
---------------------------do ie=1,5
doses(ie)=doses(ie)+dose(ie)dose2s(ie)=dose2s(ie)+dose(ie)*dose(ie)dose(ie)=0.d0
end dodoseps=doseps+dosep
dosep2s=dosep2s+dosep*dosep
dosep=0.d0
do ie=1,50
specs(ie)=specs(ie)+spec(ie)spec2s(ie)=spec2s(ie)+spec(ie)*spec(ie)spec(ie)=0.D0
end doncount = ncount + 1!
!

! Count total number of actual cases========================
if(iwatch.gt.0) call swatch(-1,iwatch)========================! ----------------------! End of CALL SHOWER loop
! -----------------------

end do!
!

========================
if(iwatch.gt.0) call swatch(-88,iwatch)========================
call plotxyz(99,0,0,0.D0,0.D0,0.D0,0.D0,0,0.D0,0.D0)
write(39,fmt="(’9’)")300! Set end of batch for CG View

tt=etime(tarray)tt1=tarray(1)cputime=tt1-tt0
write(6,300) cputime
format(’ Elapsed Time (sec)=’,G15.5)!----------------------------------------------------------------------! Step 9: Output-of-results!----------------------------------------------------------------------write(6,310) ncount,ncases,totke
310
FORMAT(/,’ Ncount=’,I10,’ (actual cases run)’,/,*’ Ncases=’,I10,’ (number of cases requested)’,/,*’ TotKE =’,G15.5,’ (total KE (MeV) in run)’)

320

330

if (totke .le. 0.D0) thenwrite(6,320) totke,availke,ncount
FORMAT(//,’ Stopped in MAIN with TotKE=’,G15.5,/,*’ AvailKE=’,G15.5, /,’ Ncount=’,I10)
stop
end iftdet=10.0
rdet=5.0
tcase=0.2
rcase=0.2
write(6,330) tdet,rdet,tcase,rcase
FORMAT(/’ Detector length=’,G15.5,’ cm’/*
’ Detector radius=’,G15.5,’ cm’/*
’ Polyetylene case thickness=’,G10.2,’ cm’/*
’ Polyetylene case side thickness=’,G10.2,’ cm’/)340

write(6,340) ekein
FORMAT(’ Results for ’,G15.5,’MeV photon’/)350

write(6,350) (media(i,2),i=1,24)
FORMAT(’ Cover material is ’,24A1/)

30

360

write(6,360) (media(i,1),i=1,24)
FORMAT(’ Inside chamber is ’,24A1/)

365

write(6,365) zin
FORMAT(’ Cs-137 exists at

!
!
!
!
!
!
!

370
!
!
!
380

z=’,F8.2,’ cm, x=y=0.0’/)----------------------------------Calculate average and its deviation----------------------------------Absorbed energy (MeV) devided by weight of air inside chamber (dweight)
MeV/g
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
tdet=10.0
rdet=5.0
vol=pi*rdet*rdet*tdet! volume of air cm3
weight=vol*rhom(med(1))! weight of air inside chamber g! rhom(1) is a density of air (medium 1)
write(6,370) vol,weight
FORMAT(’ Volume of chamber : ’,G15.5,’ cm3’/*’ weight of chamber : ’,G15.5,’ g’)--------------Absorbed dose--------------write(6,380)
FORMAT(/’ Result from absorbed energy in air inside chamber’)write(6,’(A)’) ’ Air absorbed dose from ’do ie=1,5
avdose = doses(ie)/ncountsigdose=dose2s(ie)/ncountsigdose=dsqrt((sigdose-avdose*avdose)/ncount)
avdose = avdose/weight*1.602E-10
sigdose = sigdose/weight*1.602E-10
if(ie.eq.1) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron produced in air =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.2) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from front wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.3) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from bottom wall =’,*avdose,’+-’,sigdose,’ Gy/incident photon’
elseif(ie.eq.4) thenwrite(6,’(A,G15.5,A,G15.4,A)’)*’ Electron from side wall=’,*
avdose,’+-’,sigdose,’ Gy/incident photon’
else
write(6,’(A,G15.5,A,G15.4,A)’)*’ Total=’,*avdose,’+-’,sigdose,’ Gy/incident photon’
end ifend do

400

!
!
!
!
!
!
!

write(6,400)
FORMAT(/’ Result from absorbed energy in air’/*’using track length of photon inside chamber and’,*’ mass energy absorption coefficient of air’)Track lenght (cm) gives photon fluence deviding by volume
Conversion from air absorbed dose in MeV cm2/g to that in Gy
Unit of mass energy absorption coefficient mu_en is cm2/g
Ausgab scores energy (MeV) times mu_en in unit of MeV cm2/g.
Main routine converts this into Gy (J/kg).
1MeV=1.602E-13J, 1kg=1000g
1MeV/g=1.602E-13(J/MeV)*1000(g/kg)=1.602E-10 Gy
avdose = doseps/ncountdosep2s=dosep2s/ncountsigdose=dsqrt((dosep2s-avdose*avdose)/ncount)
avdose = avdose/vol*1.602E-10
sigdose = sigdose/vol*1.602E-10

31

410
!
!
!
420

write(6,410) avdose,sigdose
FORMAT(’ Air absorbed dose =’,G15.4,’+-’,G15.4,*’ Gy/source photon’)---------------------------------------Average photon spectrum inside chamber.
---------------------------------------write(6,420)
FORMAT(/’ Average photon spectrum inside chamber’/)do ie=1,50
elow=deltae*(ie-1)
eup=deltae*ie
avspe = specs(ie)/ncountspec2s(ie)=spec2s(ie)/ncountsigspe=dsqrt((spec2s(ie)-avspe*avspe)/ncount)avspe=avspe/vol/deltae
sigspe= sigspe/vol/deltae

430

write(6,430) eup,avspe,sigspe
FORMAT(’ Upper energy’, G15.5,’ MeV--’,G12.5,’+-’,G12.5,*’ photons/cm2/MeV per source photon’)end donlist=1!
!
!
!

=============================
call ecnsv1(nlist,nreg,totke)
call ntally(nlist,nreg)=================================================
call counters_out(1)====================
stop
end

!-------------------------last line of main code-----------------------!-------------------------------ausgab.f-------------------------------! Version:
080708-1600! Reference: SLAC-265 (p.19-20, Appendix 2)!----------------------------------------------------------------------!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12! ---------------------------------------------------------------------! Required subroutine for use with the EGS5 Code System! ---------------------------------------------------------------------! A AUSGAB to:!
!
1) Score energy deposition!
2) Score particle information enter to detector from outside!
3) Print out particle transport information!
4) call plotxyz if imode=0! ---------------------------------------------------------------------subroutine ausgab(iarg)
implicit none
include ’include/egs5_h.f’
include ’include/egs5_epcont.f’include ’include/egs5_misc.f’include ’include/egs5_stack.f’
include ’include/egs5_useful.f’! Main EGS "header" file
! COMMONs required by EGS5 code

include ’auxcommons/aux_h.f’! Auxiliary-code "header" fileinclude ’auxcommons/etaly1.f’include ’auxcommons/lines.f’include ’auxcommons/ntaly1.f’include ’auxcommons/watch.f’! Auxiliary-code COMMONs

common/totals/! Variables to score32

* dose(5),dosep,deltae,spec(50),maxpict
real*8 dose,dosep,deltae,spec
integer maxpict
integer* iarg! Arguments

real*8* edepwt,dcon,decon,encoea! Local variables

integer* ie,iql,irl!
!
!

-----------------------Set some local variables-----------------------irl = ir(np)iql = iq(np)edepwt = edep*wt(np)if(iql.eq.-1.and.irl.ne.irold.and.irold.eq.2) thenif(z(np).le.0.0) thenlatch(np)=1
elseif(z(np).ge.10.0) thenlatch(np)=2
else
latch(np)=3
end ifend if!
!
!

----------------------------------------------------------Keep track of energy deposition (for conservation purposes)----------------------------------------------------------if (iarg .lt. 5) thenesum(iql+2,irl,iarg+1) = esum(iql+2,irl,iarg+1) + edepwt
nsum(iql+2,irl,iarg+1) = nsum(iql+2,irl,iarg+1) + 1
end if!
!
!
!

----------------------------------------------------------------Print out particle transport information (if switch is turned on)----------------------------------------------------------------========================
if (iwatch .gt. 0) call swatch(iarg,iwatch)========================
if(iarg .ge. 5) return!

if (irl .eq. 1) then ! inside ion chamberdose(5) = dose(5) + edepwt
if(latch(np).ge.0.and.latch(np).le.3) thendose(latch(np)+1) = dose(latch(np)+1) + edepwt
end if!
!
!

-------------------------------------------------------Score particle information if it proceeds inside chamber-------------------------------------------------------if (iql .eq. 0) then! photon
dcon=encoea(e(np))
dosep=dosep+wt(np)*e(np)*dcon*ustep
ie = e(np)/deltae +1
if(ie .gt. 50) ie = 50
spec(ie) = spec(ie) + wt(np)*e(np)*dcon*ustep
end ifend if!
!
!

---------------------------------------------------------------Print out stack information (for limited number cases and lines)---------------------------------------------------------------if (ncount .le. nwrite .and. ilines .le. nlines) thenilines = ilines + 1
write(6,100) e(np),x(np),y(np),z(np),u(np),v(np),w(np),*iql,irl,iarg
100
FORMAT(7G15.7,3I5)
end if!
!
!

-----------------------------------Output particle information for plot-----------------------------------if (ncount.le.maxpict) thencall plotxyz(iarg,np,iq(np),x(np),y(np),z(np),e(np),ir(np),33*wt(np),time(np))
end ifreturnend!--------------------------last line of ausgab.f-----------------------!-------------------------------encoea.f-------------------------------! Version:
030831-1300! Reference: SLAC-265 (p.19-20, Appendix 2)!----------------------------------------------------------------------!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12!
!
double precision function encoea(energy)!
Function to evaluate the energy absorption coefficient of air.
!
(Tables and Graphs oh photon mass attenuation coefficients and!
energy-absorption coefficients for photon energies 1 keV to!
20 MeV for elements Z=1 to 92 and some dosimetric materials,!
S. M. Seltzer and J. H. Hubbell 1995, Japanese Society of!
Radiological Technology)!---------------------------------------------------------------------double precision function encoea(energy)
real hnu(38)/0.001,0.0015,0.002,0.003,0.0032029,0.0032029,*
0.004,0.005,0.006,0.008,0.01,0.015,0.02,0.03,0.04,*
0.05,0.06,0.08,0.10,0.15,0.2,0.3,0.4,0.5,0.6,0.8,1.0,*
1.25,1.5,2.0,3.0,4.0,5.0,6.0,8.0,10.0,15.0,20.0/
real enmu(38)/3599., 1188., 526.2, 161.4, 133.0, 146.0,*
76.36, 39.31, 22.70, 9.446, 4.742, 1.334, 0.5389,*
0.1537,0.06833,0.04098,0.03041,0.02407,0.02325,0.02496,*
0.02672,0.02872,0.02949,0.02966,0.02953,0.02882,0.02789,*
0.02666,0.02547,0.02345,0.02057,0.01870,0.01740,0.01647,*
0.01525,0.01450,0.01353,0.01311/;
real*8 energy,enm1,hnu1,ene0,slope;
integer iif (energy.gt.hnu(38)) thenencoea=enmu(38)returnend ifif (energy.lt.hnu(1)) thenencoea=enmu(1)returnend ifdo i=1,38
if(energy.ge.hnu(i).and.energy.lt.hnu(i+1)) thenenm1=alog(enmu(i+1))
enm0=alog(enmu(i))
hnu1=alog(hnu(i+1))
hnu0=alog(hnu(i))
ene0=dlog(energy)
slope=(enm1-enm0)/(hnu1-hnu0)
encoea=exp(enm0+slope*(ene0-hnu0))
returnend ifif(energy.eq.hnu(i+1)) thenencoea=enmu(i+1)
returnend ifend do! If sort/interpolation cannot be made, indicate so by writing! a comment and stopping here.
write(6,100) energy
100
FORMAT(///,’ *****STOPPED IN ENCOEA*****’,/,’ E=’,G15.5,///)
returnend!--------------------------last line of encoea.f-----------------------!-------------------------------howfar.f-------------------------------! Version:
070627-1600! Reference: T. Torii and T. Sugita, "Development of PRESTA-CG! Incorporating Combinatorial Geometry in EGS4/PRESTA", JNC TN1410 2002-201,! Japan Nuclear Cycle Development Institute (2002).

34

! Improved version is provided by T. Sugita. 7/27/2004!----------------------------------------------------------------------!23456789|123456789|123456789|123456789|123456789|123456789|123456789|12! ---------------------------------------------------------------------! Required (geometry) subroutine for use with the EGS5 Code System! ---------------------------------------------------------------------! This is a CG-HOWFAR.
! ---------------------------------------------------------------------subroutine howfar
implicit nonecinclude ’include/egs5_h.f’! Main EGS "header" fileinclude ’include/egs5_epcont.f’ ! COMMONs required by EGS5 code
include ’include/egs5_stack.f’
include ’auxcommons/geom_common.f’ ! geom-common fileccinteger i,j,jjj,ir_np,nozone,jty,kno
integer irnear,irnext,irlold,irlfg,itvlfg,ihitcg
double precision xidd,yidd,zidd,x_np,y_np,z_np,u_np,v_np,w_np
double precision tval,tval0,tval00,tval10,tvalmn,delhow
double precision atvaltmp
integer iq_npc
ir_np = ir(np)iq_np = iq(np) + 2cif(ir_np.le.0) thenwrite(6,*) ’Stopped in howfar with ir(np) <=0’
stop
end ifc
if(ir_np.gt.izonin) thenwrite(6,*) ’Stopped in howfar with ir(np) > izonin’
stop
end ifc
if(ir_np.EQ.izonin) thenidisc=1
returnend ifc
tval=1.d+30
itvalm=0ccbody check
u_np=u(np)v_np=v(np)w_np=w(np)x_np=x(np)y_np=y(np)z_np=z(np)

c

c

c

c

c

c

do i=1,nbbody(ir_np)
nozone=ABS(nbzone(i,ir_np))
jty=itblty(nozone)
kno=itblno(nozone)
rpp checkif(jty.eq.ityknd(1)) thenif(kno.le.0.or.kno.gt.irppin) go to 190
call rppcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
sph checkelseif(jty.eq.ityknd(2)) thenif(kno.le.0.or.kno.gt.isphin) go to 190
call sphcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
rcc checkelseif(jty.eq.ityknd(3)) thenif(kno.le.0.or.kno.gt.irccin) go to 190
call rcccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
trc checkelseif(jty.eq.ityknd(4)) thenif(kno.le.0.or.kno.gt.itrcin) go to 190
call trccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
tor checkelseif(jty.eq.ityknd(5)) then

35

c

c

c

c

c

c

c

c

c

if(kno.le.0.or.kno.gt.itorin) go to 190
call torcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
rec checkelseif(jty.eq.ityknd(6)) thenif(kno.le.0.or.kno.gt.irecin) go to 190
call reccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
ell checkelseif(jty.eq.ityknd(7)) thenif(kno.le.0.or.kno.gt.iellin) go to 190
call ellcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
wed checkelseif(jty.eq.ityknd(8)) thenif(kno.le.0.or.kno.gt.iwedin) go to 190
call wedcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
box checkelseif(jty.eq.ityknd(9)) thenif(kno.le.0.or.kno.gt.iboxin) go to 190
call boxcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
arb checkelseif(jty.eq.ityknd(10)) thenif(kno.le.0.or.kno.gt.iarbin) go to 190
call arbcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
hex checkelseif(jty.eq.ityknd(11)) thenif(kno.le.0.or.kno.gt.ihexin) go to 190
call hexcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
haf checkelseif(jty.eq.ityknd(12)) thenif(kno.le.0.or.kno.gt.ihafin) go to 190
call hafcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
tec checkelseif(jty.eq.ityknd(13)) thenif(kno.le.0.or.kno.gt.itecin) go to 190
call teccg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)
gel checkelseif(jty.eq.ityknd(14)) thenif(kno.le.0.or.kno.gt.igelin) go to 190
call gelcg1(kno,x_np,y_np,z_np,u_np,v_np,w_np)cc**** add new geometry in herecend if
190 continue
end docirnear=ir_np
if(itvalm.eq.0) thentval0=cgeps1
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
310continue
if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) goto 320
tval0=tval0*10.d0
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
go to 310320
continuecwrite(*,*) ’srzone:1’
call srzone(xidd,yidd,zidd,iq_np,ir_np,irnext)cif(irnext.ne.ir_np) thentval=0.0d0
irnear=irnext
else
tval00=0.0d0
tval10=10.0d0*tval0
irlold=ir_np
irlfg=0330
continue
if(irlfg.eq.1) go to 340
tval00=tval00+tval10
if(tval00.gt.1.0d+06) thenwrite(6,9000) iq(np),ir(np),x(np),y(np),z(np),&u(np),v(np),w(np),tval00
9000 format(’ TVAL00 ERROR : iq,ir,x,y,z,u,v,w,tval=’,&2I3,1P7E12.5)

36

340

stop
end ifxidd=x_np+tval00*u_np
yidd=y_np+tval00*v_np
zidd=z_np+tval00*w_np
call srzold(xidd,yidd,zidd,irlold,irlfg)
go to 330
continuectval=tval00
do j=1,10xidd=x_np+tval00*u_np
yidd=y_np+tval00*v_np
zidd=z_np+tval00*w_npcwrite(*,*) ’srzone:2’
call srzone(xidd,yidd,zidd,iq_np,irlold,irnext)
if(irnext.ne.irlold) thentval=tval00
irnear=irnext
end iftval00=tval00-tval0end doif(ir_np.eq.irnear) thenwrite(0,*) ’ir(np),tval=’,ir_np,tval
end ifend ifelse
do j=1,itvalm-1
do i=j+1,itvalm
if(atval(i).lt.atval(j)) thenatvaltmp=atval(i)atval(i)=atval(j)atval(j)=atvaltmp
endif
enddo
enddo
itvlfg=0
tvalmn=tval
do jjj=1,itvalm
if(tvalmn.gt.atval(jjj)) thentvalmn=atval(jjj)
end ifdelhow=cgeps2
tval0=atval(jjj)+delhow
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
410
continue
if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) go to 420
delhow=delhow*10.d0
tval0=atval(jjj)+delhow
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
go to 410
420
continuecwrite(*,*) ’srzone:3’
call srzone(xidd,yidd,zidd,iq_np,ir_np,irnext)
if((irnext.ne.ir_np.or.atval(jjj).ge.1.).and.
&tval.gt.atval(jjj)) THEN
tval=atval(jjj)
irnear=irnext
itvlfg=1
goto 425
end ifend do425
continueif(itvlfg.eq.0) thentval0=cgmnst
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
430
continue
if(x_np.ne.xidd.or.y_np.ne.yidd.or.z_np.ne.zidd) go to 440
tval0=tval0*10.d0
xidd=x_np+tval0*u_np
yidd=y_np+tval0*v_np
zidd=z_np+tval0*w_np
go to 430
440
continue37

if(tvalmn.gt.tval0) thentval=tvalmn
else
tval=tval0
end ifend ifend ifihitcg=0
if(tval.le.ustep) thenustep=tval
ihitcg=1
end ifif(ihitcg.eq.1) THEN
if(irnear.eq.0) THEN
write(6,9200) iq(np),ir(np),x(np),y(np),z(np),&u(np),v(np),w(np),tval
9200 format(’ TVAL ERROR : iq,ir,x,y,z,u,v,w,tval=’,2I3,1P7E12.5)
idisc=1
itverr=itverr+1
if(itverr.ge.100) thenstop
end ifreturnend ifirnew=irnear
if(irnew.ne.ir_np) thencall rstnxt(iq_np,ir_np,irnew)
endif
end ifreturnend!--------------------last line of subroutine howfar---------------------

38