PowerPoint プレゼンテーション

PHITS
Multi-Purpose Particle and Heavy Ion Transport code System
PHITS 講習会 基礎実習(III):
計算条件の設定
2015年9月改訂
title
1
はじめに
• PHITS計算は,[Parameters]セクションで定義され
る様々なパラメータ(物理モデルの選択など)により
コントロールされています。
• 全てのパラメータには初期設定値があり,基本的に
は変更する必要はありません。
• ただし,最適な計算結果を得るためには,いくつか
のパラメータを変更する必要があります。
本実習では,パラメータの設定方法について学習します
2
本実習の目標
PHITSの便利な機能を活用し,統計数や核反応モデル
を変化させ,最適な計算結果が得られるようになる
宿題体系における初期設定での
陽子(上)・中性子(下)フルエンス
ヒストリー数,切断エネルギー,核データを適切
に設定した陽子(上)・中性子(下)フルエンス
3
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
4
計算モードの選択
通常輸送計算
反応なしモード
体系確認モード
5
計算モードの選択
体系確認モード
まずは、基礎実習(II)で
行ったジオメトリを確認
するタリーを用いて、サ
ンプルインプットの体系
を確認してみましょう。
6
体系の確認([t-3dshow])
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 11
[t-3dshow]を実行
maxcas = 100
maxbch = 10
file(6) = phits.out
[T-3Dshow]
output = 3
material = -1
6
x0 = 0
y0 = 0
z0 = 0
e-the = 70 $ eye
e-phi = 20
e-dst = 80
l-the = 20 $ light
l-phi = 0
l-dst = 100
w-wdt = 50 $ window
w-hgt = 50
w-dst = 25
heaven = z
line = 1
shadow = 2
1層5cm幅の玉ねぎ構造
file = 3dshow.out
title = Check onion structure using [T-3dshow] tally
epsout = 1
・・・・・・
3dshow.eps
7
体系の確認(gshowオプション)
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 11
8
修正して実行
maxcas = 100
maxbch = 10
file(6) = phits.out
[T-Track]
mesh = xyz 各領域のセル番号と満たさ
x-type = 2 れている物質名をチェック!
nx = 100
xmin = -50.
xmax = 50.
y-type = 1
ny = 1
-5.0 5.0
z-type = 2
nz = 100
zmin = -50.
zmax = 50.
・・・・・・
axis = xz
file = track_xz.out
title = Check source direction using [T-track] tally
gshow = 3
epsout = 1
track_xz.eps
8
線源の確認([t-track])
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 5
maxcas = 100
maxbch = 10
file(6) = phits.out
修正して実行
[T-Track]
mesh = xyz
x-type = 2
nx = 100
xmin = -50.
xmax = 50.
y-type = 1
ny = 1
-5.0 5.0
z-type = 2
nz = 100
zmin = -50.
zmax = 50.
e-type = 1
ne = 1
0. 200.
unit = 1
axis = xz
file = track_xz.out
title = Check source direction using [T-track] tally
epsout = 1
9
線源の飛跡確認
track-xz.eps
全ての物質がvoidとなるため,まっすぐに飛んでいく
(線源の発生位置・方向が確認できる)
10
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
11
外部ファイルの挿入
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 5
maxcas = 100
maxbch = 10
file(6) = phits.out
set: c1[10]
[Source]
・・・・・・
•
onion.inp
[Material]
・・・・・・
[Mat Name Color]
・・・・・・
[Surface]
10 so
500.
11 so
5.
12 so
10.
13 so
15.
14 so
20.
15 so
25.
[Cell]
100 -1
10
infl: {onion.inp}[1-33]
101
1 -19.32
-11
102
2 -1.
11 -12
103
3 -8.93
12 -13
インプットファイルとは別のファイルを取り込む場合は”infl:”コマンドを使う
104
4 -1.
13 -14
{}でファイル名,[n1-n2]で取り込む行の範囲を指定する
105
5 -0.9
14 -15
106
6 -1.20e-3
15 -10
インプットファイルの一行目に“file = インプットファイル名”を書く(重要)
•
• “off”したセクションの下にある場合、無視されるので注意
12
ユーザー定義定数と数式の利用
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
・・・・・・
set: c1[10]
[Source]
s-type = 1
proj = proton
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 150
totfact = pi*c1**2
ユーザー定義定数の利用:
同じ値を使う場合は,予め定義しておくと便利!
PHITSでは定数を定義して、インプットファイル内でその定数を参照
することが可能。
“set: ci[x]”
•名前はc1からc99まで使用可能
•書いた箇所より下の部分で有効となる。また、何度も同じ名前を
定義できる
•“off”で飛ばしたセクションの下にある場合、無視されるので注意
数式の利用:
PHITS入力ファイルでは,FORTRAN形式の数式が
利用可能
• 円周率(pi)は定義しなくても使える
• べき乗は「**」なので注意
• expやcosなどの関数も使える
13
規格化について
lec03.inp
file = lec03.inp
[Title]
・・・・・・
[Parameters]
・・・・・・
set: c1[10]
[Source]
s-type = 1
proj = proton
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 150
totfact = pi*c1**2
• PHITSのタリー結果は,発生する線源数*(試行
回数)当たりで出力される
(*厳密には,線源粒子のウェイト当たり)
• 測定値などと比較するには,Bqや/cm2/sで表され
る線源の発生頻度で規格化する必要がある
• タリー結果を定数倍(totfactで定義)することによ
り,直接, /cm2/sやmGy/hなどの単位で出力でき
るようになる
• 今回の場合,規格化しないと線源のフルエンスは
1/πr2 (/cm2/source)となるため,タリー結果にπr2を
乗じることでフルエンスが1(/cm2)となるように規
格化している
s-type=1でz0=z1とすると、半径 r0
r0の円の面が線源領域となる
14
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
15
体積、面積計算の原理
モンテカルロ積分: 乱数を用いて定積分の近似解を求める
方法。積分範囲の内か外かを判定することにより、積分値を
確率的に評価する。計算精度は試行回数に依存する。
体積計算にこの考え方を適用すると…
線の長さの合計 [cm]
= 線の数(試行回数) [source]
× 平均track length [cm/source]
求める体積 [cm3]
= 線の長さの合計 ÷ 線の面密度
= 平均track length ÷ 平均flux
= [t-track]の結果 × 線源面積
試行回数を増やすこ
とで精度が向上する
線の面密度(間隔) [1/cm2]
= 線の数(試行回数) [source]
× 平均flux [1/cm2/source]
あらかじめ
予測可能
(今回はπr2 )
16
課題1: 体系の体積を計算してみよう
lec03.inp
[Parameters]
icntl = 5
maxcas = 100
maxbch = 10
file(6) = phits.out
set: c1[30]
[Source]
s-type = 1
proj = proton
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 150
totfact = pi*c1**2
[T-Track]
mesh = reg
reg = 101 102 103 104 105
・・・・・・
・・・・・・
file = volume.out
修正して実行
体系の体積や面積を計算する場合:
• icntl=5,
• r0を体系が全部入る大きさにする
17
体積計算の結果
モンテカルロ積分の要領で体系の体積を求めることが可能
volume.out
・・・・
#num
1
2
3
4
5
・・
reg
101
102
103
104
105
volume
1.0000E+00
1.0000E+00
1.0000E+00
1.0000E+00
1.0000E+00
flux
3.6914E+02
3.5499E+03
1.0295E+04
1.9822E+04
3.1834E+04
r.err
0.2362
0.0965
0.0577
0.0380
0.0247
各層の体積は?
•4π(5)3/3=524 cm3
•4π(10)3/3 - 524=3665 cm3
•4π(15)3/3 - 4189=9948 cm3
•4π(20)3/3 - 14137=19373 cm3
•4π(25)3/3 - 33510=31940 cm3
内側の結果ほど一致していない
→ 統計量が十分ではない
→ 統計量を増やすには?
18
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
19
試行回数(ヒストリ数)の設定
モンテカルロ計算の統計精度は、シミュレーションの試行回
数に依存しています。したがって、信頼できる計算結果を得る
ためには、十分な数の試行回数を設定する必要があります。
maxcas
(D=10)
1バッチのヒストリー数
maxbch
(D=10)
バッチ数
rseed
(D=0.0)
rseed < 0
rseed = 0
rseed > 0
irskip
(D=0)
乱数のコントロール
irskip > 0 irskip回ヒストリーをスキップする
irskip < 0 Irskip回乱数をスキップする
maxcas ×maxbch =全試行回数
• 初期設定では,常に
初期乱数オプション
同じ乱数を使う
計算開始時間から初期乱数を決定
• 毎回違う結果を得た
デフォルトの初期乱数
い場合は,rseed < 0
Rseedの値を初期乱数とする
とする
20
ヒストリーとバッチの関係
Q. ヒストリー数とは?
[source]セクションで指定した線源を発生させる回数→ぶどうの粒の数
Q. バッチ数とは?
ある一定のヒストリー数(maxcas)を束(バッチ)にして,その束を実行する回数(maxbch)
→ maxcas: 1房当たりのぶどうの粒の数 maxbch:ぶどうの房の数
Q. 各バッチの終了時には何をするの?
タリーの結果を取りまとめて平均値や統計誤差を導出し,途中経過を出力する (itall=1)
→収穫したぶどうの房の味見をする
Q. なぜ束(バッチ)に分けるの?
一気に全てのヒストリー数を実行すると,何か間違っていたときに再計算が大変!
→味見もせずに全てのぶどうを収穫すると,収穫時期を間違えたときに大損害!
各ヒストリー終了時に統計処理をすると,計算時間が掛かりすぎてしまう
→1粒1粒,味見をしながら収穫すると,収穫に時間が掛かりすぎてしまう
統計処理の時間が気にならない程度に,maxcasとmaxbchを調整する
例)1バッチ当たりの計算時間を2~3分にする
21
試行回数と統計誤差の関係
試行回数を増やすと統計誤差が小さくなる
lec03.inp
[Parameters]
icntl = 5
maxcas = 100
1000
10000
maxbch = 10
file(6) = phits.out
set: c1[30]
[Source]
s-type = 1
proj = proton
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 150
totfact = pi*c1**2
相対統計誤差
volume.out
・・・・
#num
1
2
3
4
5
・・
reg
101
102
103
104
105
volume
1.0000E+00
1.0000E+00
1.0000E+00
1.0000E+00
1.0000E+00
flux
3.6914E+02
4.8308E+02
5.2339E+02
3.5499E+03
3.6685E+03
3.6356E+03
1.0295E+04
1.0019E+04
9.9282E+03
1.9822E+04
1.9150E+04
1.9319E+04
3.1834E+04
3.1346E+04
3.1999E+04
r.err
0.0638
0.2362
0.0199
0.0965
0.0299
0.0095
0.0577
0.0184
0.0058
0.0380
0.0123
0.0039
0.0247
0.0080
0.0025
各層の体積は?
•4π(5)3/3=524 cm3
•4π(10)3/3 - 524=3665 cm3
•4π(15)3/3 - 4189=9948 cm3
•4π(20)3/3 - 14137=19373 cm3
•4π(25)3/3 - 33510=31940 cm3
22
再開始計算
lec03.inp
統計が足りない過去のタリー結果を読み込んで計算の再開が可能
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 5
maxcas = 10000
maxbch = 10
file(6) = phits.out
istdev = -1
再開始計算が実行された場合、コンソール画面に
メッセージが表示されます。
set: c1[30]
[Source]
s-type = 1
proj = proton
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 150
totfact = pi*c1**2
計算が終了した各バッチの情報も出力されます。
23
バッチ計算の活用
PHITSでは、タリーの結果を各バッチが終了する度に出力さ
せることができ、それを見ながらモンテカルロ計算を中断する
ことができます。
itall
(D=0) タリーの出力をバッチ毎に行うオプション
= 0 出力なし
= 1 同じファイルに出力
batch.now
1 <--- 1:continue, 0:stop
計算の途中で1を0に書き換え保存すると
その次のバッチで計算をやめる。
------------------------------------------------------------------------------bat[ 560] ncas =
560. : rijk = 151264979546685.
low neutron =
0. : ncall/s = 4.000000000E+00
cpu time = 0.288 s.
date = 2012-05-02
time = 15h 08m 25
24
lec03.inp
課題2: batch.nowを使って
計算を中断させてみましょう
file = lec03.inp
[Title]
・・・・・・
[Parameters]
icntl = 5
itall = 1
maxcas = 10000
maxbch = 100
file(6) = phits.out
$ istdev = -1
• maxbchを増やして長時間の計算を実行
• itallを使って途中結果を出力させる
(結果が刻々と変化するのを確認する)
• batch.nowを開いて中断
track_xz.eps
batch.now
01 <--- 1:continue, 0:stop
-------------------------------------------------------------------------------
途中のバッチで終わっているか
phits.outで確認してみましょう。
25
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
26
切断エネルギーとは?
• 注目していない、或いは重要でないと思われる粒子の輸送や核
反応の発生をストップし、数値計算の量を減らすことが可能です。
• 初期設定では輸送しない粒子(低エネルギー中性子,光子,電子
など)を輸送するように変更できます。
粒子
ID
切断エネルギー
陽子
emin(1)
1 MeV
中性子
emin(2)
1 MeV
π,μ
emin(3~8)
1 MeV
電子・陽電子
emin(12,13)
109 MeV
光子
emin(14)
109 MeV
原子核
emin(15~19)
1 MeV/u
27
通常の輸送計算を実行
lec03.inp
[Parameters]
icntl = 5
0
itall = 1
maxcas = 10000
1000
maxbch = 100
10
file(6) = phits.out
$ istdev = -1
set: c1[30]
[Source]
s-type = 1
proj = proton
・・・・・・
・・・・・・
・・・・・・
e0 = 150
輸送計算を実行する
修正して実行
[ T - C r o s s ] off
mesh = reg
reg = 5
r-in r-out area
102 101 1.0
103 102 1.0
104 103 1.0
105 104 1.0
106 105 1.0
e-type = 2
ne = 100
emin = 0.
emax = 200.
axis = eng
unit = 1
output = flux
file = cross.out
epsout = 1
玉ねぎ構造の各層を通過する粒子
をタリーする。
28
切断エネルギーの設定
lec03.inp
50MeV以下の中性子を輸送させない。
[Parameters]
icntl = 0
itall = 1
maxcas = 1000
maxbch = 10
emin(2) = 50
file(6) = phits.out
$ istdev = -1
修正して実行
set: c1[30]
cross.eps
50MeV未満の中性子はカット
→ 輸送されず各層を通過していない
29
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
30
断面積データライブラリとは?
• 光子・電子・陽電子並びに低エネルギー中性子(20MeV以下)の
断面積は,ターゲットやエネルギーに複雑に依存するため,画一
的なモデルでは表現できない → 断面積データライブラリが必要
• 断面積ライブラリを使った計算は,計算時間が掛かる場合がある
ため,デフォルトでは,それらの粒子の切断エネルギーは高く設
定されている
正しい輸送計算結果を得るため
には,断面積データライブラリを
使うように[parameters]セクショ
ンでeminなどを設定する必要が
ある
113Cdの中性子反応断面積(JENDL-4.0より)
31
データライブラリの設定
1. データライブラリファイルの確認
c:/phits/XS/neu
中性子用核データライブラリ
2. アドレス(xsdir)ファイルの確認
c:/phits/data/xsdir.jnd
アドレスファイル
3. アドレスファイルの1行目の確認
datapath=c:/phits/XS
データライブラリを含むフォルダ名
4. [Parameters]セクションにおいて“emin(i)”
“dmax(i)” ,“file(7)”パラメータを設定
以上はWindowsの場合です。Macの場合は、/Users/ユーザ名/phits/***の
ようになりますので、***の部分をそれぞれの内容に置き換えてください。
32
課題3: 中性子用核データライブラリの利用
lec03.inp
修正して実行
[Parameters]
icntl = 0
itall = 1
maxcas = 1000
中性子のデータを使う
maxbch = 10
(emin < dmaxを確認)
emin(2)= 1.0e-10
dmax(2) = 20
file(6) = phits.out
file(7) = c:/phits/data/xsdir.jnd
$ istdev = -1
アドレスファイル名
[T-Track]
mesh = reg
reg = (101 102 103 104 105)
e-type = 2
ne = 100
emin = 0
emax = 20.
axis = eng
unit = 1
part = proton neutron photon alpha
file = track_eng.out
epsout = 1
• 各領域の合計値を出力したい場
合は,合計したい領域を()で囲む
• 中性子の輸送を10-10(MeV)=0.1 (meV)まで計算する
• track_eng.eps から、20MeV以下の中性子が輸送されていることを確認する
33
課題3: 中性子用核データライブラリの利用
track_eng.eps
phits.out
--------------------------------------------------------------CPU time and number of event called in PHITS
--------------------------------------------------------------・・・・・・
・・・・・・
・・・・・・
dklos =
0.
hydro =
255.
n-data =
46556.
h-data =
0.
p-data =
0.
e-data =
0.
p-egs5 =
0.
e-egs5 =
0.
・・・・・・
・・・・・・
track_eng.epsで、中性子のスペクトルが見えた!
phits.outで、核データを使ったことが確認できる
課題4: 光子・電子用データライブラリの利用
• 電子・陽電子・光子の切断エネルギーemin(12~14) を設定する
→計算時間が長くなるので,ここでは1MeVに設定
• 電子・陽電子・光子のライブラリ上限エネルギーdmax(12~14)
を設定する → 全て1GeVに設定
PHITSを実行
• track_eng.epsで、光子スペクトルを確認
• phits.outにあるe-data,p-dataを確認
演習
[t-track](track_xz.out)のpartを変えて電子・光
子・陽電子の飛跡を確認してみよう
注)中性子の設定emin(2)&dmax(2)は,そのまま残しておいて下さい
35
課題5: EGS5の利用
lec03.inp
phits.out
[Parameters]
icntl = 0
itall = 1
maxcas = 1000
maxbch = 10
emin(2)= 1.0e-10
dmax(2) = 20
emin(12) = 1.0
emin(13) = 1.0
emin(14) = 1.0
dmax(12) = 1000.0
dmax(13) = 1000.0
EGS5用ライブラリ
dmax(14) = 1000.0
格納フォルダ名
file(6) = phits.out
file(7) = c:/phits/data/xsdir.jnd
file(20) = c:/phits/XS/egs
negs = 1
$ istdev = -1
光子・電子の輸送に
--------------------------------------------------------------CPU time and number of event called in PHITS
--------------------------------------------------------------・・・・・・
・・・・・・
・・・・・・
dklos =
112.
hydro =
252.
n-data =
43719.
h-data =
0.
p-data =
0.
e-data =
0.
p-egs5 =
1382.
e-egs5 =
29835.
・・・・・・
・・・・・・
EGS5を使う
(eminからdmaxまで)
36
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
37
γ脱励起オプション
igamma: 核反応で生成される放射性核種からのγ線放出
を考慮するためのオプション。デフォルトはγ線放出を考慮
しないので注意が必要
igamma
(D=0)
=0
=1
=2
=3
残留核のγ崩壊オプション
γ崩壊を考慮しない
γ崩壊を考慮する
γ崩壊をEBITEMモデルを用いて考慮する
γ崩壊とアイソマー生成をEBITEMモデルを
用いて考慮する
Version 2.64以降の奨励値は igamma = 2
38
荷電粒子のビームライン設計オプション
nspred & nedisp: 荷電粒子の物質中での角度・エネルギー
分散を考慮するためのオプション。加速器ビームのビーム
ライン設計の際,不可欠となる
nspred
(D=0)
=0
=1
=2
= 10
荷電粒子のクーロン散乱(角度分散)オプション
クーロン散乱を考慮しない
クーロン散乱をNMTCモデルを用いて考慮する
クーロン散乱をLynchの式を用いて考慮する
クーロン散乱をATIMAの式を用いて考慮する
nedisp
(D=0)
=0
=1
= 10
荷電粒子のエネルギー分散オプション
エネルギー分散を考慮しない
エネルギー分散をLandau Vavilovの式を用いて考慮する
エネルギー分散をATIMAの式を用いて考慮する
ビームライン設計時の奨励値は nspred = 2 & nedisp = 1
39
核反応モデルの変更
PHITSには、Bertini, JAM, JQMD, INCL, INC-ELFという
核反応モデルがあり、状況に応じてユーザーが使い分け
ることができます。
inclg
(D=1)
ejamnu
(D=20.)
JAMへの切り替えエネルギー(MeV)
ejampi
(D=20.)
パイオン入射反応のJAMへの切り替えエネル
ギー (MeV)
eqmdnu
(D=20.)
JQMDへの切り替えエネルギー (MeV)
eqmdmin
(D=10.)
JQMD適用の下限エネルギー (MeV/u)
ejamqmd
(D=3500.)
原子核反応のJQMDからJAMQMDへの切り替
えエネルギー (MeV/u)
incelf
(D=0)
INC-ELFを使用する場合の切り替えオプション
dmax(i)
(D=emin(i))
i-th粒子のライブラリー利用の上限エネルギー
INCLを使用する場合の切り替えオプション
40
核反応モデルの切替エネルギー
(1MeV)
emin(i)
Nucleon
(=emin)
(3.0GeV)
dmax(i)
核データ
einclmax
INCL (inclg=1)
(1MeV)
(3.0GeV)
emin(i)
Pion
einclmax
INCL (inclg=1)
ejamqmd
eqmdmin
Kaon, Hyperon
JAM
(3.5GeV/u)
(10MeV/u)
Nucleus
(d, t, 3He, α)
JAM
JQMD
INCL (inclg=1)
JAMQMD
JAM
41
イベントジェネレーターモード
核データ+統計崩壊モデルの利用により、残留核からの
放出粒子まで考慮し、かつエネルギーと運動量の保存則
を満たした核反応イベントを生成するモード。
イベントジェネレータモードを使った方がよい場合
• 低エネルギー中性子が放出する陽子やα粒子スペクトルが見たい
• 残留核の情報(反跳エネルギーなど)が知りたい
• イベント毎の情報が知りたい(検出器応答関数の計算など)
イベントジェネレータモードを使わない方がよい場合
• 中性子と光子の情報だけ知りたい(遮へい計算など)
• 中性子の透過率を計算したい
• 熱中性子の正確な挙動が見たい
42
イベントジェネレーターモードの設定
• 使用する方法は、核データを使用する設定を施した上で、
[parameters]セクションにおいて“e-mode=2”とするだけです
(指定しない場合は“e-mode=0”となっています)
• この場合“igamma”パラメータは自動的に2になります
e-mode
(D=0)
=0
=1
=2
Event generator modeオプション
通常のモード
Event generator mode version 1
Event generator mode version 2
43
課題6: イベントジェネレーターモードの利用
lec03.inp
[Parameters]
icntl = 0
itall = 1
maxcas = 1000
maxbch = 10
emin(2)= 1.0e-10
dmax(2) = 20
emin(12) = 1.0
emin(13) = 1.0
emin(14) = 1.0
dmax(12) = 1000.00000
dmax(13) = 1000.00000
dmax(14) = 1000.00000
file(6) = phits.out
…
修正して実行
[Source]
s-type = 1
proj = neutron
dir = 1.0
r0 = c1
z0 = -30.
z1 = -30.
e0 = 20.0
totfact = pi*c1**2
e-mode = 0
• まずはイベントジェネレーターモードを使わずに計算してみる。(デフォルト設定)
• 中性子の影響を見るため,線源は20MeV中性子とする
• 球全体の粒子フルエンスを出力して,陽子・α粒子が検出されたか確認する
44
課題6: イベントジェネレーターモードの利用
lec03.inp
修正して実行
[Parameters]
icntl = 0
itall = 1
maxcas = 1000
maxbch = 10
emin(2)= 1.0e-10
dmax(2) = 20
emin(12) = 1.0
emin(13) = 1.0
emin(14) = 1.0
dmax(12) = 1000.00000
dmax(13) = 1000.00000
dmax(14) = 1000.00000
file(6) = phits.out
file(7) = c:/phits/data/xsdir.jnd
...
e-mode = 2
track_eng.eps
陽子やα粒子のスペクトルが見えた!
45
入出力ファイルの設定
サマリーの出力ファイル名。指定しない時は、
標準出力
File(6)
(D=phits.out)
File(7)
(D=c:/phits/data/xsdir.jnd)
File(15)
(D=dumpall.dat)
File(18)
(D=voxel.bin)
Ivoxel=1,2を指定した際のバイナリデータ用
ファイル名
File(20)
(D=c:/phits/XS/egs/)
EGS用断面積データを格納したフォルダ名
(negs=1の際に必要)
File(21)
断面積のアドレスファイル名
Dumpall=1を指定した時の出力ファイル名
(D=c:/phits/dchain-sp/data/) DCHAIN-SP用データを格納したフォルダ名
46
実習内容
 計算モードの選択
 便利な入力補助機能
 統計的に信頼できる結果を得るための設定
 モンテカルロ積分を使った体積・面積計算
 試行回数と統計誤差
 物理的に信頼できる結果を得るための設定
 切断エネルギーの設定
 断面積データライブラリの利用
 物理モデルの選択
 まとめ
47
まとめ
• PHITSでは、[Parameters]セクションの設定を変更することで、
様々な計算モードや物理モデルの選択を行うことができる。
• 各計算モード(icntl)を使い分けることにより、ジオメトリの確
認,ソースのチェック、粒子の輸送計算を行う。
• 統計精度と試行回数(maxcas, maxbch)は密接に関連して
おり、適切な数を設定する必要がある。
• 低エネルギーの中性子や光子,電子,陽電子等は,切断エ
ネルギー(emin)や核データの上限エネルギー(dmax)を変
更することで精度の高い計算結果を得ることができる。その
際,file(7)でアドレスファイルを指定する必要がある
• 状況に応じてイベントジェネレーターモード(e-mode)など物理
モデルの設定を変更する必要がある。
最適な設定が分からない場合は,奨励設定ファイル(recommendation)をご参照ください
48
宿題
• 宿題ファイルで,熱中性子まで考慮し,
20MeV以下の中性子に対して核データラ
イブラリを使う
• イベントジェネレータモードを使ってみる
• 線量-深さ分布の相対標準偏差2%以内を
目指す(maxcas, istdev, batch.nowな
どを活用する)
49
宿題(解答例)
陽子(上)・中性子(下)フルエンス
円柱中心(上)と端側(下)における
付与エネルギーの深さ分布
50