CG体系を用いたEGS5入門

基礎から検出器応答まで

2025-9-2 IWASE, SUGIHARA

元となるファイル

以下の三つが必要

simple.f

simple.inp

simple.data

(+ egs5フォルダ一式)

フォルダ構成

egs5/

      VERSION auxcommons/ cgview-307/ docs/ egs5run* extra_ucodes/ pegs/ run5again* slac730.pdf uccg-210721/
      auxcode/ calc/ data/ egs/ pegscommons/ samplecodes/ tutorcodes/

egs5/calc/

      egs5run*  simple/

egs5/calc/simple/

      simple.f  simple.inp  simple.data

simple.f

CG体系を用いた最小構成プログラム

メインルーチン(主プログラム)と2つのサブルーチン(副プログラム)からなる

program simple

subroutine ausgab

subroutine howfar

注意!

スライド中ソースのコピペに際して

スライド中にあるソースはコピーされる前提で載せていますが

ページ内に多めの情報を表示させるため、fortranソース行頭スペースが排除されたもの多数です

一方EGS5は基本fortran77書式のため、各行頭に7つスペースが必要です

自分のソースへペーストする際は、必要に応じて行頭7スペースを追加願います

ただしこのスペース領域に番号や継続行記号が入ることがあります


!     ~~~~~~~ Pegs COPTIONs INPUT ~~~~~~~
      IRAYL  = 1    ! Rayleigh scattering
 100  format(a3, i5, 3f9.5, 2x, 3f9.5, 2x, f9.5)
      write(*,*)'----------------- something comments -----------------
     $----------------'

上から、コメント文字、通常文、制御番号、継続行記号の例です

7スペースと、これらの制御記号に注意してコピペして下さい

注意2!

本スライドはキーボードの左キー・右キーでページめくりできます

(または、マウススクロール強め）

本スライドの右上にページ番号を表示しています

ページをリロードすると、ページ表示がズレます

ページ番号がおかしい場合、表紙まで戻ってからリロードして下さい

simple.f一式のダウンロード

egs5/calc/simpleに同梱されていない場合、以下からダウンロード

simple.fの内容確認その１

  nmed = 1               ! number of medium: 物質数

  medarr(1) = 'AL      ' ! 物質名を任意で入力（24文字)
                         ! 物質番号1にALという名称の物質を定義した
  chard(1) = 7.62        ! 物質1のcharacteristic dimension

  call pegs5             ! simple.inpを元に物質毎のテーブルを作成

  iphter(i) = 0 ...      ! EGSオプション

  iqin = 0 ...           ! 線源入力

simple.fの内容確認その2


call rluxinit         ! Initialize the Ranlux random-number generator

! source
      iqin  = 0             ! 線種：ガンマ線
      ekein = 1.253         ! 線源エネルギー
      xin   =  0.0          ! 線源座標 (x,y,z)
      yin   =  0.0
      zin   = -5.0
      uin   =  0.0          ! 方向ベクトル (u,v,w)
      vin   =  0.0
      win   =  1.0
      irin  =  0            ! initial region (0: Automatic search in CG)
      wtin  = 1.0           ! Weight = 1 since no variance reduction used

      if(irin.le.0.or.irin.gt.nreg) then

simple.fの内容確認その3

  ncases = 1000          ! 計算回数 (線源発生回数) 

  maxpict = 50           ! CGVIEW用に出力する線源数 

  call shower            ! EGSの実効的メインルーチンである shower の call

simple.inp の内容

.inpファイル = 物質情報と、物質ごとの計算パラメータの設定


ELEM
 &INP IRAYL=1 /END
AL                            AL               
AL
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 /END
PWLF
 &INP  /END
DECK
 &INP  /END

&INPの行: PEGSオプションの指定行、この場合IRAYL=1

最初のAL:物質名(任意の名称:左詰24文字)、simple.fで引用される

2番目のAL:密度効果フラグ(31文字目から入力)、

3番目のAL:元素記号(左詰、2文字+スペース)

AE:電子陽電子の下限運動エネルギーに0.511を足した値(MeV)

AP:光子の下限エネルギー(MeV)

UE:電子陽電子の上限運動エネルギーに0.511を足した値(MeV)

UP:光子の上限エネルギー(MeV)

PWLF, DECK は入れておく

