PHITSを用いたフォトンカウンティングCTシステムのシミュレーション

令 和 元 年 度 修 士 論 文

PHITS を用いた

フォトンカウンティング CT システムの

シミュレーション

指導教員 櫻井 浩 教授

群馬大学大学院理工学府 理工学専攻

電子情報・数理教育プログラム

阿久津 和音

目次

1 序論 ... 1 1.1 研究背景 ... 1 1.2 フォトンカウンティング CT システムの開発 ... 2 1.3 線減弱係数の測定結果 ... 6 1.4 バックグラウンドスペクトル𝑰_𝑩𝑮(𝑬)の測定結果 ... 8 1.5 シミュレーションに使用するソフトウェアについて ... 9 2 本研究の目的 ... 11 3 シミュレーションシステムの構築 ... 12 3.1 シミュレーション条件 ... 12 3.2 ジオメトリ ... 13 3.3 単エネルギーX 線の散乱シミュレーション ... 15 4 シミュレーションによるバックグラウンド要因の検討 ... 16 4.1 コリメータ穴に関する検証 ... 16 4.2 検出器カバー材質に関する検証 ... 17 4.3 サンプルによる散乱に関する検証 ... 19 4.4 パイプの効果に関する検証 ... 20 5 シミュレーションによるバックグラウンド再現の試み ... 21 5.1 10[keV]ごとの散乱シミュレーション ... 21 5.2 重み付け加算 ... 24 5.3 シミュレーションによるバックグラウンド ... 25 5.4 実測バックグラウンドとの比較 ... 26 5.5 コンプトン散乱の検討 ... 28 5.6 真空雰囲気でのバックグラウンドシミュレーション ... 32 6 まとめ ... 35 7 今後の課題 ... 36 8 参考文献 ... 37 9 謝辞 ... 38 10 補足 ... 39 10.1 ジオメトリの作成について ... 39 10.2 PHITS 入力データについて ... 42 10.3 参考文献（補足） ... 51

1

1 序論

1.1 研究背景

放射線治療では、放射線を照射する位置や線量を決定するための治療計画 が、治療の良し悪しを決定する重要な手順である。治療計画では放射線が通過 する組織の電子密度分布が重要となる。これは、放射線（X 線や荷電粒子）が 物質内でエネルギーを損失する主なメカニズムが電磁相互作用（光電効果や 弾性・非弾性散乱など）であることによる。このことから、組織内の電子密度 の情報が治療計画に重要であることが分かる。 以前より、生体内の電子密度分布と実効原子番号を得る試みとして、２つの エネルギーの X 線を用いた CT スキャン（Dual-Energy X-ray CT (DXCT)） が提唱されてきた。この方法を更に敷衍すれば、CT スキャンにおいて X 線光 子一つ一つのエネルギーを識別することにより、電子密度や実効原子番号を 定量的に求めることが可能である。この定量測定の基礎は、X 線に対する物質 （被写体）の減弱係数を定量性よく測定することにある。近年 Photon Counting CT（PCCT）と呼ばれる装置が医療現場に展開されつつあり、それ を実現できる素地が整えられつつある。 当研究室では PCCT に着目し、実験室規模の CT システムを用いて様々な 物質の減弱係数を高い精度で測定しようと研究を進めてきた。より高い精度 で減弱係数を測定できれば、治療計画に使用する CT の精確さが向上し、最終 的に重粒子線治療の精度向上に繋がるためである。

2

1.2 フォトンカウンティング CT システムの開発

当研究室にて開発したフォトンカウンティング CT システムの全体図を、 次に示す。システム全体は鉛ハッチの内部に配置されており、一台の PC で X 線源、検出器、精密ステージを全てコントロール可能となっている。測定した スペクトルは、補正ソフトにより補正を加える。 図 1.1 PCCT システムの概略図 図 1.2 鉛ハッチ写真

3 X 線管球は浜松ホトニクス製マイクロフォーカス X 線源（L12161-07）[3]_で ある。L12161-07 の仕様を表 1.1 に示す。また、実験によって得た 1 次スペク トルを、図 1.3 に示す。 表 1.1 L12161-07 の仕様[4] 図 1.3 L12161-07 の 1 次スペクトル測定結果 0 500000 1000000 1500000 2000000 2500000 3000000 3500000 4000000 4500000 0 50 100 150 200 250 Inte nsit y [C o unts] Energy [keV]

4

また、検出器は Amptek 社製の CdTe 検出器（XR-100T-CdTe）[5]_である。 XR-100T-CdTe の仕様および検出効率を次に示す。

表 1.2 XR-100T-CdTe の仕様（一部抜粋）[5]

図 1.4 XR-100T-CdTe の相互作用確率の対数プロット[5] (1 keV ～ 1 MeV)

5

図 1.5 XR-100T-CdTe の相互作用確率の線形プロット[5] (10 keV ～ 250 keV)

6

1.3 線減弱係数の測定結果

