Multi-Purpose Particle and Heavy Ion Transport code System PHITS Multi-Purpose Particle and Heavy Ion Transport code System 奨励設定及びツールを使った演習 PHITS講習会　入門実習 2013年4月改訂 title 1

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 2

Recommendation Settings 奨励設定とは？ FAQ マニュアルを読んでも自分が必要とする計算でどのようなパラメータを設定して良いかわからない 典型的な例題に対して，その奨励設定インプットファイルを作成し，アウトプットファイルと共にWeb上で公開 Answer http://phits.jaea.go.jp/examples.html Recommendation Settings 3

Recommendation Settings 奨励設定の内容 DetectorResponse.inp 　　　　検出器の応答関数計算の例題。付与エネルギーのヒストリー間分散を計算 Shielding.inp 加速器遮へい設計の例題。フルエンスより実効線量を導出 ParticleTherapy.inp 粒子線治療計画の例題。水中の線量，線量当量，LET, y分布を計算 PhotonTherapy.inp 光子治療計画の例題。電子加速器から発生するX線や中性子の線量を計算 SemiConductor.inp 半導体ソフトエラー発生率評価の例題。ミクロ領域での吸収線量を計算 NuclearReaction.inp 核反応断面積を計算する例題。２次粒子のフルエンスを計算 H10multiplier.inp [multiplier]セクションを使ってH*(10) を計算 Counter.inp [Counter]を使う例題。1次粒子と2次粒子を識別して吸収線量を計算 Recommendation Settings 4

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 5

ParticleTherapy.inp Tally 水中での吸収線量-深度分布 水中での線量当量-深度分布 Water Phantom 10 cm Carbon 200 MeV/u 10 cm Tally [t-deposit] (file=dose) 水中での吸収線量-深度分布 [t-deposit] (file=dose-equivalent.out) 水中での線量当量-深度分布 [t-LET] 水中でのLET分布（確率密度分布） [t-SED] 水中でのLineal-energy (y) 分布（確率密度分布） ParticleTherapy 6

吸収線量と線量当量の違い ParticleTherapy 7 Dose.eps Dose-equivalent.eps [ T - Deposit ] mesh = r-z r-type = 2 rmin = 0.000000 rmax = 1.000000 nr = 1 z-type = 2 zmin = 0.000000 zmax = 10.00000 nz = 100 dedxfnc = 0 [ T - Deposit ] dedxfnc = 1 R 線量は，usrdfn1.fに書かれた関数により重み付けして出力される Z Q(L) 関数が初期設定で組み込まれている r-z mesh 吸収線量が線量当量に変換され出力 ParticleTherapy 7

LET of C(200 MeV/u) = 16 keV/mm LET & y 分布 [T-LET] [T-SED] LET-distribution.eps (page 1: front surface) y-distribution.eps (page 1: front surface) LET of C(200 MeV/u) = 16 keV/mm yは単色でも確率密度分布を持つ Sharp peak Broad peak ParticleTherapy 8

入射粒子を変更 ParticleTherapy 9 Water Phantom Carbon Proton 200 MeV/u Z Water Phantom X 10 cm Carbon 200 MeV/u Proton (0,0,-20) 10 cm [ S o u r c e ] s-type = 1 proj = 12C e0 = 200.00 r0 = 0.0000 x0 = 0.0000 y0 = 0.0000 z0 = -20.000 z1 = -20.000 dir = 1.0000 [ T - Deposit ] … [ T – L E T ] ... [ T – S E D ] [ T - Deposit ] … [ T - Deposit ] off [ T – L E T ] off ... [ T – S E D ] off Proton Execute Dose.eps 200 MeV 陽子は10cmの水を透過!!! ParticleTherapy 9

