PHITS

PHITS
Multi-Purpose Particle and Heavy Ion Transport code System
総合実習（I）:
奨励設定及びツールを使った演習
2013年4月改訂
title
1
実習内容
• 奨励設定
• ParticleTherapy
• H10mupliplier
• 便利なツール
• 粒子飛跡のアニメーション化
• SimpleGEO
• まとめと演習問題
Contents
2
奨励設定とは？
FAQ
マニュアルを読んでも自分が必要とする計算でどのよう
なパラメータを設定して良いかわからない
Answer
典型的な例題に対して，その奨励設定インプットファイル
を作成し，アウトプットファイルと共にWeb上で公開
http://phits.jaea.go.jp/examples.html
Recommendation Settings
3
奨励設定の内容
 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
Carbon 200 MeV/u
Water
Phantom
10 cm
ParticleTherapy.inp
10 cm
Tally
 [t-deposit] (file=dose)
水中での吸収線量-深度分布
 [t-deposit] (file=dose-equivalent.out)
水中での線量当量-深度分布
 [t-LET]
水中でのLET分布（確率密度分布）
 [t-SED]
水中でのLineal-energy (y) 分布（確率密度分布）
ParticleTherapy
6
吸収線量と線量当量の違い
Dose.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
Dose-equivalent.eps
[ T - Deposit ]
dedxfnc = 1
R
線量は，usrdfn1.fに書かれた関数
により重み付けして出力される
Z
r-z mesh
Q(L) 関数が初期設定で組み込まれ
ている
吸収線量が線量当量に変換され出力
ParticleTherapy
7
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
入射粒子を変更
Z
X
(0,0,-20)
Water
Phantom
Carbon
Proton
200 MeV/u
10 cm
Y
10 cm
[Source]
s-type = 1
proj = 12C
Proton
e0 = 200.00
r0 = 0.0000
x0 = 0.0000
y0 = 0.0000
z0 = -20.000
z1 = -20.000
dir = 1.0000
[ T - Deposit ]
…
[ T - Deposit ] off
…
Execute
[ T – L E T ] off
...
[ T – S E D ] off
...
200 MeV 陽子は10cmの水を透過!!!
ParticleTherapy
Dose.eps
9
体系を変更
X
(0,0,-20)
Water
Phantom
Carbon
Proton
200 MeV/u
10 cm
10 cm
Z
10 cm
Y
30 cm
[Cell]
1
1 -1.0000000E+00 1 -2 -3
2
0
#1 -999
3
-1
999
[Surface]
1
pz
0.0000000E+00
2
pz
1.0000000E+01
3.0000000E+01
3
cz
5.0000000E+00
999
so
1.0000000E+02
Execute
タリー領域の変更忘れ!!!
General description
Dose.eps
10
タリー領域の変更
X
(0,0,-20)
Carbon
Proton
200 MeV/u
Water
Phantom
10 cm
Z
10 cm
Y
30 cm
[Cell]
1
1 -1.0000000E+00 1 -2 -3
2
0
#1 -999
3
-1
999
set: c1[30.0]
[Surface]
1
pz
0.0000000E+00
c1
2
pz
2.0000000E+01
3
cz
5.0000000E+00
999
so
1.0000000E+02
[ 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
c1
zmax = 10.00000
nz = 100
ParticleTherapy
統計が
不十分!!
Dose.eps
11
接続計算
X
(0,0,-20)
Carbon
Proton
200 MeV/u
Water
Phantom
10 cm
Z
10 cm
Y
30 cm
[Parameters]
icntl =
0
maxcas =
100
maxbch =
102
istdev = -1
Execute
200MeV陽子入射に対する線量-深度分布
を十分な統計で導出することができた!!
ParticleTherapy
Dose.eps
12
実習内容
• 奨励設定
• ParticleTherapy
• H10mupliplier
• 便利なツール
• 粒子飛跡のアニメーション化
• SimpleGEO
• まとめと演習問題
Contents
13
Air
400 MeV Neutron
Concrete
100 cm
H10multiplier.inp
100 cm
「線量」を計算する３つの方法
 [t-track]
フルエンスに線量換算係数を乗じて導出する方法。この場合は
[Multiplier]セクションで定義したH*(10)に対する換算係数を利用
 [t-heat]
電子を除く荷電粒子の電離損失と，光子・中性子のKerma
ファクターを用いて領域内の吸収線量を計算
 [t-deposit]
電離損失より領域内の吸収線量を計算（中性子・光子の寄与なし）
H10multiplier
14
タリー結果
[T-Track]
[T-Heat]
[T-Deposit]
 [t-track] フルエンスと線量換算係数を乗じてH*(10)を計算

*
H
(10)

→ コンクリートと空気の間にギャップがない
0  (E)d (E)dE
 [t-heat] 荷電粒子の電離損失と光子・中性子に対するKerma近似
→コンクリートと空気の間にギャップがある
 [t-deposit] 荷電粒子による電離損失のみ
→ コンクリートと空気の間にギャップがある。空気線量の統計誤差が大きい
H10multiplier
15
H*(10)及び実効線量の計算
[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-Track]
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)
multiplier = all
part = photon
emax = 1000.0
mat
mset1 mset2
all ( 1.0 -251) (1.0 -114)
H*(10)
Effective Dose
(Page 1)
(Page 2)
Track.eps
ほとんど同じにしか見えない !
•
中性子 (-250) 及び光子 (-251)に対する
[multiplier]セクションで定義したH*(10)換算係数
• 中性子 (-102) 及び光子 (-114)に対する
内蔵の実効線量換算係数
H10multiplier
16
Axisの変更
[T-Track]
title = H*(10) in xy
part = ( neutron photon )
mesh = xyz
x-type = 2
xmin = -60.0
xmax = 60.0 Execute
nx = 60 1
y-type = 2
ymin = -60.0
ymax = 60.0
ny = 1
z-type = 2
H*(10)
Effective Dose
たくさんグラフが
zmin = 0.0
(Page 1)
(Page 2)
出力されてしまう
zmax = 120.0
のを防ぐため
nz = 60
e-type = 1
両グラフの軸が一致していないので
ne = 1
1.00000E-10 1.00000E+03
どちらが大きいかよく分からない
unit = 1
material = all
2Dプロット（等高線図）から
axis = xzz
1Dプロット（ヒストグラム）へ変更
file = track.out
Track.eps
H10multiplier
17
縦軸の調整
[T-Track]
title = H*(10) in xy
part = ( neutron photon )
mesh = xyz
x-type = 2
xmin = -60.0
xmax = 60.0 Execute
nx = 601
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
z
axis = xz
file = track.out
angel = ymin(5.e-05) ymax(5.e-4)
H*(10)
Effective Dose
(Page 1)
(Page 2)
Track.eps
H*(10) < Effective dose
ANGELパラメータを加える
（縦軸の最小＆最大値）
H10multiplier
18
グラフを加工する（ANGELの実行）
Track.outをTrack2.outとしてコピーして編集
#------------------------------------------------ 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
H10
# 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
（PHITS入力
…
ファイルと同じ形式）
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
rn
Effective
# z-lower
z-upper
flux
r.err
0.0000E+00 1.0000E+00 1.9351E-04 0.0174
...
Track2.outを右クリック → 送る →
ANGEL
Track2.eps
レジェンド変更
「(」や「)」は使わない
色を追加（r:赤，b：青, g:緑 etc)
hh0lとの間にスペースを空けない
H10multiplier
19
どのエネルギー領域が支配的？
[T-Track]
title = H*(10) in xy
part = ( neutron photon )
mesh = xyz
x-type = 2
xmin = -60.0
xmax = 60.0 Execute
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 = 12
1.00000E-10 1.00000E+03
20.0 1.0E+03
unit = 1
material = all エネルギービン
を２つに分割
axis =
z
file = track.out
angel = ymin(1.e-05) ymax(2.e-4)
Low Energy
(page1)
Low Energy
(page3)
≅
High Energy
(page2)
High Energy
(page4)
<
H*(10)
H10multiplier
Effective Dose
20
実習内容
• 奨励設定
• ParticleTherapy
• H10mupliplier
• 便利なツール
• 粒子飛跡のアニメーション化
• SimpleGEO
• まとめと演習問題
Contents
21
便利なツールの内容
 Animation