図 1.1 の CT 装置を用い、小林が線源弱係数の測定実験を行った[2]_。その結 果を次に示す。なお、理論値は NIST (アメリカ国立標準技術研究所)が開発し た元素データベースである XCOM[6]_{を利用した。} 図 1.6 小林による線源弱係数の測定結果（試料：Al, H2O,C, Mg）[2] また、理論値との差を検討するために、実験値/理論値の計算を行った結果 を次に示す。値が１に近いほど誤差が小さいことを示している。 図 1.7 小林による理論値との差の検討[2]

7 図 1.6 より、実験値は理論値と比較して 2～10％程度小さいことが分かっ た。特に高エネルギー側で差が大きくなる傾向にあることも分かった。 我々は線源弱係数が一致しない原因として、『背景放射線（バックグラウンド） の寄与』と『多重散乱の寄与』の２つを考えた。特に、『バックグラウンドの 寄与』について、説明する。 線源弱係数は、下記の式にて説明できる。なお、𝜇(E)は線源弱係数、𝐼(E)は 物質通過後の X 線強度、𝐼₀(𝐸)は照射された初期の X 線強度である。

𝜇(𝐸) = −ln⁡(

𝐼(𝐸) 𝐼0(𝐸)

)

(1) 我々は、これにバックグラウンド成分𝐼𝐵𝐺(𝐸)が含まれているのではないかと 考えた。すなわち、式(1)は、次の式(2)のようになると考えられる。

𝜇(𝐸) = −ln⁡(

𝐼(𝐸)+𝐼𝐵𝐺(𝐸) 𝐼0(E)+𝐼𝐵𝐺(𝐸)

)

(2)

8

1.4 バックグラウンドスペクトル𝑰

_𝑩𝑮

(𝑬)の測定結果

線源弱係数の理論値との不一致の理由をバックグラウンドにあると仮定し、 実際にバックグラウンドスペクトル

𝐼

_𝐵𝐺

(𝐸)

の測定を行った。その結果を、次 に示す。 図 1.8 バックグラウンドスペクトル

_𝐼

_𝐵𝐺

_(𝐸)

の測定結果 図 1.8 は、Pb-𝐾_𝛽線の強度が Pb-𝐾_𝛼線よりも大きくなっていることを示して いる。通常は Pb-𝐾𝛼線のほうが Pb-𝐾𝛽線よりも強度が大きいはずである。こ の原因は、Pb-𝐾𝛽線近傍何らかの散乱が重なっているためであると推測される。 私はこの散乱について、シミュレーションでより詳細に原因を検討できな いかと考えた。更に、CT 装置のシミュレーションシステムを構築することに より、将来的には遮蔽計算や装置設計、データの分析等に活かせるのではない かと考え、CT 装置の散乱シミュレーションに関する研究を行うこととした。

9

1.5 シミュレーションに使用するソフトウェアについて

シミュレーションには放射線輸送計算に特化したモンテカルロ計算コード である PHITS（Particle and Heavy Ion Transport code System）[7]_{を利用する} こととした。PHITS は、任意の体系中における様々な放射線の挙動を各反応 モデルや核データを用いて模擬するモンテカルロ計算コードである。JAEA （日本原子力研究開発機構）、RIST（高度情報科学技術研究機構）、KEK（高 エネルギー加速器研究機構）、九州大学など、日本の研究機関が中心となって 開発が行われている。加速器遮蔽設計や放射線治療研究、放射線防護研究、宇 宙・地球惑星科学など幅広い分野で利用されている。 図 1.9 PHITS について[8]

10 PHITS を利用する利点としては、シミュレーションプログラムの作成が簡 単である点が挙げられる。ユーザーはシミュレーション条件やジオメトリ等 を記載した入力データを用意すればよい。それを PHITS に読み込ませるだけ で、自動でシミュレーションが実行される。また、シミュレーション結果は数 値のみならず自動でグラフ化される。これは PHITS にグラフィックソフトウ ェア ANGEL が同梱されているためである。また、日本の研究機関が中心と なって開発されたため、日本語でのサポートが充実している点も利点の一つ である。以上より、PHITS は他の放射線輸送計算コードと比較して、初心者 に対する敷居が比較的低いと言える。 PHITS と他の放射線輸送計算コードとの比較を、次に示す。 表 1.3 PHITS と他の放射線輸送計算コードとの比較

11

2 本研究の目的

本研究の目的は、PHITS を用いたシミュレーションシステムの構築するこ と、および、シミュレーションを用いて線源弱係数が理論値と一致しない原因 を検討することである。後者については、特にバックグラウンドの検討を行う こととし、その第一歩として、まずはバックグラウンドをシミュレーションに おいて再現、その評価を行うこととした。

12

3 シミュレーションシステムの構築

3.1 シミュレーション条件

シミュレーション条件は、次のようにした。

13

3.2 ジオメトリ

ジオメトリは、実際の CT 装置を再現したものを作製した。概略図を次に示 す。 図 3.1 実際の CT 装置を再現したジオメトリ（横から見た図） シミュレーションはジオメトリ全体を空気雰囲気として行った。 また、タリー（PHITS における検出器のようなもの）は検出器カバー内部 全体に指定した。寸法は表 3.1 のとおりである。 実際の実験で使用している検出器の検出素子は CdTe であるため、実際の 検出素子に合わせた寸法の CdTe およびタリーを設定することが望ましい。 だが、今回は簡便のために、タリーを上記のように設定し、内部での散乱の影 響を最小限に抑えるために、タリー内部を真空雰囲気とした。 そ れ ぞ れ の 材 質 の 密 度 情 報 は 、 次 の よ う に し た 。 空 気 に つ い て は (窒素) : (酸素) ＝ 8 : 2 の混合気体とした。 表 3.2 それぞれの材質の密度