体系を変更 General description 10 Water Phantom Carbon Proton 200 MeV/u Y Z Water Phantom X 10 cm 10 cm Carbon 200 MeV/u (0,0,-20) Proton 30 cm 10 cm [ C e l l ] 1 1 -1.0000000E+00 1 -2 -3 2 0 #1 -999 3 -1 999 [ S u r f a c e ] 1 pz 0.0000000E+00 2 pz 1.0000000E+01 3 cz 5.0000000E+00 999 so 1.0000000E+02 Execute 3.0000000E+01 タリー領域の変更忘れ!!! Dose.eps General description 10

タリー領域の変更 ParticleTherapy 11 Water Phantom Carbon Proton 200 MeV/u 統計が Z Water Phantom X 10 cm 10 cm Carbon 200 MeV/u (0,0,-20) Proton 30 cm [ T - Deposit ] title = depth-dose d mesh = r-z x0 = 0.000000 y0 = 0.000000 r-type = 2 rmin = 0.000000 rmax = 1.000000 nr = 1 z-type = 2 zmin = 0.000000 zmax = 10.00000 nz = 100 [ C e l l ] 1 1 -1.0000000E+00 1 -2 -3 2 0 #1 -999 3 -1 999 [ S u r f a c e ] 1 pz 0.0000000E+00 2 pz 2.0000000E+01 3 cz 5.0000000E+00 999 so 1.0000000E+02 set: c1[30.0] 統計が 不十分!! c1 c1 Dose.eps ParticleTherapy 11

接続計算 ParticleTherapy 12 Water Phantom Carbon Proton 200 MeV/u Execute Z Water Phantom X 10 cm 10 cm Carbon 200 MeV/u (0,0,-20) Proton 30 cm [ P a r a m e t e r s ] icntl = 0 maxcas = 100 maxbch = 　2 Execute 10 istdev = -1 200MeV陽子入射に対する線量-深度分布を十分な統計で導出することができた!! Dose.eps ParticleTherapy 12

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 13

H10multiplier.inp 「線量」を計算する３つの方法 電子を除く荷電粒子の電離損失と，光子・中性子のKerma Air 100 cm Concrete 400 MeV Neutron 100 cm 100 cm 「線量」を計算する３つの方法 [t-track] フルエンスに線量換算係数を乗じて導出する方法。この場合は [Multiplier]セクションで定義したH*(10)に対する換算係数を利用 [t-heat] 電子を除く荷電粒子の電離損失と，光子・中性子のKerma 　　　ファクターを用いて領域内の吸収線量を計算　　　 [t-deposit] 電離損失より領域内の吸収線量を計算（中性子・光子の寄与なし） H10multiplier 14

タリー結果 H10multiplier 15 [T-Track] [T-Heat] [T-Deposit] [t-track] フルエンスと線量換算係数を乗じてH*(10)を計算 　　→ コンクリートと空気の間にギャップがない [t-heat] 荷電粒子の電離損失と光子・中性子に対するKerma近似 　　　→コンクリートと空気の間にギャップがある [t-deposit] 荷電粒子による電離損失のみ 　　　 → コンクリートと空気の間にギャップがある。空気線量の統計誤差が大きい H10multiplier 15

H*(10)及び実効線量の計算 Track.eps H10multiplier 16 [t-track]で計算したフルエンスは，内蔵もしくは[multiplier] セクションで定義した複数の関数と掛け合わせることができる [ Multiplier ] number = -250 interpolation = log ne = 140 1.13000E-08 9.14000E+00*3600*1.0e-6 1.42000E-08 9.44000E+00*3600*1.0e-6 1.79000E-08 9.83000E+00*3600*1.0e-6 ... [ T - T r a c k ] title = H*(10) in xyz mesh in uSv/h multiplier = all part = neutron emax = 1000.0 mat mset1 　mset2 all ( 1.0 -250 ) (1.0 -102) part = photon mat mset1 mset2 all ( 1.0 -251) (1.0 -114) H*(10) (Page 1) Effective Dose (Page 2) Track.eps ほとんど同じにしか見えない ! 中性子 (-250) 及び光子 (-251)に対する[multiplier]セクションで定義したH*(10)換算係数 中性子 (-102) 及び光子 (-114)に対する 　　内蔵の実効線量換算係数 H10multiplier 16