PHITSでシミュレーションした粒子飛跡の動画作成方法
 Rotate3dshow
[t-3dshow]で表示したジオメトリを回転させる動画作成方法
 SimpleGEO
SimpleGEOで作ったジオメトリをPHITSに組み込む方法
 Autorun
少しずつ計算条件を変えてPHITSを連続的に実行する方法
Animation
Rotate3dshow
Utilities
SimpleGEO
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に変更して保存する
%%PageBoundingBox:
27
36
571
803
%%PageBoundingBox:
0
0
571
803
3. EPSファイルをGIFアニメーションに変換する
① コマンドプロンプトを開いて，EPSファイルのあるフォルダに移動する
② ‘convert -dispose background -rotate 90 XXX.eps XXX.gif‘
Animation
24
動画時間分解能の向上
animation.inp
[T-Track]
C
-- Contour figure Tally -mesh = xyz
...
t-type = 2
nt = 2060
# 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
20 frame
フレーム数を増やす
PHITSを実行
EPSを編集
EPSからGIFに変換
60 frame
Animation
25
実習内容
• 奨励設定
• ParticleTherapy
• H10mupliplier
• 便利なツール
• 粒子飛跡のアニメーション化
• SimpleGEO
• まとめと演習問題
Contents
26
SimpleGEOを用いた体系作成方法
 必要ソフトウェア
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での体系変更
doll.pht
[Cell]
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
[Surface]
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.inp
[ 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 = -10.00000
0.000000
# minimum value of z-mesh points
zmax = 20.00000
# maximum value of z-mesh points
nz = 40
# number of z-mesh points
Execute
タリー領域の拡張
SimpleGEO
1. ‘Macros -> Load Plugins -> DaVis3D’を選択
2. ‘DaVis3D‘ボタンを押す
3. ‘deposity-xy.out’を選択し，‘Load data’を押す
SimpleGEO
29
実習内容
• 奨励設定
• ParticleTherapy
• H10mupliplier
• 便利なツール
• 粒子飛跡のアニメーション化
• SimpleGEO
• まとめと演習問題
Contents
30
まとめ
• PHITSの入力ファイルを全て自分で作るのは大変
なので，奨励設定から近い例題を選んで，それを
編集して入力ファイルを作成する
• 追加してほしい例題・ツールなどがありましたら，
PHITS事務局までご連絡下さい
１GeV陽子をPHITS形状の水ファントムに入射したときの吸収線量分布
http://phits.jaea.go.jp
31
演習問題
1. 250MeV/uの11Bビームを水ファントムに入射したときの吸収線量の深さ分
布を，様々な半径方向に対して計算する
2. 計算した吸収線量の深さ分布を，寄与する粒子別（中性子，陽子，He, Li,
Be, B）に出力する
3. ファントムを２つの領域に分割し，片方をアルミ製 (2.7 g/cm3) にする
4. 試行回数を増やす，もしくは接続計算して，統計誤差を小さくする
ヒント
•
•
•
•
•
•
奨励設定「ParticleTherapy」にある１つ目の[t-deposit]タリーを使う
[t-deposit]タリーの半径方向に対するメッシュ数(nr)を増やす
計算時間を短くするため，不要なタリーを「off」コマンドを用いてコメントアウト
[source]セクションにて線源を変更
[t-deposit]タリーの粒子(part)を変更する
新しい物質，サーフェス，セルを定義して，アルミファントムを加える
Summary
32
演習問題回答の例
r < 1 cm
1 cm < r < 2 cm
半径1cm以内及び1 cmから2 cmの範囲内における吸収線量の深さ分布
考えてみよう
• ブラッグピークより後方では，どのような粒子が線量に寄与しているか？
• なぜ深さ5cmのところで，吸収線量が急激に大きくなるのか？
• なぜこの計算で中性子の寄与はまったくないのか?
Summary
33

Download Report