14

線源から射出される X 線は、次のような頂角 43 度のコーンビームである。

15

3.3 単エネルギーX 線の散乱シミュレーション

初めに、入射 X 線のエネルギーを 120[keV]とし、図 3.1 のジオメトリを用 いてシミュレーションを行った。その結果を、次に示す。 図 3.3 120[keV]の X 線を入射させた場合の散乱シミュレーション （対数プロット） まず、バックグラウンドの要因、特に『どこからの散乱 X 線がバックグラ ウンドへの寄与が強いのか』について検討を行うために、ジオメトリを様々に 変更してシミュレーションを行い、その結果を図 3.3 と比較を行うこととし た。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [/cm^ 2/so urc e] Energy [keV]

16

4 シミュレーションによるバックグラウンド要因の検討

4.1 コリメータ穴に関する検証

コリメータの穴を塞げば、シミュレーションでバックグラウンドを得られ るのではないかと考えた。そこで、コリメータの穴を塞ぎ、シミュレーション を行った。その結果を次に示す。 図 4.1 コリメータ穴に関する検証（対数プロット） 図 4.1 より、コリメータの穴の有無によるスペクトルの変化はほぼ見られ ないことが分かる。その原因としては、『散乱 X 線がコリメータ（材質： W）を透過している』こと、あるいは、『散乱 X 線が検出器カバー（材質： Pb）を透過している』ことが考えられる。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [/cm^ 2/so urc e] Energy [keV] コリメータ穴あり コリメータ穴なし

17

4.2 検出器カバー材質に関する検証

前項において、『散乱 X 線がコリメータ（材質：W）を透過している』こ と、『散乱 X 線が検出器カバー（材質：Pb）を透過している』ことが示唆さ れた。私は後者の『散乱 X 線が検出器カバー（材質：Pb）を透過している』 可能性に着目した。検出器カバーの材質をタングステンに変更し、シミュレ ーションを行った。その結果を、次に示す。 図 4.2 検出器カバーに関する検証（対数プロット） 図 4.2 より、検出器カバー材質をタングステンに変更すると、検出器内部 に入る散乱 X 線が減少することが分かる。これは、検出器カバーからの透過 X 線があったことを示している。 本研究ではシミュレーションによるバックグラウンドの再現が目的である ため、以降も検出器カバーの材質は鉛のままでシミュレーションを行うこと とした。なお、実験に際しては、検出器カバーの材質をタングステンに変更 したほうがよい可能性がある。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [/cm^ 2/so urc e] Energy [keV] 検出器カバー材質：Pb 検出器カバー材質：W

18 次に、XCOM[6]_{による Pb および W の線減弱係数のエネルギー依存性を示} す。W は Pb より線減弱係数が大きく、放射線遮蔽に適していることが分か る。 図 4.3 Pb および W の線減弱係数のエネルギー依存性（対数プロット） 0.1 1 10 100 1000 0 50 100 150 200 線減弱係数 [cm^ -1] エネルギー [keV] Pb W

19

4.3 サンプルによる散乱に関する検証

次に、サンプル（図 3.1 の Al 棒）からの散乱の影響の程度を検証するため に、サンプルを排除した散乱シミュレーションを行った。その結果を次に示す。 図 4.4 サンプルによる散乱に関する検証（対数プロット） 図 4.4 より、サンプルの有無による散乱の変化はあまり見られない事が分 かる。よって、サンプルからの散乱の影響は無視できる程度のものであると 考えられる。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [a /c m ^3/sour ce ] Energy [keV] サンプルあり サンプルなし

20

4.4 パイプの効果に関する検証

次に、試料-検出器間に挿入されているパイプにより、どの程度散乱を減ら せているかについて検証を行うため、パイプを排除した散乱シミュレーショ ンを行った。その結果を次に示す。 図 4.5 パイプの効果に関する検証（対数プロット） 図 4.5 より、パイプの有無によるスペクトルの変化がみられたものの、傾向 としては軽微なものであることが分かる。よって、バックグラウンドには、パ イプの外側を通る散乱 X 線の寄与が大きいと推測される。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [/cm^ 2/so urc e] Energy [keV] パイプなし パイプあり

21

5 シミュレーションによるバックグラウンド再現の試み

5.1 10[keV]ごとの散乱シミュレーション

バックグラウンドを再現するにあたって、まず、散乱 X 線を 30, 40, 50, 60, 70, 80, 90, 100, 110, 120[keV]と 10[keV]刻みに変更し、散乱シミュレーショ ンを行った。その結果を次に示す。

(a) E=30 [keV] (b) E=40 [keV]

(c) E=50 [keV] (d) E=60 [keV] 図 5.1 10[keV]ごとの散乱シミュレーション (1) （対数プロット） 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV]