Axisの変更 Track.eps H10multiplier 17 Execute H*(10) Effective Dose [ T - T r a c k ] title = H*(10) in xy part = ( neutron photon ) mesh = xyz x-type = 2 xmin = -60.0 xmax = 60.0 nx = 60 y-type = 2 ymin = -60.0 ymax = 60.0 ny = 1 z-type = 2 zmin = 0.0 zmax = 120.0 nz = 60 e-type = 1 ne = 1 1.00000E-10 1.00000E+03 unit = 1 material = all axis = xz file = track.out Execute 1 H*(10) (Page 1) Effective Dose (Page 2) たくさんグラフが出力されてしまうのを防ぐため Track.eps 両グラフの軸が一致していないので どちらが大きいかよく分からない 2Dプロット（等高線図） から 1Dプロット（ヒストグラム）へ変更 z H10multiplier 17

H*(10) < Effective dose 縦軸の調整 [ T - T r a c k ] title = H*(10) in xy part = ( neutron photon ) mesh = xyz x-type = 2 xmin = -60.0 xmax = 60.0 nx = 60 y-type = 2 ymin = -60.0 ymax = 60.0 ny = 1 z-type = 2 zmin = 0.0 zmax = 120.0 nz = 60 e-type = 1 ne = 1 1.00000E-10 1.00000E+03 unit = 1 material = all axis = xz file = track.out Execute 1 H*(10) (Page 1) Effective Dose (Page 2) Track.eps H*(10) < Effective dose z ANGELパラメータを加える （縦軸の最小＆最大値） angel = ymin(5.e-05) ymax(5.e-4) H10multiplier 18

Track.outをTrack2.outとしてコピーして編集 Track2.outを右クリック → 送る → ANGEL #------------------------------------------------ 80行目 #newpage: # no. = 1 ie = 1 ix = 1 iy = 1 … x: z[cm] y: Flux [1/cm^2/source] p: xlin ylog afac(0.8) form(0.9) p: ymin(5.e-05) ymax(5.0e-4) h: n x y(p1-group),hh0l n # z-lower z-upper flux r.err 0.0000E+00 1.0000E+00 1.2809E-04 0.0263 #------------------------------------------------- 232行目 newpage: # no. = 2 ie = 1 ix = 1 iy = 1 0.0000E+00 1.0000E+00 1.9351E-04 0.0174 ... H10 # コメントアウト （PHITS入力 ファイルと同じ形式） Track2.eps レジェンド変更 「(」や「)」は使わない Effective r 色を追加　（r:赤，b：青, g:緑 etc) hh0lとの間にスペースを空けない H10multiplier 19

どのエネルギー領域が支配的？ ≅ < H10multiplier 20 Execute H*(10) Effective Dose [ T - T r a c k ] title = H*(10) in xy part = ( neutron photon ) mesh = xyz x-type = 2 xmin = -60.0 xmax = 60.0 nx = 1 y-type = 2 ymin = -60.0 ymax = 60.0 ny = 1 z-type = 2 zmin = 0.0 zmax = 120.0 nz = 60 e-type = 1 ne = 1 1.00000E-10 1.00000E+03 unit = 1 material = all axis = z file = track.out Low Energy (page1) Low Energy (page3) Execute ≅ High Energy (page2) High Energy (page4) < 2 20.0 1.0E+03 エネルギービンを２つに分割 angel = ymin(1.e-05) ymax(2.e-4) H*(10) Effective Dose H10multiplier 20

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 21

便利なツールの内容 PHITSでシミュレーションした粒子飛跡の動画作成方法 [t-3dshow]で表示したジオメトリを回転させる動画作成方法 Animation PHITSでシミュレーションした粒子飛跡の動画作成方法 Rotate3dshow [t-3dshow]で表示したジオメトリを回転させる動画作成方法 SimpleGEO SimpleGEOで作ったジオメトリをPHITSに組み込む方法 Autorun 　　少しずつ計算条件を変えてPHITSを連続的に実行する方法 Animation Rotate3dshow SimpleGEO Utilities 22

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 23

