対話型EGSアプリの作成

検出器編




2024-7-9 IWASE

対話型EGSアプリとは


  • 線源や計算体系の変更をパラメータ制御できるようにしたEGSアプリ(exeファイル)

  • EGSアプリを実行後、いくつかの設問(パラメータ入力)へ回答

  • 様々な条件の計算を入力ファイルの編集なしに実行できる

  • 実行頻度の少ない計算でも、ブランク復帰時の労力ゼロ





    • (参考)
  • EGSを用いた放射線教育
  • ソースコード(zip)

    • 内容



    実習A



    det.f det.inp
    det.data    による

    検出器応答計算

    の入力作成

    実習B



    appdet.f による

    対話側検出器応答計算アプリ

    の開発

    対話側EGSアプリ検出器 実習A

    検出器の応答計算



    内容

    A-1. NaIに光子を入射

    det0.f, det0.inp, det0.data

    A-2. エネルギーの集計

    det1.f, det1.inp, det1.data

    A-3. 応答関数の出力

    det4.f, det4inp, det4.data

    A-4. Geなど他の検出器を追加

    det5.f, det5inp, det5.data

    注意!

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



  • スライド中にあるソースはコピーされる前提で載せていますが
  • ページ内に多めの情報を表示させるため、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!



  • 本スライドはキーボードの左キー・右キーでページめくりできます
    • (または、マウススクロール強め）

  • 本スライドの右上にページ番号を表示しています
  • ページをリロードすると、ページ表示がズレます

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

    • 実習A-1. NaIに光子を入射 (1)

    ファイル構成
  • det0.f
  • det0.inp
  • det0.data


    • 計算体系
    体系ファイル (det0.data)
    
RCC 1 0 0 0   0 0 5   10
RCC 2 0 0 -10 0 0 20  15
END
 Z1       +1 
 Z2       +2      -1 
 Z3       -2
END
 1    0   0
    物質ファイル (det0.inp)
    
COMP
 &INP NE=2 RHO=3.667 IRAYL=1 PZ=1,1 &END
NAI                           NAI              
NA I
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END
    ユーザーコード (det0.f)の線源部分
    
!     ~~~~~~~ SOURCE INPUT ~~~~~~~
      iqin  = 0         
      ekein = 1.d0   
      xin   =  0.0      
      yin   =  0.0
      zin   = -1.0
      uin   =  0.0      
      vin   =  0.0
      win   =  1.0
      irin  =  0        
      wtin  = 1.0
    ユーザーコード (det0.f)の物質定義部分
    
!     ~~~~~~~ MATERIALs INPUT ~~~~~~~
      medarr(1) ='NAI                      '

    実習A-1. NaIに光子を入射 (2)

    入力内容を確認したら、計算を実行

    (Win) egs5run.batの冒頭部分
    @ECHO OFF
SET BASKET=..\..\egs5
SET uc=det0
    (Mac/Linux) egs5runの冒頭部分
    #!/bin/sh
BASKET="../../egs5"
MY_MACHINE=gfort
OPT_LEVEL=O
UC=det0
  • EGS5ディレクトリは本実習ディレクトリ内
  • 相対パスで指定している
  • EGS5をインストールした絶対パスで再指定可能
  • 変数UCにファイル名"det0"を先行指定
    • egs5の実行
    (Win) 
    $ ./egs5run                                                                                    
   1 file(s) copied.
   1 file(s) copied.

    (Mac/Linux) 
    $ .\egs5run                                                                                    
============================
egs5run script has started 
============================
 
working directory is /Users/iwase...

    実習A-1. NaIに光子を入射 (3)

    CGVIEWで計算結果を確認

    cgviewの起動
    Win 
    $ ..\..\egs5\cgview-307\bin\cgview.exe
    Mac/Linux 
    $ wine ../../egs5/cgview-307/bin/cgview.exe
    体系・飛跡データファイル(egs5job.pic)の読み込み
  • 「ファイル」⇒「体系・飛跡データ読込」
  • フォルダがappdet/lec-a/det0であることを確認して
  • egs5job.picを開く
  • 飛跡がプロットされたことを確認

    • det1.fのメインルーチン部を編集

    
character*24 medarr(MXMED)
integer IRAYL,IAPRIM,IB...
!integer ISSB,FPSTFL,FU...
integer IPHTE,IEDGF,IAU...
common /cmn1/maxpict,escore
real*8 escore(0:3)
  • (備考) 通常はescore(4)と定義するところをescore(0:3)とした
  • これにより、escore(0) ~ escore(3)が使用できる
    • （escore(0)を入射、escore(i)を領域iと対応付けるため) 

    