22

(e) E=70 [keV] (f) E=80 [keV]

(g) E=90 [keV] (h) E=100 [keV]

(i) E=110 [keV] (j) E=120 [keV] 図 5.2 10[keV]ごとの散乱シミュレーション (2) （対数プロット） 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/cm^ 2/so urc e] Energy [keV]

23 図 5.3 10[keV]ごとの散乱シミュレーションを重ねたもの （対数プロット） 図 5.3 は、図 5.1 および図 5.1 の結果を同スケールにて重ねたものである。 図 5.3 より、いずれも散乱も 80[keV]近傍に共通して散乱が見られる。 私はこれがバックグラウンドになっていると考え、入射 X 線をそれぞれ取 り除いた後、1 次スペクトルから重み付けをして加算し、シミュレーションに よるバックグラウンドを得ることとした。 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 20 40 60 80 100 120 140 Fl ux [/cm^ 2/so urc e] Energy [keV] E=30[keV] E=40[keV] E=50[keV] E=60[keV] E=70[keV] E=80[keV] E=90[keV] E=100[keV] E=110[keV] E=120[keV]

24

5.2 重み付け加算

X 線管球の 1 次スペクトル（図 1.3）を用い、次の式によってエネルギーご との重みを得た。

重み =

X 線強度 積分強度 (3) 図 5.4 1 次スペクトルより求めた重み 図 5.4 にエネルギーごとの重みを示す。特に、30, 40, 50, 60, 70, 80, 90, 100, 110, 120[keV]の重みを抽出し、重み付け加算を行った。 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 W eig ht Energy [keV]

25

5.3 シミュレーションによるバックグラウンド

シミュレーションによって得られたバックグラウンドを、次に示す。 図 5.5 シミュレーションによって得られたバックグラウンド （線形プロット） 0.0E+00 5.0E-12 1.0E-11 1.5E-11 2.0E-11 2.5E-11 3.0E-11 3.5E-11 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 Fl ux [/cm^ 2/so urc e] Energy [keV]

26

5.4 実測バックグラウンドとの比較

次に、シミュレーションにより得られたバックグラウンドと、実測バックグ ラウンド(図 1.8)の比較を示す。 図 5.6 シミュレーション BG と実測 BG の比較(1) （左縦軸：実測 BG／右縦軸：シミュレーション BG）

27 図 5.7 シミュレーション BG と実測 BG の比較(2) （左縦軸：実測 BG／右縦軸：シミュレーション BG） 図 5.7 は、図 5.6 の横軸レンジを狭めて拡大したものである。図 5.7 よ り、シミュレーションにおいても Pb-𝐾_𝛽線の強度が Pb-𝐾_𝛼線よりも強くなっ ていることが分かる。また、実測バックグラウンドとシミュレーションによ るバックグラウンドでは、Pb-𝐾𝛽線だと考えられるピークの位置が、わずか にずれていることが分かる。 その原因は、Pb-𝐾_𝛽線近傍に何らかの散乱が重なっていることであると推測 される。Pb-𝐾_𝛽線近傍に着目すると、実測バックグラウンドでは１つのピーク のように見えた部分に３つのピークを判別できる。

28

5.5 コンプトン散乱の検討

前項図 5.7 より、Pb-𝐾𝛽線近傍に３本のピークが重なっていることが分かっ た。私は、その原因がコンプトン散乱にあると考えた。 コンプトン散乱とは、入射 X 線が粒子として振る舞い、電子と衝突するこ とで、入射 X 線よりもエネルギーの小さい散乱 X 線を生み出す現象である。 コンプトン散乱は、次の式にて説明できる。ここで、m_eは電子の静止質量、c は光速、𝐸は入射 X 線のエネルギー、𝐸_𝐶は散乱 X 線のエネルギー、𝜑は散乱角 である。

𝐸

_𝐶

=

𝐸 1+ 𝐸 𝑚𝑒c2(1−cos𝜑) (4) 図 5.8 コンプトン散乱の概略図 更に、式(4)を変形することで、次のような式が得られる。