粒子飛跡のアニメーション化 必要ソフトウェア ImageMagick (http://www.imagemagick.org/script/binary-releases.php) 　　複数ページのEPSファイルを１つのGIFアニメーションにするためのソフト 手順 1. PHITSを実行する 　　[t-track]や[t-deposit]タリーで時間メッシュ（t-type）を導入し， 　　フラックスや発熱量の時間変化を出力する 2. EPSファイルを編集する 　　① 出力したEPSファイル（track.eps）をテキストエディタで開く 　　② 「PageBoundingBox」を２回検索する 　　③ １つ目と２つ目の数字を0に変更して保存する 3. EPSファイルをGIFアニメーションに変換する 　　① コマンドプロンプトを開いて，EPSファイルのあるフォルダに移動する 　　② ‘convert -dispose background -rotate 90 XXX.eps XXX.gif‘ %%PageBoundingBox: 27 36 571 803 %%PageBoundingBox: 0 0 571 803 Animation 24

動画時間分解能の向上 Animation 25 animation.inp 20 frame PHITSを実行 EPSを編集 [ T - T r a c k ] C -- Contour figure Tally -- mesh = xyz ... t-type = 2 nt = 20 # Number of frame tmin = 0.00 # Initial time (nsec) tmax = 40 # Final time (nsec) angel = cmin(1.e-05) cmax(1.e+00) epsout = 1 60 20 frame フレーム数を増やす PHITSを実行 EPSを編集 EPSからGIFに変換 60 frame Animation 25

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 26

SimpleGEOを用いた体系作成方法 SimpleGEO＆プラグインパッケージ（Windows専用） SimpleGEO 27 必要ソフトウェア SimpleGEO＆プラグインパッケージ（Windows専用） Download site: http://theis.web.cern.ch/theis/simplegeo/ 手順 1. SimpleGEOを用いて体系を作成する 必ずVoid（空気）の部分も定義すること。詳細はSimpleGEOのマニュアル参照 2. SimpleGEOの体系をPHITS形式で出力（File → Export → PHITS) [cell] と [surface] セクションのみ出力される 3. [cell]と[surface]セクションを除いたPHITS入力ファイルを作成する ‘infl:’ コマンドを使って，手順２で作成したExportファイルをインクルードする 4. 作成した入力ファイルを使ってPHITSを実行する 5. PHITSのmesh=xyzを使ったタリーの結果をSimpleGEOで可視化する SimpleGEOでプラグインDaVis3Dを読み込む（Macros → Load Plugins → DaVis3D) xyz-meshタリー出力を選択して，可視化する　　 SimpleGEO 27

SimpleGEOでの体系変更 SimpleGEO 28 doll.pht Export to PHITS [ C e l l ] c Body 00001 26 -1 -1 c Eyes 00002 11 -7.87 -6 : -3 c Head 00003 6 -1.9 -2 #2 c Void 00004 0 -4 +1 #2 #3 +7 +8 c Outervoid 00005 -1 -5 +4 c leg1 00006 26 -1 -7 c leg2 00007 26 -1 -8 [ C e l l ] c Body 00001 26 -1 -1 c Eyes 00002 11 -7.87 -6 : -3 c Head 00003 6 -1.9 -2 #2 c Void 00004 0 -4 +1 #2 #3 c Outervoid 00005 -1 -5 +4 [ S u r f a c e ] c Body 1 RCC 0.00 0.00 0.00 0.00 0.00 10.00 5.00 c Headball 2 SPH 0.00 0.00 15.00 5.00 c LeftEye 3 SPH -3.00 4.00 16.00 1.00 c OuterSphere 4 SPH 0.00 0.00 0.00 100.00 c Outmostsphere 5 SPH 0.00 0.00 0.00 100.00 c RightEye 6 SPH 3.00 4.00 16.00 1.00 [ S u r f a c e ] c Body 1 RCC 0.00 0.00 0.00 0.00 0.00 10.00 5.00 c Headball 2 SPH 0.00 0.00 15.00 5.00 c LeftEye 3 SPH -3.00 4.00 16.00 1.00 c OuterSphere 4 SPH 0.00 0.00 0.00 100.00 c Outmostsphere 5 SPH 0.00 0.00 0.00 100.00 c RightEye 6 SPH 3.00 4.00 16.00 1.00 c leg1 7 RCC -2.50 0.00 -10.00 0.00 0.00 10.00 2.00 c leg2 8 RCC 2.50 0.00 -10.00 0.00 0.00 10.00 2.00 Export to PHITS SimpleGEO 28

PHITS出力の３次元可視化 SimpleGEO 29 SimpleGEO.inp Execute SimpleGEO [ T - Deposit ] title = [t-deposit] in xyz mesh mesh = xyz # mesh type is xyz scoring mesh x-type = 2 # x-mesh is linear given by xmin, xmax and nx xmin = -5.000000 # minimum value of x-mesh points xmax = 5.000000 # maximum value of x-mesh points nx = 20 # number of x-mesh points y-type = 2 # y-mesh is linear given by ymin, ymax and ny ymin = -5.000000 # minimum value of y-mesh points ymax = 5.000000 # maximum value of y-mesh points ny = 20 # number of y-mesh points z-type = 2 # z-mesh is linear given by zmin, zmax and nz zmin = 0.000000 # minimum value of z-mesh points zmax = 20.00000 # maximum value of z-mesh points nz = 40 # number of z-mesh points -10.00000 Execute タリー領域の拡張 SimpleGEO ‘Macros -> Load Plugins -> DaVis3D’を選択 ‘DaVis3D‘ボタンを押す ‘deposity-xy.out’を選択し，‘Load data’を押す SimpleGEO 29

実習内容 奨励設定 便利なツール まとめと演習問題 ParticleTherapy H10mupliplier 粒子飛跡のアニメーション化 SimpleGEO 「実習」　前の基本的な話から Contents 30

１GeV陽子をPHITS形状の水ファントムに入射したときの吸収線量分布 まとめ PHITSの入力ファイルを全て自分で作るのは大変 なので，奨励設定から近い例題を選んで，それを 編集して入力ファイルを作成する 追加してほしい例題・ツールなどがありましたら， PHITS事務局までご連絡下さい １GeV陽子をPHITS形状の水ファントムに入射したときの吸収線量分布 http://phits.jaea.go.jp 31

演習問題 250MeV/uの11Bビームを水ファントムに入射したときの吸収線量の深さ分 布を，様々な半径方向に対して計算する 計算した吸収線量の深さ分布を，寄与する粒子別（中性子，陽子，He, Li, Be, B）に出力する ファントムを２つの領域に分割し，片方をアルミ製 (2.7 g/cm3) にする 試行回数を増やす，もしくは接続計算して，統計誤差を小さくする ヒント 奨励設定「ParticleTherapy」にある１つ目の[t-deposit]タリーを使う [t-deposit]タリーの半径方向に対するメッシュ数(nr)を増やす 計算時間を短くするため，不要なタリーを「off」コマンドを用いてコメントアウト [source]セクションにて線源を変更 [t-deposit]タリーの粒子(part)を変更する 新しい物質，サーフェス，セルを定義して，アルミファントムを加える Summary 32

演習問題回答の例 考えてみよう Summary 33 半径1cm以内及び1 cmから2 cmの範囲内における吸収線量の深さ分布 r < 1 cm 1 cm < r < 2 cm 半径1cm以内及び1 cmから2 cmの範囲内における吸収線量の深さ分布 考えてみよう ブラッグピークより後方では，どのような粒子が線量に寄与しているか？ なぜ深さ5cmのところで，吸収線量が急激に大きくなるのか？ なぜこの計算で中性子の寄与はまったくないのか? Summary 33