DO i=1,ncases  
  etot = ekein + iabs(iqin)*RM
 escore(0) = escore(0) + etot


    
 call shower (iqin,etot,xin,...

 ncount = ncount + 1
 if(iwatch.gt.0) call...
ENDDO

!/////////////////////////..
!     POST PROCESSING
!/////////////////////////..

write(6,*)"escore results"
write(6,*)escore(0)
write(6,*)escore(1)
write(6,*)escore(2)
write(6,*)escore(3)

    det1.fのausgabを編集

    
!//////////////////////////////////////////////////////////////////////      
      subroutine ausgab(iarg)
!//////////////////////////////////////////////////////////////////////      
      implicit none

      include 'include/egs5_h.f'         ! Main EGS "header" file

      include 'include/egs5_epcont.f'    ! COMMONs required by EGS5 code
      include 'include/egs5_misc.f'
      include 'include/egs5_stack.f'
      include 'include/egs5_useful.f'
      include 'auxcommons/aux_h.f'       ! Auxiliary-code "header" file
      include 'auxcommons/etaly1.f'      ! Auxiliary-code COMMONs
      include 'auxcommons/lines.f'
      include 'auxcommons/ntaly1.f'
      include 'auxcommons/watch.f'

integer iarg,maxpict
common/cmn1/maxpict,escore
real*8 escore(0:3)

! Print out particle transpo...
if (iwatch .gt. 0) call swatch(iarg,iwatch)

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

    計算を実行

    egs5の実行
    (Win) 
    $ ./egs5run                                                                                    
   1 file(s) copied.
   1 file(s) copied.

    (Mac/Linux) 
    $ .\egs5run                                                                                    
============================
egs5run script has started 
============================
 
working directory is /Users/iwase...


    結果を確認
     escore results
   10000.000000000000     
   4939.4334043026493     
   0.0000000000000000     
   5060.5665957998553
    エネルギーが保存されていることを確認

    実習A-3 応答関数の出力


    手順
  • 応答関数とは、各ヒストリごとの検出器内dEのヒストグラムである
  • ヒストグラムの下限値と上限値、およびビン数またはビン幅を決める
  • 各ヒストリごとに検出器内のdEを集計 (ausgab にて de = de + edep)
  • 各ヒストリの終了後にdeの値がどのビンに属しているか調べる
  • de値が属するビンが決まったら、ヒストグラムを+1 する (実際には + wt)
  • ヒストリ計算の後処理として、ヒストグラムを書き出す

    • ヒストグラムに必要な変数の追加

    det2.fのメインルーチンを変更
    
integer IPHTE,IEDGF,IAUGE,LPOLA
common /cmn1/maxpict,escore,de
real*8 escore(0:3)
real*8 de,emin,emax,ebin,e1,e2,dehist(10000)
integer nbin,ie
!     ~~~~~~~ INITIALIZATION ~~~~~~

    
open(ifto,FILE='egs5job.pic',STATUS='unknown')
open(95,FILE='dehist.out',STATUS='unknown')
call counters_out(0)  ! set zero for all counters

  • call showerの前にdeをリセット
  • 作成するヒストグラムのビン数および下限値上限値を設定
  • ビン幅を計算
  • deが値を持つとき
  • de値をビン幅で割った値に+1したものがビン番号なので
  • そのビン番号に対応するヒストグラムに１カウント加える(実際にはウエイト)
  • ヒストリ計算終了後に、ヒストグラムをファイル出力
    • 
!     ~~~~~~~ EGS5 MAIN ROUTINE ~~~~~~~
de = 0d0
call shower (iqin,etot,xin,yin,zin,uin,vin,

ncount = ncount + 1    ! Count total number
if(iwatch.gt.0) call swatch(-1,iwatch)

! HISTOGRAM dE
nbin = 256
emin = 0
emax = ekein+abs(iqin)*2*RM
ebin = (emax - emin)/real(nbin - 1)

if (de .gt. 0d0 ) then
 ie = de/ebin +1
 if (ie .gt. nbin) ie = nbin
! dehist(ie) = dehist(ie) + 1
 dehist(ie) = dehist(ie) + wtin
endif

ENDDO
!/////////////////////////////////
!     POST PROCESSING
!/////////////////////////////////
! WRITE histogram
e1 = emin
e2 = e1 + ebin
do i=1,nbin+1
write(95,*) e1, e2, dehist(i)/real(ncases)
e1 = e2
e2 = e1 + ebin
enddo
write(6,*)"escore results"
    det2.fのausgabを変更
    

integer iarg,maxpict
common/cmn1/maxpict,escore,de
real*8 escore(0:3)
real*8 de

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

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

if(iarg .ge. 5) return
    dehist.outを確認
    
   0.0000000000000000        3.9215686274509803E-003   1.4000000000000000E-003
   3.9215686274509803E-003   7.8431372549019607E-003   1.4000000000000000E-003
   7.8431372549019607E-003   1.1764705882352941E-002   1.8000000000000000E-003
   1.1764705882352941E-002   1.5686274509803921E-002   1.4000000000000000E-003
   1.5686274509803921E-002   1.9607843137254902E-002   2.3999999999999998E-003
:

  • dehist.out は e1, e2, dehist で書き出されている
  • 応答関数のデータが得られた
  • dehist.outをグラフソフトでプロットすると次のようになる

    • エクセサイズ

  • det2.fにおいて、ekein の値を 1 MeV から変化させて計算し
  • 異なる入射エネルギーに対する応答関数をプロットする


  • 例えば 1 MeV, 2 MeV, 0.5 MeV で計算し
  • グラフソフトでプロットすると次のようになる

    • 実習A-4 Geなど他の検出器を追加

  • NaI以外のシンチレータや半導体検出器を使えるようにする
  • 検出器の外側領域を空気とする


    • 物質ファイル(.inpファイル)を編集
  • det4.inp ⇒ NaIを使用するために必要な入力
  • 使用したい物質だけ追記する

  • Ge半導体検出器を追加してみましょう
    • det3.fの中身 
    
COMP
 &INP NE=2 RHO=3.667 IRAYL=1 PZ=1,1 &END
NAI                           NAI              
NA I
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END

    det3.inpを編集

    物質にGeを使用するときの.inp 
    
ELEM
 &INP IRAYL=1 &END
GE                            GE               
GE
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END
  • 単体物質なので、"ELEM" で記述
  • 単体物質なので、組成などの入力不要
  • その他はNaIの場合と同様

  • det3.inpオリジナルにこの部分をペースト ⇒
    • det3.f(編集後) 
    
COMP
 &INP NE=2 RHO=3.667 IRAYL=1 PZ=1,1 &END
NAI                           NAI              
NA I
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END
ELEM
 &INP IRAYL=1 &END
GE                            GE
GE
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END

    体系ファイル (.data) のおさらい

    
  RCC    1     0 0 0     0 0 5    10
  RCC    2     0 0 -10   0 0 20   15
  END
   Z1       +1 
   Z2       +2      -1 
   Z3       -2
  END
   1    0   0
  • 境界1 としてRCC(円筒), 底面中心が原点、高さ5cm, 半径10cm
  • 境界2 としてRCC(円筒), 底面中心がz=-10、高さ20cm, 半径15cm

  • 領域1は 境界1の内側
  • 領域2は 境界2の内側で境界1の外側
  • 領域3は 境界2の外側

  • 領域1は物質は物質番号１
  • 領域2の物質は真空
  • 領域3の物質は真空
  • 領域1は検出器
  • 領域2は検出器外領域
  • 領域3は計算外領域

    • det3.fにおいて、medarr(1)='NAI'を'GE'へ書き換える

    
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
!  ~~~~~~~ MATERIALs INPUT ~~~~~~~
medarr(1) ='GE                       '
     
! ~~~~~~~ PEGS COPTIONs INPUT ~~~~~~~
IRAYL  = 1    ! Rayleigh scattering
IAPRIM = 1    ! Emrirical correction in brems XS
  • 検出器領域は物質１を使うことになっていた
  • 物質１をNAIからGEに変更
  • 検出器はGEに変更された

    • 計算を実行

  • Geの1 MeV光子に対する応答関数が得られた

  • ncase=10000では計算精度は不十分
  • 検出器固有の分解能再現にはdehistへ外部関数が必要

    • 空気をdet3.inpに追加

    
MIXT
 &INP RHO=1.2050d-3 NE=3 RHOZ=0.75575,0.23143,0.01282 GASP=1.0 IRAYL=1 &END
AIR AT NTP                    AIR-GAS                       
N  O  AR
ENER
 &INP AE=0.521,AP=0.010,UE=2.511,UP=2.0 &END
PWLF
 &INP  &END
DECK
 &INP  &END
  • .inpではELEM(単体), COMP(化学式)、MIXT(混合物)から選択
  • MIXTの場合、元素を重量比で入力

    • det3.fで空気を使えるようにする

    使用する物質数nmedを1⇒2へ変更 
    
! ~~~~~~~ PARAMETERs INPUT ~~~~~~~
ncases = 10000      
maxpict= 50
nmed = 2
luxlev = 1 
inseed = 1

    物質2を追加して、.inpで定義した空気を入力 
    
! ~~~~~~~ MATERIALs INPUT ~~~~~~~
medarr(1) ='GE                       '
medarr(2) ='AIR AT NTP               '

    det3.datで領域2を空気で定義

    
  RCC    1     0 0 0     0 0 5    10
  RCC    2     0 0 -10   0 0 20   15
  END
   Z1       +1 
   Z2       +2      -1 
   Z3       -2
  END
   1    2   0
    領域1は物質1, 領域2は物質2 と定義

    計算を実行

    
escore results
   10000.000000000000     
   5932.9017078206480     
  0.83643223900442398
  • Geの1 MeV光子に対する応答関数が得られた

    • (備考) ncase=10000では計算精度は不十分
    検出器固有の分解能再現にはdehistへ外部関数が必要

    その他の検出器を追加



  • det3.inpに以下の物質を追加
    • BGO 
    
COMP
&INP RHO=7.13 NE= 3 PZ=4,3,12 IRAYL=1 &END
BGO                           BGO                     
BI GE O 
ENER
&INP AE=0.521 AP=0.01 UE=1.521 UP=1 &END
PWLF
  &INP  &END
DECK
  &INP  &END
    CSI 
    
COMP
&INP RHO=4.51 NE= 2 PZ=1,1 IRAYL=  1 &END
CSI                           CSI                     
I  CS
ENER
&INP AE=0.521 AP=0.01 UE=1.521 UP=1 &END
PWLF
  &INP  &END
DECK
  &INP  &END
    PbWO4 
    
COMP
&INP RHO=8.24 NE= 3 PZ=4,1,1 IRAYL=  1 &END
PBWO4                         PBWO4                   
O  W  PB
ENER
&INP AE=0.521 AP=0.01 UE=1.521 UP=1 &END
PWLF
  &INP  &END
DECK
  &INP  &END
    Plastic Scint. 
    
MIXT
&INP RHO=1.032 NE=2 RHOZ=0.085,0.915 IRAYL=1 &END
PLASTIC SCINT.                PLASTIC SCINT.          
H  C 
ENER
&INP AE=0.521 AP=0.01 UE=1.521 UP=1 &END
PWLF
  &INP  &END
DECK
  &INP  &END

    medarr(1)='GE'を'別の検出器'へ書き換える

    
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
!  ~~~~~~~ MATERIALs INPUT ~~~~~~~
medarr(1) ='PLASTIC SCINT.           '
     
! ~~~~~~~ PEGS COPTIONs INPUT ~~~~~~~
IRAYL  = 1    ! Rayleigh scattering
IAPRIM = 1    ! Emrirical correction in brems XS
  • 物質１をGEからPLASTIC SCINT.(プラシン)に変更
  • 検出器はプラシンに変更された

    • 計算を実行


  • プラシンの1 MeV光子に対する応答関数が得られた

    • 以上、実習A


  • 検出器の応答計算を行う入力を作成
  • det3.f, det3,inp, det3.data
  • リンク先はこちら (ファイルUP後リンクします)

    • 対話側EGSアプリ検出器 実習B

    「検出器の応答」対話型アプリ



    内容

    B-1. det3.fにdet3.inpとdet3.dataを組み入れる

    det.f

    B-2. いくつかの変数をターミナル入力へ変更

    appdet.f

    実習B-1-1 det3.fにdet3.dataを組み入れる

  • egs5run内部で.dataファイルをegs5job.inpと改名して扱っている
  • egs5job.inpのファイル識別番号は4
    • 
open(99,FILE='egs5job.out',STATUS='unknown')
open(4,FILE='egs5job.inp')
open(25,FILE='pgs5job.pegs5inp')
open(39,FILE='egs5job.pic',STATUS='unknown')
open(95,FILE='dehist.out',STATUS='unknown')
call counters_out(0)  ! set zero for all counters
  • det.fに下を追記
  • 下記の内容を4番ファイル(上でegs5job.inpと定義)にwrite
    • 
ifti = 4
! Write CG input data
write(ifti,*)"RCC", 1, 0.,0.,0., 0.,0.,5, 10
write(ifti,*)"RCC", 2, 0.,0.,-10, 0.,0., 20, 15
write(ifti,*)"END"
write(ifti,*)"Z1 +1"
write(ifti,*)"Z2 +2 -1"
write(ifti,*)"Z3 -2"
write(ifti,*)"END"
write(ifti,*)"1 2 0"
100 format(a3, i5, 3f9.5, 2x, 3f9.5, 2x, f9.5)
rewind(ifti)

    実習B-1-2 det3.fにdet3.inpを組み入れる

    方針検討

  • det3.inpを眺めていると、ENERより上の部分と下の部分で分けられそうである
  • ELEM, COMP, MIXT の定義部分でひとかたまり
  • ENERの定義部分でひとかたまり、と考えてみる

  • 前半部は、物質の定義とオプション定義
  • 後半部は、エネルギー上限と下限値の定義

  • それぞれを変数化し、整形してファイル出力するようプログラムを書いてみる

    • NaIで書いてみる検討

    こう変数化しておいて 
    
subst(1)='COMP'
medarr(1)='NAI                     '
denarr(1)=medarr(1)
rho(1) = 3.667 
ne(1) = 2
pz(1,1) = 1.0  
pz(1,2) = 1.0
elem(1,1)='NA'
elem(1,2)='I '

    
AE = ecut(1)+RM
AP = pcut(1)
UE = ekein*1.01+RM  
UP = ekein
    こうファイルへ書きこんでいけば、 .inp が作れそうである 
    
write(25,1)"&INP"
write(25,1,advance='no')"RHO="
write(25,2) rho(m)
if(subst(m).eq.'COMP')then
write(25,1,advance='no')"NE="
write(25,3) ne(m)
write(25,1,advance='no')"PZ="
write(25,2)(pz(m,i),i=1,ne(m))
elseif(subst(m).eq.'MIXT')then
write(25,1,advance='no')"NE="
write(25,3) ne(m)
write(25,1,advance='no')"RHOZ="
write(25,2)(rhoz(m,i),i=1,ne(m))

    
write(25,1)"ENER"
write(25,1)"&INP"
write(25,1,advance='no')"AE="
write(25,2)AE
write(25,1,advance='no')"AP="
write(25,2)AP
write(25,1,advance='no')"UE="
write(25,2)UE
write(25,1,advance='no')"UP="
write(25,2)UP

    書いてみる


    左がその結果で、右は左のファイルの改行やスペースを省いたもの (ekein=20MeVの場合) 
    
COMP
&INP
RHO=3.6670E+00
NE=  2
PZ=1.0000E+00 1.0000E+00
&END
NAI                           NAI                     
NA I 
ENER
&INP
AE=5.2100E-01
AP=1.0000E-02
UE=2.0711E+01
UP=2.0000E+01
&END

    
COMP
&INP RHO=3.667 NE=2 PZ=1,1 &END
NAI                           NAI                     
NA I 
ENER
&INP AE=0.521 AP=0.01 UE=20.711 UP=20 &END
    &INP &END間は namelist というフォートランの機能で
    カンマやスペースや複数行など寛容性の高い入力形式
    左の情報でもEGS5動作を確認

    PEGSパラメータ取り扱いの検討

    以下のように変数化しておいて 
    
! ~~~~~~~ Pegs COPTIONs INPUT ~~~~~~~
      IRAYL  = 1  
      IAPRIM = 1  
      IBOUND = 1  
      INCOH  = 0  
      ICPROF = -3 
      IMPACT = 0
    以下のように途中へ追記すればよさそうである 
    
write(25,1,advance='no')"IRAYL="
write(25,3)IRAYL
write(25,1,advance='no')"IAPRIM="
write(25,3)IAPRIM
write(25,1,advance='no')"IBOUND="
write(25,3)IBOUND
write(25,1,advance='no')"INCOH="
write(25,3)INCOH 
write(25,1,advance='no')"ICPROF="
write(25,3)ICPROF
write(25,1,advance='no')"IMPACT="
write(25,3)IMPACT

    その結果


    
COMP
&INP
RHO=3.6670E+00
NE=  2
PZ=1.0000E+00 1.0000E+00
IRAYL=  1
IAPRIM=  1
IBOUND=  1
INCOH=  0
ICPROF= -3
IMPACT=  0
&END
NAI                           NAI                     
NA I 
ENER
&INP
AE=5.2100E-01
AP=1.0000E-02
UE=2.0711E+01
UP=2.0000E+01
&END
    PEGSパラメータがCOMPの&INP と &ENDの中に組み込まれた

    .inp内製を実装

    
real*8 totke,rnnow,etot,...
character*24 medarr(MXMED)
integer NELEM
parameter(NELEM = 20)
character*24 denarr(MXMED)
character*2  elem(MXMED,NELEM)
character*4  subst(MXMED)   ! substance (ELEM, COMP, or MIXT)
integer IRAYL,IAPRIM,IBOUND,INCOH,ICPROF,IMPACT,IGSDST

    
integer nbin,ie
real*8  rho(MXMED),gasp(MXMED),pz(MXMED,NELEM),rhoz(MXMED,NELEM)
integer ne(MXMED)
! ~~~~~~~ INITIALIZATION ~~~~~~
npreci=3                  ! PICT data mode for CGView in free format
do m=1,MXMED
  gasp(m) = 0d0          
  ne(m)   = 1
enddo

    
      write(ifti,*)"1 2 0"
 100  format(a3, i5, 3f9.5, 2x, 3f9.5, 2x, f9.5)
      rewind(ifti)

      subst(1)='COMP'
      medarr(1)='NAI                     '
      denarr(1)=medarr(1)
      rho(1) = 3.667        ! DENSITY (g/cm^3)
      ne(1) = 2
      pz(1,1) = 1.0         ! ATOM FRACTIONs
      pz(1,2) = 1.0
      elem(1,1)='NA'
      elem(1,2)='I '

      subst(2)='MIXT'
      medarr(2) ='AIR AT NTP              '
      denarr(2) ='AIR-GAS                 '
      rho(2) = 1.2050d-3
      gasp(2) = 1.0
      ne(2) = 3
      rhoz(2,1) = 0.75575 
      rhoz(2,2) = 0.23143
      rhoz(2,3) = 0.01282
      elem(2,1)='N '
      elem(2,2)='O '
      elem(2,3)='AR'

!/////////////////////////////////////////////////////////////
!     PREPARE CALCULATION
!/////////////////////////////////////////////////////////////
      AE = ecut(1)+RM
      AP = pcut(1)
      UE = ekein*1.01+RM  ! just e+RM is ideal but when e >> 1, RM may disapear in PEGS input
      UP = ekein

!     ~~~~~~~ FOR PEGS ~~~~~~~
      do m = 1, nmed
         write(25,1)subst(m)
         write(25,1)"&INP"
         if(subst(m).eq.'COMP')then
            write(25,1,advance='no')"RHO="
            write(25,2) rho(m)
            write(25,1,advance='no')"NE="
            write(25,3) ne(m)
            write(25,1,advance='no')"PZ="
            write(25,2)(pz(m,i),i=1,ne(m))
         elseif(subst(m).eq.'MIXT')then
            write(25,1,advance='no')"RHO="
            write(25,2) rho(m)
            write(25,1,advance='no')"NE="
            write(25,3) ne(m)
            write(25,1,advance='no')"RHOZ="
            write(25,2)(rhoz(m,i),i=1,ne(m))
         endif
         if(gasp(m).ne.0)then
            write(25,1,advance='no')"GASP="
            write(25,2) gasp(m)
         endif
         write(25,1,advance='no')"IRAYL="
         write(25,3)IRAYL
         write(25,1,advance='no')"IAPRIM="
         write(25,3)IAPRIM
         write(25,1,advance='no')"IBOUND="
         write(25,3)IBOUND
         write(25,1,advance='no')"INCOH="
         write(25,3)INCOH 
         write(25,1,advance='no')"ICPROF="
         write(25,3)ICPROF
         write(25,1,advance='no')"IMPACT="
         write(25,3)IMPACT
         write(25,1)"&END"

         write(25,4)adjustl(medarr(m)),adjustl(denarr(m))
         write(25,5)(elem(m,i),i=1,ne(m))

         write(25,1)"ENER"
         write(25,1)"&INP"
         write(25,1,advance='no')"AE="
         write(25,2)AE
         write(25,1,advance='no')"AP="
         write(25,2)AP
         write(25,1,advance='no')"UE="
         write(25,2)UE
         write(25,1,advance='no')"UP="
         write(25,2)UP
         write(25,1)"&END"
         write(25,1)"PWLF"
         write(25,*)" &INP  &END"
         write(25,1)"DECK"
         write(25,*)" &INP  &END"
      enddo
      
 1       format(a)
 2       format(50(1pe10.4,x))
 3       format(i3)
 4       format(a24,6x,a24)
 5       format(50(a2,x))

      close(25)
! ~~~~~~~ PREPARE CALCULATION ~~~~~~~
      if(nmed.gt.MXMED) then

    実行する



    注意!
  • .inpファイルと.dataファイルは内製するので必要ないはず
  • しかしegs5run内で実行前のファイルチェックがあるためエラーとなる
  • .inpと.dataは空ファイルでよいので作成しておく


  • ここまでで実習B-1 .inpと.dataの組入れの基本は完成
  • ELEMの入力例も次に示す

    • .inp内製でELEMを入力する例

    
         subst(1)='ELEM'
         medarr(1)='GE                      '
         denarr(1)=medarr(1)
         elem(1,1)='GE'
    COMPやMIXTと比べるとシンプル

    以上.inpファイルと.dataファイルを内製するサンプルプログラムは こちら(リンク後日)

    実習B-2 対話型アプリの開発

  • 実習B-2-1 体系のパラメータ化と対話入力

  • 実習B-2-2 物質の対話入力

    • 実習B-2-1 体系のパラメータ化と対話入力



    計算体系のいくつかの寸法を図のようにパラメータ化すると


    計算体系は次のようにかける 
    
      b = 1.0 + d + a + 1.0
      s = r + 1.0
      write(ifti,*)"RCC", 1, 0.,0.,0.,     0.,0., a,  r
      write(ifti,*)"RCC", 2, 0.,0.,-(d+1), 0.,0., b,  s
    このとき、対話型アプリで入力させるパラメータはa,r,dとなる

    ソースを変更

    
      integer ne(MXMED)
      real*8 d,a,r,b,s ! Geometry
!     ~~~~~~~ INITIALIZATION ~~~~~~


    
! ~~~~~~~ GEOMETRY INPUT ~~~~~~~
      write(*,*)'//////////////////////////////'
      write(*,*)'////// PROGRAM APPDET ////////'
      write(*,*)'//////////////////////////////'
      write(*,*)''
      write(*,*)'RADIUS of the detector (cm)?'
      read(*,*)r
      write(*,*)'LENGTH of the detector (cm)?'
      read(*,*)a
      write(*,*)'DISTANCE between source and deterctor (cm)?'
      read(*,*)d

      b = 1.0 + d + a + 1.0
      s = r + 1.0

      write(ifti,*)"RCC", 1, 0.,0.,0., 0.,0.,a, r
      write(ifti,*)"RCC", 2, 0.,0.,-(d-1), 0.,0., b, s
      write(ifti,*)"END"

    編集できたら実行してみましょう

    実習B-2-2 物質の対話入力



    実習Aで準備した検出器をプログラム内で準備し、対話選択できるようにする

    次のように編集

    
      integer ne(MXMED)
      real*8 d,a,r,b,s ! Geometry
      integer nd       ! Detector
!     ~~~~~~~ INITIALIZATION ~~~~~~

    
!     ~~~~~~~ PEGS (MATERIAL) INPUT ~~~~~~~
      write(*,*)'which DETECTOR?'
      write(*,*)'1:NaI, 2:Ge:, 3:BGO, 4:CsI, 5:PbWO4'
      write(*,*)'6:Pl.Scint., 7:GLASS,LEAD'
      read(*,*)nd

      if     (nd.eq.1) then
         subst(1)='COMP'
         medarr(1)='NAI                     '
         denarr(1)=medarr(1)
         rho(1) = 3.667        ! DENSITY (g/cm^3)
         ne(1) = 2
         pz(1,1) = 1.0         ! ATOM FRACTIONs
         pz(1,2) = 1.0
         elem(1,1)='NA'
         elem(1,2)='I '
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.2) then
         subst(1)='ELEM'
         medarr(1)='GE                      '
         denarr(1)=medarr(1)
         elem(1,1)='GE'
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.3) then
         subst(1)='COMP'
         medarr(1)='BGO                     '
         denarr(1)=medarr(1)
         rho(1) = 7.13         ! DENSITY (g/cm^3)
         ne(1) = 3
         pz(1,1) =  4.         ! ATOM FRACTIONs
         pz(1,2) =  3.
         pz(1,3) = 12.
         elem(1,1)='BI'
         elem(1,2)='GE'
         elem(1,3)='O '
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.4) then
         subst(1)='COMP'
         medarr(1)='CSI                     '
         denarr(1)=medarr(1)
         rho(1) = 4.51         ! DENSITY (g/cm^3)
         ne(1) = 2
         pz(1,1) =  1.         ! ATOM FRACTIONs
         pz(1,2) =  1.
         elem(1,1)='I '
         elem(1,2)='CS'
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.5) then
         subst(1)='COMP'
         medarr(1)='PBWO4                   '
         denarr(1)=medarr(1)
         rho(1) = 8.24         ! DENSITY (g/cm^3)
         ne(1) = 3
         pz(1,1) =  4.         ! ATOM FRACTIONs
         pz(1,2) =  1.
         pz(1,3) =  1.
         elem(1,1)='O '
         elem(1,2)='W '
         elem(1,3)='PB'
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.6) then
         subst(1)='MIXT'
         medarr(1)='PLASTIC SCINT.          '
         denarr(1)=medarr(1)
         rho(1) = 1.032        ! DENSITY (g/cm^3)
         ne(1) = 2
         rhoz(1,1) =  0.085    ! WEIGHT FRACTIONs
         rhoz(1,2) =  0.915
         elem(1,1)='H '
         elem(1,2)='C '
         chard(1) = a          ! STEP-length autoset (cm)
      elseif (nd.eq.7) then
         subst(1)='MIXT'
         medarr(1)='GlASS, LEAD             '
         denarr(1)=medarr(1)
         rho(1) = 6.22         ! DENSITY (g/cm^3)
         ne(1) = 5
         rhoz(1,1) =  0.156453 ! WEIGHT FRACTIONs
         rhoz(1,2) =  0.080866
         rhoz(1,3) =  0.008092
         rhoz(1,4) =  0.002651
         rhoz(1,5) =  0.751938
         elem(1,1)='O '
         elem(1,2)='SI'
         elem(1,3)='TI'
         elem(1,4)='AS'
         elem(1,5)='PB'
         chard(1) = a          ! STEP-length autoset (cm)
      endif

      subst(2)='MIXT'
      medarr(2) ='AIR AT NTP              '
      denarr(2) ='AIR-GAS                 '

    入射粒子の入力


    入射粒子種やエネルギーも対話入力できるようにする 
    
      elem(2,3)='AR'
      chard(2) = 30