𝜑 = arccos (1 −

𝐸−𝐸𝐶 𝐸×𝐸_𝐶

m𝑐

2

₎

₍₅₎ 式(5)の意義は、入射 X 線および散乱 X 線のエネルギーから散乱角を算出 できることにある。すなわち、特定のエネルギーを有する散乱 X 線の起源を 調査することが可能なのである。

29 前項にて判別した３つのピークを、次のように ROI で分割した。 図 5.9 ３つのピークの ROI 分割 表 5.1 それぞれの ROI 領域の詳細 まず、それぞれの ROI 領域において、どのエネルギーからの寄与が強いの かを検討するため、次のような式を用いて寄与率の計算を行った。その結果を 表 5.2 に示す。

寄与率[%] =

各エネルギーごとの ROI 積分強度 ROI 積分強度

(6)

30 表 5.2 それぞれのピークにおける各エネルギーの寄与率 表 5.2 より、いずれのピークも 90～110[keV]の寄与が強いことが分かる。 よって、散乱 X 線を 90, 100, 110[keV]として、式(5)、および表 5.1 の各ピー クの Peak Energy を用いて、散乱角の計算を行った。 なお、計算に用いた定数は、次の通りである。また、式(5)はエネルギーの 単位がジュール[J]であるため、キロ電子ボルト[keV]からジュール[J]への変 換は、式(7)にて行った。 電子の静止質量 m_e = 9.109 × 10−31 [kg] 真空中の光速 𝑐 = 299792458 [m/s] 電気素量 e = 1.602 × 10−19_[C] 𝐸[J] = 1000 × e × 𝐸[keV] (7)

31 各ピークの散乱角を計算した結果を、次に示す。 表 5.3 各ピークの散乱角 表 5.3 より、散乱角は特定の角度を有していないことが分かる。したがっ て、バックグラウンドには、検出器カバー近傍の空気によるコンプトン散乱 の寄与があると推測される。その概略図を、図 5.8 に示す。 図 5.10 検出器カバー近傍の空気によるコンプトン散乱の概略図

32

5.6 真空雰囲気でのバックグラウンドシミュレーション

真空雰囲気におけるバックグラウンドシミュレーションの結果を、図 5.11 に示す。また、30～120 [keV]それぞれの真空雰囲気における散乱シミュレ ーションの結果を、図 5.12 および 5.13 に示す。 なお、真空雰囲気におけるバックグラウンドシミュレーションは、空気雰 囲気のものと同様に、入射Ｘ線を取り除いたうえで、1 次スペクトルによっ て重み付け加算を行った。 真空雰囲気では、空気雰囲気と比較してスペクトルのばらつきが減少して いる。このことから、バックグラウンドには空気散乱の寄与が大きいことが わかる。 図 5.11 真空雰囲気におけるバッググラウンドシミュレーション

33

(a) E=30 [keV] (b) E=40 [keV]

(c) E=50 [keV] (d) E=60 [keV]

(e) E=70 [keV] (f) E=80 [keV] 図 5.12 10[keV]ごとの散乱シミュレーション 真空雰囲気(1) （対数プロット） 1.000E-11 1.000E-10 1.000E-09 1.000E-08 0.0 50.0 100.0 150.0 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV]

34

(g) E=90 [keV] (h) E=100 [keV]

(i) E=110 [keV] (j) E=120 [keV] 図 5.13 10[keV]ごとの散乱シミュレーション 真空雰囲気(2) （対数プロット） 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV] 1.0E-11 1.0E-10 1.0E-09 1.0E-08 0 50 100 150 Fl ux [/ cm ^ 2 /so u rc e] Energy [keV]

35

6 まとめ

1. PHITS を用いたシミュレーションシステムを構築した。 2. シミュレーションによってバックグラウンドを再現することができた。 3. バックグラウンドの寄与は検出器近傍の空気のコンプトン散乱であると 考えられる。

36

7 今後の課題

1. シミュレーションから𝐼𝐵𝐺(𝐸)を求め、次の式により線源弱係数𝜇(𝐸)を求め ること。

𝜇(𝐸) = −ln⁡(

𝐼(𝐸)+𝐼𝐵𝐺(𝐸) 𝐼₀(𝐸)+𝐼_𝐵𝐺(𝐸)

)

(2) 2. 多重散乱についても検討すること。 3. 検出器カバー内部に実際の寸法に合わせた CdTe を配置し、タリーをそれ に設定すること。更に、実際の実験装置により近づけたジオメトリを作成 すること。

37

8 参考文献

[1] 取越正己, 星和志, 大野由美子, 小林結貴, 阿久津和音, 原澤陽介, 森本一 成, 櫻井浩 (2018). 光子計測 X 線 CT 開発のための基礎研究, 2018 GHMC Physics Division Report, p.55-59

[2] 小林結貴,(2018). 一素子 CdTe 検出器を用いたフォトンカウンティング CT システムの開発, 平成 30 年度群馬大学大学院修士論文 [3] https://www.hamamatsu.com/jp/ja/product/type/L12161-07/index.html [4] https://www.hamamatsu.com/resources/pdf/etd/L12161-07_TXPR1023E.pdf [5] https://www.amptek.com/products/cdte-x-ray-and-gamma-ray-detectors/xr-100cdte-x-ray-and-gamma-ray-detector

[6] Berger, M.J., Hubbell, J.H., Seltzer, S.M., Chang, J., Coursey, J.S., Sukumar, R., Zucker, D.S., and Olsen, K. (2010), XCOM: Photon Cross Section Database (version 1.5). [Online] Available: http://physics.nist.gov/xcom

[2020, February 19]. National Institute of Standards and Technology, Gaithersburg, MD.

[7] Tatsuhiko Sato, Yosuke Iwamoto, Shintaro Hashimoto, Tatsuhiko Ogawa, Takuya Furuta, Shin-ichiro Abe, Takeshi Kai, Pi-En Tsai, Norihiro Matsuda, Hiroshi Iwase, Nobuhiro Shigyo, Lembit Sihver and Koji Niita, Features of Particle and Heavy Ion Transport code System (PHITS) version 3.02, J. Nucl. Sci. Technol. 55, 684-690 (2018)

