egs5の線源形状の変化（等方、放射状、面の統合編（２次元））

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーションにおいて、線源の形状を変化させ、デフォルトのペンシルビームから等方線源、放射状線源、面線源のいずれかにする方法を紹介します。（egs5公式はこちら）

具体的には、既に公開してあるegs5の線源形状の変化（等方線源編、放射状線源編、面線源編）を一つにまとめた統合版を紹介します。

手順

ucnaicgv.fのSTEP8のdo文内のSelect incident angleの下に

        selectmode = 3                  !1:isotropy,2:radial,3:surface
        if(selectmode .eq. 1) then
          call randomset(rnnow)
          win = 2.0*rnnow-1.0
          call randomset(rnnow)
          phai = PI*(2.0*rnnow-1.0)
          uin = dsqrt(1.0-win*win)*cos(phai)
          vin = dsqrt(1.0-win*win)*sin(phai)
        end if
        if(selectmode .eq. 2) then
          win =cos(30*PI/180)
          call randomset(rnnow)
          win = win+rnnow*(1.0-win)
          call randomset(rnnow)
          phai = PI*(2.0*rnnow-1.0)
          uin = dsqrt(1.0-win*win)*cos(phai)
          vin = dsqrt(1.0-win*win)*sin(phai)
        end if
        if(selectmode .eq. 3) then
          call randomset(rnnow)
          xin = 6*(rnnow-0.5)
          call randomset(rnnow)
          yin = 6*(rnnow-0.5)
        end if

と入力します。
selectmodeが１のとき等方線源、２のとき放射状線源、３のとき面線源になります。

次に、使用した変数(selectmodeとphai)の定義です。
STEP1の以下の場所に

,phai

,selectmode

と入力します。

