Multi-Purpose Particle and Heavy Ion Transport code System

Presentation on theme: "Multi-Purpose Particle and Heavy Ion Transport code System"— Presentation transcript:

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

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

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

4 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

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

6 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

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

8 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

9 入射粒子を変更 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 = r0 = x0 = y0 = z0 = z1 = dir = [ 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

10 体系を変更 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 ] E #1 -999 [ S u r f a c e ] pz E+00 pz E+01 cz E+00 so E+02 Execute E+01 タリー領域の変更忘れ!!! Dose.eps General description 10

11 タリー領域の変更 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 = y0 = r-type = 2 rmin = rmax = nr = 1 z-type = 2 zmin = zmax = nz = 100 [ C e l l ] E #1 -999 [ S u r f a c e ] pz E+00 pz E+01 cz E+00 so E+02 set: c1[30.0] 統計が 不十分!! c1 c1 Dose.eps ParticleTherapy 11

12 接続計算 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 = maxcas = maxbch = 　2 Execute 10 istdev = -1 200MeV陽子入射に対する線量-深度分布を十分な統計で導出することができた!! Dose.eps ParticleTherapy 12

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

14 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

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

16 H*(10)及び実効線量の計算 Track.eps H10multiplier 16
[t-track]で計算したフルエンスは，内蔵もしくは[multiplier] セクションで定義した複数の関数と掛け合わせることができる [ Multiplier ] number = -250 interpolation = log ne = 140 E E+00*3600*1.0e-6 E E+00*3600*1.0e-6 E E+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 = mat mset1 　mset2 all ( ) ( ) part = photon mat mset1 mset2 all ( ) ( ) H*(10) (Page 1) Effective Dose (Page 2) Track.eps ほとんど同じにしか見えない ! 中性子 (-250) 及び光子 (-251)に対する[multiplier]セクションで定義したH*(10)換算係数 中性子 (-102) 及び光子 (-114)に対する 　　内蔵の実効線量換算係数 H10multiplier 16

17 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 = xmax = nx = 60 y-type = 2 ymin = ymax = ny = 1 z-type = 2 zmin = zmax = nz = 60 e-type = 1 ne = 1 E E+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

18 H*(10) < Effective dose
縦軸の調整 [ T - T r a c k ] title = H*(10) in xy part = ( neutron photon ) mesh = xyz x-type = 2 xmin = xmax = nx = 60 y-type = 2 ymin = ymax = ny = 1 z-type = 2 zmin = zmax = nz = 60 e-type = 1 ne = 1 E E+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

19 Track.outをTrack2.outとしてコピーして編集 Track2.outを右クリック → 送る → ANGEL
# 行目 #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 E E # 行目 newpage: # no. = 2 ie = 1 ix = 1 iy = 1 0.0000E E E ... H10 # コメントアウト （PHITS入力 ファイルと同じ形式） Track2.eps レジェンド変更 「(」や「)」は使わない Effective r 色を追加　（r:赤，b：青, g:緑 etc) hh0lとの間にスペースを空けない H10multiplier 19

20 どのエネルギー領域が支配的？ ≅ < 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 = xmax = nx = 1 y-type = 2 ymin = ymax = ny = 1 z-type = 2 zmin = zmax = nz = 60 e-type = 1 ne = 1 E E+03 unit = 1 material = all axis = z file = track.out Low Energy (page1) Low Energy (page3) Execute ≅ High Energy (page2) High Energy (page4) < 2 E+03 エネルギービンを２つに分割 angel = ymin(1.e-05) ymax(2.e-4) H*(10) Effective Dose H10multiplier 20

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

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

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

24 粒子飛跡のアニメーション化 必要ソフトウェア ImageMagick ( 　　複数ページの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: %%PageBoundingBox: Animation 24

25 動画時間分解能の向上 Animation 25 animation.inp 20 frame PHITSを実行 EPSを編集
[ T - T r a c k ] C Contour figure Tally -- mesh = xyz ... t-type = 2 nt = # Number of frame tmin = # Initial time (nsec) tmax = # 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

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

27 SimpleGEOを用いた体系作成方法 SimpleGEO＆プラグインパッケージ（Windows専用） SimpleGEO 27
必要ソフトウェア SimpleGEO＆プラグインパッケージ（Windows専用） Download site: 手順 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

28 SimpleGEOでの体系変更 SimpleGEO 28 doll.pht Export to PHITS [ C e l l ]
c Body c Eyes : -3 c Head #2 c Void #2 # c Outervoid c leg1 c leg2 [ C e l l ] c Body c Eyes : -3 c Head #2 c Void #2 #3 c Outervoid [ S u r f a c e ] c Body 1 RCC c Headball 2 SPH c LeftEye 3 SPH c OuterSphere 4 SPH c Outmostsphere 5 SPH c RightEye 6 SPH [ S u r f a c e ] c Body 1 RCC c Headball 2 SPH c LeftEye 3 SPH c OuterSphere 4 SPH c Outmostsphere 5 SPH c RightEye 6 SPH c leg1 7 RCC c leg2 8 RCC Export to PHITS SimpleGEO 28

29 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 = # x-mesh is linear given by xmin, xmax and nx xmin = # minimum value of x-mesh points xmax = # maximum value of x-mesh points nx = # number of x-mesh points y-type = # y-mesh is linear given by ymin, ymax and ny ymin = # minimum value of y-mesh points ymax = # maximum value of y-mesh points ny = # number of y-mesh points z-type = # z-mesh is linear given by zmin, zmax and nz zmin = # minimum value of z-mesh points zmax = # maximum value of z-mesh points nz = # number of z-mesh points Execute タリー領域の拡張 SimpleGEO ‘Macros -> Load Plugins -> DaVis3D’を選択 ‘DaVis3D‘ボタンを押す ‘deposity-xy.out’を選択し，‘Load data’を押す SimpleGEO 29

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

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

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

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