!     ~~~~~~~ COMMAND ~~~~~~~
      write(*,*)'which RADIATION?'
      write(*,*)'0:Photon, -1:Electron, 1:Positron'
      read(*,*)iqin
      write(*,*)'kinetic ENERGY in MeV?'
      read(*,*)ekein

!/////////////////////////////////////////////////////////////
!     PREPARE CALCULATION
!/////////////////////////////////////////////////////////////

      AE = ecut(1)+RM
      AP = pcut(1)

    計算回数の入力


    計算回数 (ncases) も対話入力する 
    
!//////////////////////////////////////////////////////////////////////
!     SHOWER MAIN-LOOP
!//////////////////////////////////////////////////////////////////////      
      write(*,*)''
      write(*,*)'NUMBER of calculation (ncases)?'
      read(*,*)ncases

      DO i=1,ncases  

         etot = ekein + iabs(iqin)*RM      ! Incident total energy (MeV)
         escore(0) = escore(0) + etot      ! Score indicent energy (MeV)

    対話式EGSアプリ「検出器」


  • ここまでで一通りのことができるようになりました
  • ファイルはこちら（後日アップします）


  • これをベースにアルミケースの追加や、等方線源の追加などに挑戦ください


  • 実習は以上です、お疲れさまでした
    • gfortran for windows links
    tdm-gcc webpage
    exe file direct link