実行結果１（CGVIEWの表示：maxpict=200、selectmode=1(等方））

実行結果２（CGVIEWの表示：maxpict=200、selectmode=2(放射状））

実行結果３（CGVIEWの表示：maxpict=200、selectmode=3(面））

egs5の線源形状の変化（面線源編）

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーションにおいて、線源の形状を変化させ、デフォルトのペンシルビームから面線源にする方法を紹介します。（egs5公式はこちら）

手順

ucnaicgv.fのSTEP8のdo文内のSelect incident angleの下に

        call randomset(rnnow)
        xin = 6*(rnnow-0.5)
        call randomset(rnnow)
        yin = 6*(rnnow-0.5)

と入力します。call randomset(rnnow)はrnnowに0~1の値をランダムで格納するための呼び出しです。
ここでは6cm×6cmの正方形の面線源を定義してみました。

実行結果１（CGVIEWの表示：maxpict=200、6cm×6cm）

面線源になっていることがわかります。

見やすいようにmaxpictを50から200に変更しています。

実行結果２（CGVIEWの表示：maxpict=200、6cm×1.5cm）

egs5の線源形状の変化（放射状線源編）

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーションにおいて、線源の形状を変化させ、デフォルトのペンシルビームから放射状線源にする方法を紹介します。（egs5公式はこちら）

手順

ucnaicgv.fのSTEP8のdo文内のSelect incident angleの下に

        win =cos(30*PI/180)
        call randomset(rnnow)
        win = win+rnnow*(1.0-win)
        call randomset(rnnow)
        phai = PI*(2.0*rnnow-1.0)
        uin = dsqrt(1.0-win*win)*cos(phai)
        vin = dsqrt(1.0-win*win)*sin(phai)

と入力します。call randomset(rnnow)はrnnowに0~1の値をランダムで格納するための呼び出しです。
最初の行の30は中心軸からの角度で、任意で変更することができます。180にすれば等方線源になります。

次にSTEP1に上記で使用したphaiという変数を定義するために、zi0のあとに

,phai

と入力します。

実行結果１（CGVIEWの表示：maxpict=200、30°）

放射状線源になっていることがわかります。

見やすいようにmaxpictを50から200に変更しています。

実行結果２（CGVIEWの表示：maxpict=200、60°）

egs5の線源形状の変化（等方線源編）

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーションにおいて、線源の形状を変化させ、デフォルトのペンシルビームから等方線源（全方向(４π)に同一の確率で放射線を出す線源）にする方法を紹介します。（egs5公式はこちら）

手順

ucnaicgv.fのSTEP8のdo文内のSelect incident angleの下に

        call randomset(rnnow)
        win = 2.0*rnnow-1.0
        call randomset(rnnow)
        phai = PI*(2.0*rnnow-1.0)
        uin = dsqrt(1.0-win*win)*cos(phai)
        vin = dsqrt(1.0-win*win)*sin(phai)

と入力します。call randomset(rnnow)はrnnowに0~1の値をランダムで格納するための呼び出しです。

次にSTEP1に上記で使用したphaiという変数を定義するために、zi0のあとに

,phai

と入力します。

実行結果（CGVIEWの表示：maxpict=200）

等方線源になっていることがわかります。

見やすいようにmaxpictを50から200に変更しています。

egs5でエネルギー付与をcsvファイル出力する方法（パルス波高分布の表示）

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーション中において、光子がエネルギー付与した値などをcsvファイル出力する方法を紹介します。（egs5公式はこちら）

具体的には、エネルギー1.253MeVの光子（γ線）をucnaicgv.dataによって定義されたNaI検出器の体系に向かって放射したときに、NaIに付与されるエネルギーをcsvファイルに出力してパルス波高分布（γスペクトル）として見てみたいということです。

なお、デフォルトの『ucnaicgv』のシミュレーションで出力されるegs5job.outからPulse height distributionとしてパルス波高分布がわかりますが、光子数で割った値となっています。ここでは、横軸エネルギー(KeV)、縦軸count(計数値)のパルス波高分布を見てみたいと思います。

前準備と環境

・egs5が動作する環境があれば大丈夫です。

・OSはwindows10です。

・今回使用しているコンパイラはgfortranです。g77は試してないのでわかりませんが、可能だと思います。（できなかったらご指摘ください。）

手順１

以下のようにucnaicgvシミュレーションが動作する環境を作成してください。

『ucnaicgv.f』内のプログラムを書き換えていきます。

コマンドプロンプトの画面

フォルダと作成されたファイル

手順２

解析状況を出力するためのcsvファイルを作成するためSTEP1のファイル定義の部分に

      open(51,FILE='Pulse height distribution.csv',STATUS='replace')

と入力します。数字やファイル名は適当でいいですが、数字はあとのプログラムでも同じ数字に変更する必要があります。

手順３

作成したcsvファイルにNaIに付与されたエネルギーを出力するため、if (depe .gt. 0.D0) thenの直後に

          write(51,500) depe*1000
500       FORMAT(F12.6)

と入力します。

これによって、ausgabにより出力されたdepe(NaIへの付与エネルギー)の値が0より大きいとき、depeの値を作成したcsvファイルに出力するようになります。depe*1000としているのはMeVをkeVに変換するためです。（MeVで表示させたい場合は*1000を消して利用してください。）

実行結果の見方

実行最初に、ucnaicgv.fがあるフォルダに『Pulse height distribution.csv』というテキストファイルが作成されます。

ファイルを開くと下図のようになっていると思います。

A1を選択してShift+Ctrl+↓で全選択し、挿入のグラフのヒストグラムをクリックすることでパルス波高分布を表示することができます。（横軸のビン数を変更することができます。）

ucnaicgv.dataで作成されたNaI検出器のNaIの領域は結構大きいので、ほとんど入射するエネルギーである1253keVのcountが大きいです。また、1050~1120keVは他のエネルギーよりcountが小さいです。これは980~1050keVの間にコンプトンエッジがあるからだと考えられます。

egs5でシミュレーション解析状況を知る方法

本記事の紹介

本記事は、放射線輸送計算コードegs5の公式より配布されている『ucnaicgv』のシミュレーション中において、シミュレーションの解析状況を表示するようにしたときのものです。（egs5公式はこちら）

egs5シミュレーションにおいて、計算回数(ncases)を何億回にするなど、大きな数にすると、何時間もシミュレーションをし続けることになります。この間、デフォルトの『ucnaicgv』のシミュレーションでは、進捗度はわからないため、いつシミュレーションが終わるのかわかりません。

今回、その問題を解決するためにプログラムを改造し、下の図のようにテキストファイルにシミュレーションの解析状況を表示させるようにしました。

シミュレーション解析状況

前準備と環境

・egs5が動作する環境があれば大丈夫です。

・OSはwindows10です。

・今回使用しているコンパイラはgfortranです。g77は試してないのでわかりませんが、可能だと思います。（できなかったらご指摘ください。）

手順１

以下のようにucnaicgvシミュレーションが動作する環境を作成してください。

『ucnaicgv.f』内のプログラムを書き換えていきます。

コマンドプロンプトの画面

フォルダと作成されたファイル

手順２

作成したテキストファイルに項目を表示するため、STEP8のdo文の前に

      write(0,*) 'シミュレーション解析状況'
      write(0,*) '    進捗度  ','  計算回数  ','  時,分,秒,ミリ秒' 

と入力します。

手順３

手順２で作成した項目の内容を表示するためSTEP8のdo文の中のncount = ncount + 1の後に

        if (ncount/(ncases/100) .gt. bin) then
            call date_and_time(times(0), times(1), times(2),value)
            write(0,*) ncount/(ncases/100),'%',ncount,'     ',times(1)
            call flush(0)
            bin = bin + 1
        end if

※コンパイラによってはtimesの変数でエラーが発生することがあります。このときの対処として、

call date_and_time(times(0), times(1), times(2),value)を

call date_and_time(times(1), times(2), times(3),value)に変更

また、

write(0,*) ncount/(ncases/100),'%',ncount,'     ',times(1)を

write(0,*) ncount/(ncases/100),'%',ncount,'     ',times(2)に変更してみてください。

手順４

シミュレーションに要した時間を出力するために、STEP8の最後付近に

      write(0,299) cputime
299   format('シミュレーション時間',G15.5,'秒')

と入力します。

手順５

手順１～４で使用した変数（bin,times,value）の定義と初期化をします。

①binは整数,valueは整数配列とするため、STEP1のintegerのntypeの後に

,bin,value(8)

と入力します。

②timesは文字配列とするため、STEP1のcharacter*24 medarr(MXMED)の後に

      character*10 times(3)

と入力します。

③binの初期化のため、STEP7の最後に

      bin = 0

と入力します。

実行結果の見方

実行最初に、ucnaicgv.fがあるフォルダに『シミュレーション解析状況.txt』というテキストファイルが作成されます。

シミュレーション実行中に、進捗度が１％上がるごとに更新されます。

更新されたときの進捗度（シミュレーションが何％完了したか）、計算回数（ncases）、時間（下の図では、19時16分45.414秒）を知ることができます。

補足ですが、１％と２％の時間から１００％になる、つまりシミュレーションが完了するまで時間を計算することができます。（ここでは出力していません。）

最後にシミュレーションにかかった時間が出力されます。

補足

①bin = bin + 1の1の値を変更することで、更新頻度を変更できます。

例：bin = bin + 5のとき

※少数にするときは、定義を変更する必要があります。

②if (ncount/(ncases/100) .gt. bin) thenをif (ncount*100/ncases .gt. bin) thenとしていないのは、ncasesが億など大きいときにncount*100がintegerの値の範囲を超えてしまうことを考慮してです。

③デフォルトのプログラムのシミュレーション時間はegs5job.outによると

7.7188秒であったため、プログラム変更によるシミュレーション時間の影響は小さいかほぼないと考えられます。