Technical Material
Fortran 90/95 Code for Computing the G(E) Function
2025 Volume 24 Issue 4 Pages 99-104
Details
2025 Volume 24 Issue 4 Pages 99-104
We present a Fortran 90/95 code for computing the G(E) function, which transforms a gamma-ray pulse height spectrum to the ambient dose equivalent [H*(10)] rate, to make the G(E) function more easily accessible to any interested party. Our program code uses as input data pulse height spectra at various energies for a specific gamma-ray detector and the values of H*(10) per unit fluence at the corresponding energies, both of which are obtained with the PHITS code.
環境中の線量率は，様々な自然科学分野の研究対象であり，放射線防護の目的でも測定されている。特に2011年の東電1F事故（東京電力福島第一原子力発電所事故）によって環境中に放出された放射性同位元素は，福島県を含む広範な地域において環境中の線量率を増加させる要因になった。環境土壌に沈着した放射性同位元素は，様々なエネルギーのガンマ線を放出するため，入射ガンマ線エネルギーに応じた出力を得られる放射線測定器（検出器と容器などから構成される）が線量測定に活用されている。ただし，単一エネルギーのガンマ線に対しても放射線測定器の出力である波高スペクトルは1つのピークからなる単純な構造をしていないので，様々なエネルギーのガンマ線が存在する場合，線量率を算出するのは容易ではない。線量率を精度よく算出する手法の1つとして，G(E)関数法1）と呼ばれる方法がある。この方法は，波高スペクトルから線量率を直接算出でき，東電1F事故後の環境中での広域の線量率測定に活用されている2）。波高スペクトルの各エネルギービンに蓄えられたガンマ線のカウント数に，そのエネルギービンのエネルギーに対応するG(E)値を掛け，エネルギーについて足し合わせるだけで線量率が得られる。すなわち，手間のかかるガンマ線エネルギースペクトル解析をすることなく簡単に線量率を求めることができる。
G(E)関数は，使用する放射線測定器のレスポンスと線量換算係数から計算することができるが，用いる線量換算係数によって，周辺線量当量率（Sv/h）や空気吸収線量率（Gy/h）などが得られる。Reference 1）は，照射線量率（R/h）を対象としている。G(E)関数のエネルギー依存性は，線量換算係数のエネルギー依存性を反映するが，入射するガンマ線のエネルギーが大きくなると，検出器に吸収される割合が減るため，G(E)値は大きくなる。100 keV程度以下の低エネルギーでは，検出器を覆う容器による吸収もあり，G(E)値は大きくなるが，G(E)関数のE依存性は単純ではない。また，検出器のサイズが小さくなると，同じ線量率の環境でもガンマ線のカウント数が減るため，G(E)値は大きくなる。
さて，1960年代にJAERIの森内（Ref. 1）の著者）が開発したG(E)関数プログラムコードはFORTRAN 77で作成されているが，一般公開されておらず，JAERI/JAEA内部の利用に限られていた。東電1F事故後，様々な放射線測定器が開発されており，整備の上一般公開することは有用であると考えられる。そこで今回，Fortran 90/953）でプログラムコードを書き換え，公開することとした。この資料では，公開したコードの使い方の説明を行う。
本報では，書き換え当時に評価を行っていたGAGG（Gd3(Ga, Al)5O12(Ce)）シンチレーション測定器の場合を例にとる，大きさ1辺13 mmの立方体のGAGG結晶からなる測定器について，周辺線量当量率に対するG(E)をFORTRAN 77とFortran 90/95の両プログラムコードで計算した。2つのプログラムコードによって計算されたG(E)は，おおむね2%程度内で一致していることを確認している。FORTRAN 77プログラムコードは，この技術資料の著者の1人である津田がRef. 2）の研究においても使用した実績のあるものである。
先ずは，Ref. 1）にしたがって，G(E)関数の定義と求め方を説明する。以下は周辺線量当量率について述べるが，他の線量率でも同様である。
エネルギーEJ（keV），単位フルエンス率（1 photon/(cm2·min)）の入射ガンマ線に対する周辺線量当量率（ここではD(EJ)と記す）は，G(E)関数によって以下のように表現できるものとする。
| \begin{equation} D(E_{J}) = \int_{0}^{\infty }n (E,E_{J})G(E)dE \end{equation} | (1) |
ここで，n(E, EJ)はエネルギーEJの入射ガンマ線に対して測定器が示すエネルギーEの波高スペクトルであり，単位は[(counts/min)/keV]/[1 photon/(cm2·min)]である。また，D(EJ)，G(E)の単位は，それぞれ，(nSv/h)/[1 photon/(cm2·min)]，(nSv/h)/(counts/min)である。（1）式がG(E)の定義式であり，G(E)は使用する測定器に依存する関数である。
種々のエネルギーEJのガンマ線がそれぞれϕ(EJ)のフルエンス率で混在する場合，周辺線量当量率（ここではDと記す。単位はnSv/hである）は
| \begin{equation} D = \sum\nolimits_{J}\phi (E_{J}) \cdot D(E_{J}) \end{equation} | (2) |
| \begin{equation} D = \int_{0}^{\infty }N (E)G(E)dE \end{equation} | (3) |
となる。ここで，N(E)は
| \begin{equation} N(E) = \sum\nolimits_{J}\phi (E_{J})n(E,E_{J}) \end{equation} | (4) |
で定義される混合ガンマ線による波高スペクトルである（単位：(counts/min)/keV）。（3）式は，測定器によりN(E)が得られれば，G(E)関数を使って，周辺線量当量率を測定できることを示している。Dを時間で積算した周辺線量当量はH*(10)と記されている。
G(E)関数を求めるのに，G(E)関数を次のようにKMAX個の関数の和で展開する：
| \begin{equation} G(E) = \sum\nolimits_{K = 1}^{\textit{KMAX}}A(K)\{ g(E)\}^{K - M - 1} \end{equation} | (5) |
ここで，A(K)は1，2，…，KMAXの値をとる整数Kに依存する因子，$g$(E)は波高E（keV）の常用対数（$g$(E) = log10E），Mは，通常，0，1，2，3のいずれかの値をとるパラメターである。G(E)を求めることはKMAX個のA(K)を決定することに帰着される。
（5）式を（1）式に代入する際，波高分析器のチャンネル幅（エネルギービンの幅）をΔE，チャンネル番号をIと表すと，
| \begin{align} D(E_{J})& = \sum \nolimits_{I = 1}^{\textit{IMAX}}\Bigg[\int_{\Delta E(I - 1)}^{\Delta E \cdot I}n (E,E_{J})dE \\ &\quad\times \sum\nolimits_{K = 1}^{\textit{KMAX}}A (K)\{ g(E_{I})\}^{K - M - 1}\Bigg] \\ & = \sum\nolimits_{I = 1}^{\textit{IMAX}} \sum \nolimits_{K = 1}^{\textit{KMAX}} \Bigg[A(K)\{ g(E_{I})\}^{K - M - 1} \\ &\quad\times \int_{\Delta E(I - 1)}^{\Delta E \cdot I}n (E,E_{J})dE\Bigg] \end{align} | (6) |
ここで，$E_{I} = \Delta E\left( I - \dfrac{1}{2} \right)$，IMAXは波高スペクトルのチャンネル番号の最大値である。$\displaystyle\int_{\Delta E(I - 1)}^{\Delta E \cdot I}n(E,E_{J}) dE$をN(I, J)と置くと，これは，エネルギーEJ（keV）の単位フルエンス率の入射ガンマ線に対して測定器が示す$I$番目のチャンネルにおける計数率（単位：(counts/min)/{1 photon/(cm2·min)}）になり，（6）式は
| \begin{equation} D(E_{J}) = \sum\nolimits_{I = 1}^{\textit{IMAX}}\sum \nolimits_{K = 1}^{\textit{KMAX}} [A(K) \cdot N(I,J) \cdot \{ g(E_{I})\}^{K - M - 1}] \end{equation} | (7) |
となる。さらに，$\displaystyle\sum\nolimits_{I = 1}^{\textit{IMAX}}N (I,J) \cdot \{ g(E_{I})\}^{K - M - 1} = B(J,K)$と置くと，D(EJ)は
| \begin{equation} D(E_{J}) = \sum\nolimits_{K = 1}^{\textit{KMAX}}A (K) \cdot B(J,K) \end{equation} | (8) |
となる。
A(K)を決定するためには2組の入力データが必要になる。1つは，B(J, K)を得るためのN(I, J)であり，もう1つは，対応するエネルギーでのD(EJ)である。前者は，使用する測定器をおおむね忠実にモデル化した体系を設定し，それにエネルギーEJのガンマ線を照射した場合の計算によって求めることができる。この資料で例にとるGAGGシンチレーション測定器の場合のモデル体系をFig. 1に示す。薄すぎて図中で確認できないが，1辺13 mmのGAGG立方体結晶は0.3 mm厚のTiO2箔に覆われている。このGAGG結晶は2 mm厚のAl容器に入れられ，実際の測定器の重量と一致するように，容器内の隙間はプラスティックで充填されている。そしてこれがプラスティック容器に収納されている。われわれは，N(I, J)の計算に，JAEAが開発したPHITS（Particle and Heavy Ion Transport code System）4）を用いた（用いたPHITSのバージョンは3.250である）が，他の同等のモンテカルロ輸送計算コードを使用することも可能であり，その場合，PHITSの計算結果を読み込むプログラムコードの部分を変更すればよい。この計算例では，ガンマ線を測定器の正面に対して垂直に，真空中照射（平行ビーム）する場合を取りあげた。線源は円形で，ビームサイズが測定器全体をカバーするように半径を決めた（Fig. 1）。実際にPHITSで計算されるのは，線源から放射される1 photon当たりの計数であり，時間の単位を含まない。それに応じて，後者もD(EJ)ではなく時間を含まない単位フルエンス当たりのH*(10)/Φが必要となる。ここで，Φはフルエンスであり，フルエンス率を表すϕではない。線源の面積が単位面積（PHITSでは1 cm2）でない場合，PHITSの計算結果を単位フルエンス当たりにするためには，線源の面積を掛けなければならない。H*(10)/Φの値は，ICRP 745）中のTable A.21に，10 keVから10 MeVの範囲の25個のエネルギーに対して与えられている。このエネルギー範囲の他のエネルギーに対しては，例えば，内挿することにより求めることができる。あるいは，H*(10)の定義にしたがって体系を設定すれば，任意のエネルギーに対し，PHITSを用いて計算できる。
Model structure of the GAGG detector used in the present PHITS calculation
ところで，PHITSを用いた計算は，測定器のエネルギー分解能を考慮しないが，本来，現実の測定器の有限のエネルギー分解能を考慮すべきであると考えられる。しかし，エネルギー分解能を反映させるための畳み込み操作をしても，求めるG(E)の値はほとんど変わらなかったため，PHITSにより得られた波高スペクトルをそのまま用いることにする。
A(K)の決定には，KMAX個のEJに対する2組の入力データがあればよいが，個々のデータには誤差があり内挿精度を上げるために，最小二乗法を利用する。すなわち，以下のSの値を最小にするようなA(K)を決めることとなる。
| \begin{equation} S(J) = \frac{\displaystyle\sum \nolimits_{K = 1}^{\textit{KMAX}}A(K) \cdot B(J,K)}{D(E_{J})} - 1 \end{equation} | (9) |
| \begin{equation} S = \sqrt{\sum\nolimits_{J = 1}^{\textit{JMAX}}\{ S(J)\}^{2} } \end{equation} | (10) |
ここで，JMAXは波高スペクトルの個数であり（Jは入射ガンマ線のエネルギーに対応），JMAX > KMAXである。SをA(K)で偏微分し，0と置いた以下の（11）式から（12）式が得られる。
| \begin{align} \frac{\partial S}{\partial A(K)} &= \frac{1}{S}\sum\nolimits_{J = 1}^{\textit{JMAX}}S (J)\frac{B(J,K)}{D(E_{J})} \\&= \frac{1}{S} \sum \nolimits_{J = 1}^{\textit{JMAX}}\left[ \frac{\displaystyle\sum\nolimits_{K' = 1}^{\textit{KMAX}}A (K') \cdot B(J,K')}{D(E_{J})} - 1 \right]\frac{B(J,K)}{D(E_{J})} \\&= 0 \end{align} | (11) |
| \begin{align} &\sum\nolimits_{K' = 1}^{\textit{KMAX}}\frac{\displaystyle\sum\nolimits_{J = 1}^{\textit{JMAX}}\cfrac{B(J,K)}{D(E_{J})} \cdot \cfrac{B(J,K')}{D(E_{J})}}{\displaystyle\sum\nolimits_{J = 1}^{\textit{JMAX}}\cfrac{B(J,K)}{D(E_{J})} } A(K') \\&\quad= 1\ (K = 1, 2, \ldots , \textit{KMAX}) \end{align} | (12) |
したがって，係数が$\dfrac{\displaystyle\sum\nolimits_{J = 1}^{\textit{JMAX}}\cfrac{B(J,K)}{D(E_{J})} \cdot \cfrac{B(J,K')}{D(E_{J})}}{\displaystyle\sum\nolimits_{J = 1}^{\textit{JMAX}}\dfrac{B(J,K)}{D(E_{J})} }$，未知数がA(K′)，定数項がすべて1の連立1次方程式を解けばよいことになる。
書き換え前のFORTRAN 77プログラムコードでは，実数型変数は単精度と倍精度が混在していたが，書き換え後のFortran 90/95プログラムコードでは，組み込み関数selected_real_kindを利用し，その引数を20にして，実数型変数はすべて4倍精度にした。この引数を10とすれば，倍精度になる。
1. 入力データについて入力ファイルは2種類ある：
・ Open(10, file = ‘ファイル名’)
・ Open(ndck(j), file = ‘ファイル名’)
GAGGシンチレーション測定器に対して，エネルギーが45～7,000 keVの範囲で周辺線量当量率に対するG(E)を求めたものを例にとって説明する。
ディレクトリpulseheights-13mm_gaggに置いたinput.txtという名前を付けた装置番号10のファイルにはFig. 2に示すような文字と数字が書かれている：
Data in input.txt (see the text for the details)
3行目までは，コメント文である。GAGGは検出器の名前，13mm×13mm×13mmは検出器のサイズで，2つともformatはa20である。nSv/h/cpmとkeVは，G(E)とEの単位である。Formatは，それぞれ，a10とa4である。
4行目の5.0はA(K)を決定したあと，G(E)関数を生成する際のEの間隔であり，現在の場合，5 keVきざみということである。45.0はEJの最低エネルギー 45 keV，40.0は，おのおのの波高スペクトルにおいて計算に使用する最低のエネルギーである（40 keV）。3.5はPHITSの計算に用いた円形線源の半径（cm）である。線源が円形でない場合，面積（cm2単位）を直接入力するようにプログラムコードを変更すればよい。Formatはすべて並び入力である。
5行目の26はJMAX，すなわち用いる波高スペクトルの個数である。次の20はKMAXの最小値であり，KMAXとしては，JMAXより1小さい25から20までの値を使うことになる。1392は，EJの最低エネルギー 45 keVから始めて5 keVきざみで1,391（= 1,392 − 1）個のG(E)関数を生成することを意味し，Eの最大値は7,000 keVになる。最後の3はMの最大値で，（5）式中のMとして，値を3，2，1，0と変えて計算することを意味する。Formatはすべて並び入力である。
6行目以降の各行は，keV単位の入射ガンマ線のエネルギー，そのエネルギーでのH*(10)/Φの値である。これらの値はPHITSによる計算で得たものであるが，ICRP 745）中のTable A.21に掲載されている値とよく一致している。単位はpSv cm2であり，計算プログラム中でnSv cm2に変換する。次の整数は，PHITSで計算した波高スペクトルファイルに対する装置番号ndck(j)である。第4列はエネルギービンの幅（単位：MeV），第5列は計数が0でないチャンネル番号の最大値IMAXである。Formatはすべて並び入力である。
装置番号31から56までのファイルには，PHITSで計算した波高スペクトルのデータが収められている。例えば，装置番号31は，ディレクトリpulseheights-13mm_gagg\45に置かれた入射エネルギー45 keVに対応するdeposit3.outという名のファイルである。プログラム中，input.txt中の45.0という実数値を用いて，このディレクトリ名を作成している。他の装置番号のファイル名もすべて同じdeposit3.outで，ただ，入射ガンマ線のエネルギーに対応するディレクトリ名中の数字のみが変わる。装置番号31のファイル中にはFig. 3に示すような文字と数字が書かれている：
Part of data in deposit3.out for the incident gamma energy of 45 keV
上から12行目以下にあるemin，emax，edel，neは，それぞれ，波高スペクトルの最小エネルギー，最大エネルギー，エネルギービンの幅，エネルギービンの個数（チャンネル数）であり，edel = (emax − emin)/neの関係がある。エネルギーの単位はすべてMeVである。emin = 0としているので，neは，input.txt中の第5列にあるIMAXと同じである。また，edelはinput.txt中の第4列にある数字と同じである。プログラムは最初から文字r.errのところ（Fig. 3中の矢印が示すところ）まで読みとばし，それ以降にある4列の数字のうち，3列目の数（線源から放射されるガンマ線1個当たりの計数）のみがIMAX個読み込まれ，プログラム中の実数f(j = 1, i = 1, 2, …, IMAX)を定義する。ちなみに，4列の数字のうち1列目と2列目は，各エネルギービンの最小値と最大値であり，4列目は，計数の相対誤差である。ガンマ線の入射エネルギーが50，60，…，7,000 keVの場合は，f(j, i)のjは2, 3, …, 26と変わる。fの配列の大きさはf(30, 1024)としている。
2. 出力データについて出力ファイルは3種類ある：
・ Open (11, file = 'ファイル名')
・ Open (12, file = 'ファイル名')
・ Open (13, file = 'ファイル名')
G(E)関数の値は，装置番号11のファイルにFig. 4に示すように出力される。
Part of calculated G(E) values
装置番号12のファイルには，いろいろなKMAXとMの組み合わせで計算された（10）式のSを100倍した値（S′）が出力される。Table 1の3列目にS′の値を示す。これは実数型変数が4倍精度のものであるが，比較のため，倍精度にした場合のS′の値を第4列に示す。小さなS′の値に対応するG(E)関数を選択することが適切であるが，計算例において4倍精度の場合，KMAXとMのどの組み合わせでもS′の値は同程度であり，G(E)の値も数%の範囲で同じである。Fig. 5にKMAX = 25，M = 3の場合のG(E)関数（+，×印の2曲線のうち，+のもの）を示す（S′ = 0.908）。倍精度の場合，S′の値の小さなものについてのG(E)の値は4倍精度の場合と同程度であるが，S′の値の大きなものについてのG(E)の値は異なる。Fig. 6にKMAX = 22，M = 0の場合のG(E)関数を示す（S′ = 58.147）。緒言で，Fortran 90/95とFORTRAN 77のプログラムコードは互いにおおむね2%程度内で一致するG(E)値を与えると述べたが，Fig. 5にFORTRAN 77のプログラムコードで計算したG(E)値（×の曲線，KMAX = 25，M = 3）も示す。
| KMAX | M | S′ | |
|---|---|---|---|
| Quadruple precision |
Double precision |
||
| 25 | 3 | 0.908 | 10.610 |
| 24 | 3 | 0.909 | 4.104 |
| 23 | 3 | 0.920 | 4.067 |
| 22 | 3 | 0.919 | 3.418 |
| 21 | 3 | 0.921 | 3.407 |
| 20 | 3 | 0.923 | 3.726 |
| 25 | 2 | 0.908 | 4.589 |
| 24 | 2 | 0.908 | 3.532 |
| 23 | 2 | 0.924 | 3.705 |
| 22 | 2 | 0.926 | 3.611 |
| 21 | 2 | 0.922 | 4.158 |
| 20 | 2 | 0.922 | 3.985 |
| 25 | 1 | 0.908 | 13.254 |
| 24 | 1 | 0.908 | 4.153 |
| 23 | 1 | 0.908 | 4.072 |
| 22 | 1 | 0.915 | 4.078 |
| 21 | 1 | 0.918 | 4.186 |
| 20 | 1 | 0.925 | 4.040 |
| 25 | 0 | 0.908 | 5.778 |
| 24 | 0 | 0.908 | 6.254 |
| 23 | 0 | 0.909 | 7.774 |
| 22 | 0 | 0.908 | 58.147 |
| 21 | 0 | 0.912 | 4.887 |
| 20 | 0 | 0.941 | 5.105 |
G(E) function for KMAX = 25 and M = 3 (+ symbol is for Fortran 90/95 with quadruple precision real numbers and × symbol is for FORTRAN 77)
G(E) function for KMAX = 22 and M = 0 (double precision real numbers)
装置番号13のファイルには，いろいろなKMAXとMの組み合わせに対する（5）式中のA(K)の値が出力され，また，それらの値の場合について，（12）式の左辺がどの程度右辺の1に近いかを示している。1から多少なりともずれているので，得られたA(K)の値は真の値とはいえないが，これらのA(K)の値を使って求めたG(E)関数は真のA(K)値を使って得られるものと大きく違わないと考えている。比較のため，4倍精度実数におけるKMAX = 25，M = 3の場合と倍精度実数におけるKMAX = 22，M = 0の場合の対応する出力の1部をFig. 7に示す。
Part of obtained values of A(K) and the left side of Eq. (12) in the text
本プログラム（GE_v1.f90）はweb上6）にて無償で公開している。
環境中のガンマ線の波高スペクトルから周辺線量当量率を求めるのに用いられるG(E)関数を計算するプログラムコード（Fortran 90/95）を作成し，その使い方の説明を行った。空気吸収線量率など他の線量率についても，入力データの一部を変更することで対応できる。今回公開したプログラムコードによりG(E)関数の利用が多少なりとも容易になれば幸いである。
本技術資料は，『令和6年度原子力施設等防災対策等委託費（無人航空機へ搭載するLPWA等を活用した環境放射線モニタリング機器の実現可能性検証）事業』の成果を含む。