[8] PHITS 公式ホームページ (online: https://phits.jaea.go.jp/indexj.html) [9] Norman F.M. Henry & Kathleen Lonsdale (Eds.), The International tables

38

9 謝辞

本研究の検討を進める上で大変多くのご指導とご鞭撻を賜り、本論文にお いて終始適切なご指導を頂きました、群馬大学理工学府櫻井浩教授に心より 感謝の意を表し、厚くお礼申し上げます。 本研究において、多くのご指導とご鞭撻を賜りました、群馬大学理工学府花 泉修教授、加田渉准教授に心より感謝の意を表し、厚くお礼申し上げます。 本研究の検討を進める上で大変多くのご指導とご鞭撻を賜り、本論文にお いて終始適切なご指導を頂きました、群馬大学重粒子線医学研究センター取 越正己教授に心より感謝の意を表し、厚くお礼申し上げます。 本研究において、多くのご指導とご鞭撻を賜りました、群馬県立県民健康科 学大学大野由美子准教授に心より感謝の意を表し、厚くお礼申し上げます。 本研究において、多くのご指導とご鞭撻を賜りました、群馬大学理工学府鈴 木宏輔助教授に心より感謝の意を表し、厚くお礼申し上げます。 本研究において、多くのご指導とご鞭撻を賜りました、群馬大学星和志講師 に心より感謝の意を表し、厚くお礼申し上げます。 最後に、日頃より多くのご協力と激励を頂きました群馬大学理工学部櫻井 浩研究室、古澤伸一研究室の皆様に心からお礼申し上げます。 令和 2 年 3 月 4 日 群馬大学大学院 理工学府理工学専攻電子情報・数理教育プログラム 櫻井研究室 修士 2 年 阿久津 和音

39

10 補足

10.1 ジオメトリの作成について

ジオメトリの作成には、CERN（欧州原子核研究機構）開発のモデリングソ フトウェアである SimpleGEO[A]_{を使用した。SimpleGEO は３D モデリング} ソフトウェアであるため、直感的なモデリングが可能である。また、PHITS や FLUKA、MCNP 等の放射線輸送計算コードに対応した形式にモデルを変 換できるという強みがある。

40 図 10.2 SimpleGEO の PHITS 用ジオメトリ変換オプション 具体的なジオメトリ作成の方法については、SimpleGEO 公式サイト[B] の”TUTORIALS”を参考されたい。 なお、PHITS 用に変換したジオメトリには、以下の二点の注意が必要であ る。 １． 所々の文頭にコメント文を意味する『c』が記載されているが、PHITS ver.3.17 ではこれをコメント文と認識できない。ジオメトリ利用前に『c』 を PHITS ver.3.17 がコメント文として認識できる『$』等に書き換える 必要がある。 ２． 変換されたジオメトリには密度情報が書き込まれていない。そのため、 ジオメトリのアスキーファイルに直接物質密度を書き込む必要がある。

41

42

10.2 PHITS 入力データについて

PHITS 入力データおよびジオメトリ入力データの例を、図 10.4 および 10.5 に示す。なお、これは実際に図 3.3 のシミュレーションに使用したものであ る。 図 10.4 使用した PHITS 入力データ(scatteringSimu_120keV.inp)

43 図 10.5 使用したジオメトリ入力データ(geometry.pht) ここでは、特に、PHITS 入力データ（図 10.4）について、その内容の解説 を行う。scatteringSimu_120keV.inp の全文とその内容の解説を次に示す。な お、入射 X 線のエネルギーを指定する箇所とジオメトリ以外は、私が行った すべてのシミュレーションにおいて共通である。詳しくは PHITS ユーザーズ マニュアル[C]_{を参照されたい。}

44 $OMP=4 $ メモリ共有型並列処理の指定（64bit OS のみ対応） $ $OMP=N で CPU コア N 個での並列計算が可能 $ $OMP=0 で使用できる CPU コア全てを使って並列計算を行う [ T i t l e ] $ ここではタイトルを指定する ExprDeviceTest $ 任意タイトルの指定 [ P a r a m e t e r s ] $ ここでは PHITS の動作に関するパラメータを 決定する icntl = 0 $ 計算モードの指定 $ 以下の２つはよく使う $ 0 ：通常の粒子輸送計算 $ 5 ：反応を行わないモード 入射ビームの形状確認に使う $ 8 ：幾何形状表示 ジオメトリの確認に使う maxcas = 10000 $ 1 バッチあたりのヒストリー数 maxbch = 10000 $ バッチ数 $ 大きいほど誤差が小さくなるが $ シミュレーション時間が長くなる $ PHITS では maxcas×macbch [回]シミュレーションが動作する

emin(14) = 1.0e-3 $ 光子輸送エネルギーの下限 [MeV]

dmax(14) = 1000 $ 光子輸送エネルギーの上限 [MeV]

$ 図 10.3 には記述していないが、指定しなければデフォルトでこの値 に定まる

file(1) = c:/phits $ PHITS インストールフォルダ名

file(6) = phits.out $ サマリーの出力ファイル名

[ S o u r c e ] $ ここでは線源粒子の情報を定義する

s-type = 1 $ ビームの形状の指定

45 proj = photon $ 入射粒子の指定 dir = 1.0 $ 入射粒子のｚ軸方向からの方向余弦 $ dir = 1 でｚ軸に沿って真っ直ぐ飛ぶ dom = 21.5 $ 入射粒子方向の立体角範囲 [degree] $ ビームの広がりを指定する $ 立体角なので頂角はその二倍となる $ （dom=21.5 → 頂角 43 度） r0 = 0.0 $ 円柱の半径 [cm] $ r0 = 0 でペンシルビームになる z0 = 34.0 $ 円柱のｚ軸方向の下限 [cm] z1 = 34.0 $ 円柱のｚ軸方向の上限 [cm] e0 = 0.12 $ 入射粒子のエネルギー（単色の場合） [MeV/u] [ M a t e r i a l ] $ ここでは体系を構成する物質を定義する mat[10] Al 1 $ mat[物質番号] 元素 組成比 元素 組成比… といった形式で指定する mat[17] Pb 1 mat[23] W 1 mat[28] N 8 O 2

infl: {geometry.pht} $ infs: {ファイル名}で外部ファイルを挿入する $ ここでは SimpleGEO で作製したジオメトリファイルを挿入している $ ジオメトリファイルの内容は後述するが、[Surface]セクションでセ ルの形状等を定義し、[Cell]セクションでセルの中身を定義している [ T - T r a c k ] $ ここでは粒子の飛程長やフルエンスを導出する タリーである T-track の定義を行う $ タリーセクションには他にも種類があるので、取得したいデータに 応じて変更するべし mesh = xyz $ 領域メッシュの指定 $ 『タリーの指定を xyz 軸で行う』という意味

46 $ 以下、メッシュ定義文は最初の一文字以外共通である x-type = 1 $ x 軸メッシュの指定 $ *-type=1 でデータを群数・分点で与える nx = 1 $ 群数の指定 $ 『ｘ軸を１分割する』という意味 -2.55 2.55 $ 分点の指定 [cm] $ 『タリーのｘ軸範囲を x=-2.25~2.25[cm]に定義 する』という意味 y-type = 1 $ y 軸メッシュの指定 $ *-type=1 でデータを群数・分点で与える ny = 1 $ 群数の指定 $ ここでは『ｙ軸を１分割する』という意味 -1.8 1.8 $ 分点の指定 [cm] $ 『タリーのｙ軸範囲を y=-1.8~1.8[cm]に定義す る』という意味 z-type = 1 $ z 軸メッシュの指定 $ *-type=1 でデータを群数・分点で与える nz = 1 $ 群数の指定 $ ここでは『ｚ軸を１分割する』という意味 98.7 118.5 $ 分点の指定 [cm] $ 『タリーのｚ軸範囲を z=98.7~118.5[cm]に定義 する』という意味 part = all $ タリーで観測する粒子の指定 $ part=all は全ての粒子を観測する $ 細かい指定も可能 e-type = 2 $ エネルギーメッシュの指定 $ *-type=2 でエネルギーを線形で等分割する ne = 500 $ エネルギー分割数 emin = 0.0 $ タリーのエネルギー最小値 [MeV/u] emax = 0.15 $ タリーのエネルギー最大値 [MeV/u]

47 unit = 1 $ 出力単位の指定

$ [T-track]で unit=1 を指定すると、出力単位が $ Flux [1/cm2_{/source]になる}

axis = eng $ 出力データの横軸の指定

$ axis=eng で横軸をエネルギー [MeV/u]に指定するが、axis=xy, yz, zx などに指定すると、粒子の飛程をそれぞれの面から二次元的に見るこ とが可能である

file = track_eng.out $ 出力データのファイル名の指定

title = pipe-Y/sample-Y/colhole-Y/air-vac/cone - E=120 [keV]

$ 出力データのグラフ名の指定 $ 任意名称なのでシミュレーション結果には影響しない gshow = 1 $ グラフの描画に関する指定 epsout = 1 $ グラフを出力するか否かの指定 $ epsout=1 でグラフ出力を行う $ グラフのファイル名は(file で指定した名前).eps になる $ ～（中略）～ $ この部分には、散乱状況を XY、YZ、ZX 平面から眺めた図を出力するよ う設定した[T-track]セクションを３つ設けたが、前述の[T-track]セクショ ンにおいて axis をそれぞれ xy, yz, zx としたものを３つ並べた

[ E n d ]

PHITS ではグラフを画像ファイルのほか、テキストファイルとしても出力 する。このテキストファイルの中には ANGEL がグラフを作る際に使用した 数値データが入っている。これをエクセル等の表計算ソフトで編集すれば、独 自にグラフを作ることも可能である。

48 次に、シミュレーションにおいて使用したジオメトリ（図 3.1）の外部挿入 ファイル（図 10.4）について、その内容の解説を行う。geometry.pht の全文 とその内容の解説を以下に示す。これは SimpleGEO にて作製した 3D モデル を PHITS 用に変換したものである（図 10.2 参照）。 なお、以下は SimpleGEO による変換の後、10.1 節にて述べた点を修正し たものである。詳しくは PHITS ユーザーズマニュアル[C]_{を参照されたい。}

$ PHITS exporter, SimpleGeo proprietary release

$ ********************************************************************** $ Cell Cards $ ********************************************************************** [ C e l l ] $ ここでは[Surface]セクションで定義する面を用いてセルの定義を行う $ なお、記述順は以下の通りである $ [セル番号] [物質番号] [物質密度] [セル定義文] $ 以下にそれぞれの説明を記載する $ [セル番号] [Surface]セクションで設定するセルの番号を記述する $ [物質番号] [Material]セクションで指定した物質番号を記述する 真空の場合は 0 外部ボイドの場合は-1 $ [物質密度] セルに含まれる物質の密度を与える 正の値であれば粒子密度 [1024 _atoms/cm3_] 負の値であれば質量密度 [g/cm3_] セルがボイド（0 あるいは 1）の場合は必要なし $ [セル定義文] セルの幾何形状を記述する 記述には集合代数の演算子を用いる ここでは、+(and)、-(or)、#(not)を用いている $ outerBOX 外箱の設定 00001 17 -11.34 -1 -8 +1 $ tube 遮蔽用パイプの設定 00002 17 -11.34 -1 -10 +11 $ sample 試料の設定 00003 10 -2.70 -1 -9 $ colimator コリメータの設定 00004 23 -19.25 -1 -4 +5

49 $ detector_cover 検出器カバーの設定 00005 17 -11.34 -1 -6 +7 $ void ジオメトリ雰囲気の設定 00006 28 -0.0012 -1 -12 #1 #2 +9 #4 #5 +3 $ outervoid 外部ボイドの設定 PHITS ジオメトリにはこれが不可欠 00007 -1 -2 #6 #1 #2 +9 #4 #5 +3 $ air_in_detector_cover 検出器カバー内部の雰囲気の設定 00008 0 -1 -3 $ ********************************************************************** $ SurfaceCards $ ********************************************************************** [ S u r f a c e ] $ ここでは[Cell]セクションでセルを定義するために使用する面を定義する $ なお、記述順は以下の通りである $ [面番号] [面記号] [面定義数値] $ 以下にそれぞれの説明を記載する $ [面番号] 任意番号(1~999999)を指定する $ [面記号] 面の形状を決定する記号を指定する 本研究で用いた面記号は以下の通りである BOX 任意の BOX 平面（全ての角度が 90 度） SPH 球 RCC 円柱 $ [面定義数値] 定義する面により 1~15 個の入力値を設定する 各面記号に対応する入力値は以下の通りである BOX x0, y0, z0 基点の座標 Ax, Ay, Az 基点からの第１面のベクトル Bx, By, Bz 基点からの第２面のベクトル Cx, Cy, Cz 基点からの第３面のベクトル SPH x0, y0, z0 中心座標 R 球の半径 RCC x0, y0, z0 底面の中心座標 Hx, Hy, Hz 底面の中心からの上面中心への ベクトル R 円柱の半径

50 $ BOX_internal_space 外箱の内部空間の設定 1 BOX -44.90 -43.70 0.50 89.80 0.00 0.00 0.00 69.70 0.00 0.00 0.00 150.10 $ Sphere 外部ボイドの設定 2 SPH 0.00 0.00 60.00 301.00 $ air_in_detector_cover 検出器カバー内部の雰囲気の設定 3 BOX -2.55 -1.80 98.70 5.10 0.00 0.00 0.00 3.60 0.00 0.00 0.00 19.80 $ colimator_disk コリメータの設定 4 RCC 0.00 0.00 98.50 0.00 0.00 0.20 3.50 $ colomator_hole コリメータ穴の設定 5 RCC 0.00 0.00 98.50 0.00 0.00 1.00 0.005 $ detector_cover 検出器カバーの設定 6 BOX -2.75 -2.00 98.70 5.50 0.00 0.00 0.00 4.00 0.00 0.00 0.00 20.00 $ detector_cover_internal_space 検出器カバーの内部空間の設定 7 BOX -2.55 -1.80 98.50 5.10 0.00 0.00 0.00 3.60 0.00 0.00 0.00 20.00 $ outerbox 外箱の設定 8 BOX -45.40 -44.20 0.00 90.80 0.00 0.00 0.00 70.70 0.00 0.00 0.00 150.10 $ sample 試料の設定 9 RCC 0.00 1.50 60.50 -0.00 -3.00 0.00 0.30 $ tube 遮蔽用パイプの設定 10 RCC 0.00 0.00 61.50 0.00 0.00 37.00 1.69 $ tube_internal_space 遮蔽用パイプの内部空間の設定 11 RCC 0.00 0.00 61.50 0.00 0.00 37.00 1.59 $ void ジオメトリ雰囲気の設定 12 SPH 0.00 0.00 60.00 300.00

51

10.3 参考文献（補足）

[A] Theis C., Buchegger K.H., Brugger M., Forkel-Wirth D., Roesler S., Vincke H. ,"Interactive three dimensional visualization and creation of geometries for Monte Carlo calculations", Nuclear Instruments and Methods in Physics Research A 562, pp. 827-829 (2006).

[B] http://theis.web.cern.ch/theis/simplegeo/