上記例は電子・光子とも 0.01 MeV 〜 2 MeVまでの範囲で計算可能とするデータテーブル作成を意味する

計算体系

simple.dataの内容

simple.data

  RCC    1     0 0 0     0 0 5    10
              底面z=0     5cm高   10cm半径
  RCC    2     0 0 -10   0 0 20   15
              底面z=-10  10cm高   15cm半径
  END
   Z1       +1
   Z2       +2      -1
   Z3       -2  ← 一番最後の領域
	        =「計算打ち切り領域」
		に自動設定される
  END
   1    0   0
   ↑ z1,z2,z3
   それぞれの物質番号
   0は真空

円筒を二つ定義, RCC 1 と RCC 2

RCC形状の入力１：底面中心座標(x,y,z)

RCC形状の入力２：高さベクトル(x1,y1,z1)

RCC形状の入力３：半径(r) (cm)

計算を実行してみましょう

egs5の実行

(Win)

$ .\egs5run simple\simple
   1 file(s) copied.
   1 file(s) copied.

(Mac/Linux)

$ ./egs5run simple
============================
egs5run script has starte
============================
working directory is /Users/iwase/egs5/e/egs5/calc

以下の出力が出ますがエラーではありません


  22459 |       call PWLF1(NGL,NALG,AP,UP,RMT2,EPG,ZTHRG,ZEPG,NIPG,DLOG,DEXP,AXG,B
  
      |                                                         1
Warning: Interface mismatch in dummy procedure 'xfun' at (1): 'dlog' is not a subroutine
egs5job.f:22459:62:

22459 |       call PWLF1(NGL,NALG,AP,UP,RMT2,EPG,ZTHRG,ZEPG,NIPG,DLOG,DEXP,AXG,B
      |                                                              1
Warning: Interface mismatch in dummy procedure 'xfi' at (1): 'dexp' is not a subroutine
egs5job.f:22528:68:

22528 |         call PWLF1(NGR,NALR,1.0D-3,1.0D2,1.0D-3,EPR,ZTHRR,ZEPR,NIPR,DLOG
      |                                                                    1
Warning: Interface mismatch in dummy procedure 'xfun' at (1): 'dlog' is not a subroutine
egs5job.f:22529:9:

22529 |      *  ,DEXP, AXR,BXR,100,1,AFR,BFR,RFUNS2)
      |         1
:

CGVIEWで飛跡の確認

「ファイル」「体系・飛跡データ読込」から egs5/calc/egs5job.picを開く

飛跡の確認

z=-5 から光子が50本まっすぐ円柱１に入射し、5cmのアルミ内で半数程度が散乱している

散乱せず透過した光子は軸上を直進している

エクセサイズ

simple.fにてmaxpictを50から大きめの数字にしてみる

simple.f: 150行目あたり


!     =========================
      call ecnsv1(0,nreg,totke)
      call ntally(0,nreg)
!     =========================

      ncases=10000
      maxpict=500

      write(39,fmt="('0    1')")
      if(iwatch.gt.0) call swatch(-99,iwatch)


      do i=1,ncases                         

        wtin = 1.0

→ 実行してCGVIEWで飛跡を観察してみましょう

光子500本を5cmアルミに入射したときの光子の散乱がプロットされた

maxpictの値に応じてpictファイルサイズが肥大化します

フラックス計算に必要な情報を書き出す

surface crossingの出力

z1からz2に放出される粒子を出力

irold=1 かつ ir=2 かつ z=アルミ厚のとき

注意!

領域の呼び方問題

EGSで体系を記述するとき(CG体系）、領域をzone 1 や zone 2 と定義する

プログラム内では領域を region とか ir と呼んでいる

z2 も ir = 2 も同じ意味で、領域2を示しています

(備考)

PHITSでも入力で cell と定義したり reg と呼んだりする

両方とも、領域を指しています

補足:粒子情報について

各粒子の情報は以下のように変数化

iq: 粒子識別 ( 0:光子, -1:電子, 1:陽電子)

x,y,z: 位置座標

e: 全エネルギー

u,v,w: 方向ベクトル

irold: ステップ前の領域番号

ir: ステップ後の領域番号

実際には２次粒子３次粒子と複数粒子を扱うため配列化

np: 粒子番号

iq(np): 粒子識別( 0:光子, -1:電子, 1:陽電子)

x(np),y(np),z(np): 位置座標

e(np): 全エネルギー

u(np),v(np),w(np): 方向ベクトル

irold(np): ステップ前の領域番号

ir(np): ステップ後の領域番号

surface crossing 設定例

領域1から領域2へ移動する粒子

simple.f:39行目あたり

real*8 rnnow,etot
integer i,idin,ifti,ifto, nlist,j,k,n
character*24 medarr(MXMED)
real*8 ekine

simple.f: 246行目あたり

 if (iwatch .gt. 0) call swatch(iarg,iwatch)

      if (irold.eq.1 .and. ir(np).eq.2 .and. z(np).eq.5.0d0 )then

         if(iq(np).eq.0)then
            ekine=e(np)
         else
            ekine=e(np)-RM
         endif
         write(6,100)iq(np),ekine
 100     format(i10,f10.3)
      endif

      if(iarg .ge. 5) return

結果

egs5job.outを閲覧

egs5/calcにおいて

(Win)

$ type egs5job.out

(Mac/Linux)

$ less egs5job.out

結果 (背景色をつけた部分, 前ページの write(6,100)文に相当)

   RCC                1  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  0.0000E+00  5.0000E+00  1.0000E+01
   RCC                2  0.0000E+00  0.0000E+00 -1.0000E+01  0.0000E+00  0.0000E+00  2.0000E+01  1.5000E+01
   END       
   Z1                       1
   Z2                       2         -1
   Z3                      -2
   END       
 ranlux luxury level set by rluxgo : 1     p=  48
 ranlux initialized by rluxgo from seed           1
 RAYLEIGH DATA AVAILABLE FOR MEDIUM  1 BUT OPTION NOT REQUESTED.

 EMAXE set in HATCH to MIN(UE,UP+RM), =    2.5110E+00
 EGS SUCCESSFULLY 'HATCHED' FOR ONE MEDIUM.
         0     0.843
         0     1.253
         0     1.253
         0     0.956
         0     0.664
         0     0.081
         0     1.244
         0     0.312
         0     1.253
         0     0.646
         0     1.253
         0     1.253

検出器の応答、標的の熱計算

energy deposit量の出力

領域1に付与されたエネルギーを求める

EGS5の吸収エネルギー

光子が入射する場合

非電荷粒子である光子は物質中の移動でクーロン相互作用を伴わない

光子は物質中の飛行ではエネルギーロスがない

一方で光子は物質と相互作用すると散乱や吸収される

これらコンプトン散乱や光電吸収の際に電子が放出される

電子は電荷を持つため物質中の飛行でエネルギーを失う

1光子入射を起因として生成される電子が失ったエネルギーすべてを計上する

計算の1ステップごとのdE

ステップ毎のedep

ここで、「ある距離進む」ことをステップと呼ぶ

領域に付与されたエネルギーを dE と呼ぶ

EGSでは、ステップごとのdE値は変数 "edep"に格納

光子は最初のステップで衝突点が決定された。edep=0

衝突点で生成された電子の最初のステップで edep = e1-e1'

電子の次のステップで edep = e1'-e1''

↓

ステップごとに変化する edep を全て足し上げること

= １光子入射を起因とするdEに相当

dEをとるための変更

simple.f:39行目あたり

integer i,idin,ifti,ifto, nlist,j,k,n
character*24 medarr(MXMED)
common/cmn1/de
real*8 de

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)

simple.f:193行目あたり

uf(1)=0.0
vf(1)=0.0
wf(1)=0.0  ! Needed if lpolar(i)=1

de = 0d0

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

write(93,*)de

ncount = ncount + 1

simple.f:248行目あたり

real*8 angle,ekine
common/cmn1/de
real*8 de

if (iwatch .gt. 0) call swatch(iarg,iwatch)

if( iarg.le.4 .and. ir(np).eq.1 ) then
   de = de + edep
endif

if(iarg .ge. 5) return

IARGについて

結果

fort.93の中身

0.40966530204219881     
0.0000000000000000     
0.0000000000000000     
0.29684500780318535     
0.50804680524619394     
8.9349074513715987E-003
0.94055520763261802     
0.0000000000000000     
0.60678716436476021     
1.0392183721103767     
0.0000000000000000     
0.0000000000000000     
0.78842705839405813

(備考)gnuplotで１次元データからヒストグラムを作成する例

  set style fill solid
  filter(x,y)=floor(x/y)*y
  plot 'fort.93' u (filter($1,0.02)):(1) smooth frequency with boxes

エネルギー収支を確認する

線源エネルギー(1) = 吸収エネルギー (2) + 体系外へ逃げたエネルギー (3)

エネルギー収支の取れていない計算はどこかにミスがある可能性

メインルーチンの変更点

simple.f:39 行目あたり

real*8 totke
real*8 rnnow,etot
real*8 availke,wtin,wtsum                   
integer i,idin,ifti,ifto, nlist,j,k,n
character*24 medarr(MXMED)
common/cmn1/escore(3)
real*8 escore

open(6,FILE='egs5job.out',STATUS='unknown')

simple.f: 103行目あたり

! source
iqin  = 0   
ekein = 1d0
xin   =  0.0 
yin   =  0.0
zin   = -5.0
uin   =  0.0
vin   =  0.0
win   =  1.0
irin  =  0  
wtin  = 1.0

simple.f: 158行目あたり


ncases=10000
maxpict=50

write(39,fmt="('0    1')")
if(iwatch.gt.0) call swatch(-99,iwatch)

do i=1,3
   escore(i) = 0d0
enddo

do i=1,ncases  

escore(1) = escore(1) + ekein

simple.f:207行目あたり


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

!  write(93,*)de 

  ncount = ncount + 1         ! Count total number of actual cases

  if(iwatch.gt.0) call swatch(-1,iwatch)

end do                    ! End of CALL SHOWER loop

write(6,*)'escore results'
write(6,*)escore(1), escore(2), escore(3), escore(2) + escore(3)

ausgabの変更点

integer maxpict
real*8 angle,ekine
common/cmn1/escore(3)
real*8 escore

if (iwatch .gt. 0) call swatch(iarg,iwatch)

if( iarg.le.4 ) then
   if( ir(np).eq.1) escore(2) = escore(2) + edep
   if( ir(np).eq.3) escore(3) = escore(3) + edep
endif

CG体系を用いたEGS5入門

基礎から検出器応答まで

元となるファイル

以下の三つが必要

simple.f

simple.inp

simple.data

(+ egs5フォルダ 一式)

フォルダ構成

egs5/

egs5/calc/

egs5/calc/simple/

simple.f

CG体系を用いた最小構成プログラム

メインルーチン(主プログラム)と2つのサブルーチン(副プログラム)からなる

program simple

subroutine ausgab

subroutine howfar

注意!

スライド中ソースのコピペに際して

注意2!

simple.f一式のダウンロード

egs5/calc/simpleに同梱されていない場合、以下からダウンロード

simple.fの内容確認その１

simple.fの内容確認その2

simple.fの内容確認その3

simple.inp の内容

.inpファイル = 物質情報と、物質ごとの計算パラメータの設定

計算体系

simple.dataの内容

simple.data

計算を実行してみましょう

egs5の実行

CGVIEWで飛跡の確認

「ファイル」「体系・飛跡データ読込」から egs5/calc/egs5job.picを開く

飛跡の確認

エクセサイズ

simple.fにてmaxpictを50から大きめの数字にしてみる

→ 実行してCGVIEWで飛跡を観察してみましょう

フラックス計算に必要な情報を書き出す

surface crossingの出力

z1からz2に放出される粒子を出力 irold=1 かつ ir=2 かつ z=アルミ厚 のとき

注意!

領域の呼び方問題

補足:粒子情報について

各粒子の情報は以下のように変数化

実際には２次粒子３次粒子と複数粒子を扱うため配列化

surface crossing 設定例

領域1から領域2へ移動する粒子

結果

egs5job.outを閲覧

egs5/calcにおいて

結果 (背景色をつけた部分, 前ページの write(6,100)文に相当)

検出器の応答、標的の熱計算

energy deposit量の出力

領域1に付与されたエネルギーを求める

EGS5の吸収エネルギー

光子が入射する場合

計算の1ステップごとのdE

ステップ毎のedep

dEをとるための変更

IARGについて

結果

(備考)gnuplotで１次元データからヒストグラムを作成する例

エネルギー収支を確認する

メインルーチンの変更点

ausgabの変更点

結果

(+ egs5フォルダ一式)

z1からz2に放出される粒子を出力

irold=1 かつ ir=2 かつ z=アルミ厚のとき