PHITS 講習会 基礎実習（III）: 計算条件の設定 Multi-Purpose Particle and Heavy Ion Transport code System PHITS 講習会　基礎実習（III）: 計算条件の設定 2016年8月改訂 title 1

本実習では，パラメータの設定方法について学習します はじめに PHITSの計算は，[Parameters]セクションで定義さ れる様々なパラメータ（物理モデルの選択など）によ りコントロールされています。 全てのパラメータには初期設定値があり，基本的に は変更する必要はありません。 ただし，最適な計算結果を得るためには，着目して いる輸送粒子の条件に応じて、いくつかのパラメー タを変更する必要があります。 本実習では，パラメータの設定方法について学習します 2

本実習の目標 PHITSの便利な機能を活用し，統計量や物理モデルを変化させ，最適な計算結果が得られるようになる 3 宿題体系における初期設定での 陽子（上）・中性子（下）フルエンス ヒストリー数，切断エネルギー，核データを適切 に設定した陽子（上）・中性子（下）フルエンス 3

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 4

計算モードの選択 通常輸送計算 反応なしモード 体系確認モード 5

計算モードの選択 体系確認モード まずは、基礎実習(II)で行ったジオメトリを確認するタリーを用いて、サンプルインプットの体系を確認してみましょう。 6

体系の確認([t-3dshow]) lec03.inp 7 1層5cm幅の玉ねぎ構造 3dshow.eps [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 file = 3dshow.out title = Check onion structure using [T-3dshow] tally epsout = 1 ・ ・ ・ ・ ・ ・ file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 11 maxcas = 100 maxbch = 10 file(6) = phits.out [t-3dshow]を実行 1層5cm幅の玉ねぎ構造 3dshow.eps 7

体系の確認(gshowオプション) lec03.inp 8 track_xz.eps file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 11 maxcas = 100 maxbch = 10 file(6) = phits.out file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 8 maxcas = 100 maxbch = 10 file(6) = phits.out [ T - T r a c k ] 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 9 修正して実行（反応なしモード） [ T - T r a c k ] 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 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out 9

全ての物質がvoidとなるため，まっすぐに飛んでいく 線源の飛跡確認 track_xz.eps 全ての物質がvoidとなるため，まっすぐに飛んでいく （線源の発生位置・方向が確認できる） 10

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 11

外部ファイルの挿入 lec03.inp 12 インプットファイルとは別のファイルを取り込む場合は”infl:”コマンドを使う onion.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[5] [ S o u r c e ] infl: {onion.inp}[1-33] [ M a t e r i a l ] ・ ・ ・ ・ ・ ・ [ M a t N a m e C o l o r ] [ S u r f a c e ] 10 so 500. 11 so 5. 12 so 10. 13 so 15. 14 so 20. 15 so 25. [ C e l l ] 100 -1 10 101 1 -19.32 -11 102 2 -1. 11 -12 103 3 -8.93 12 -13 104 4 -1. 13 -14 105 5 -0.9 14 -15 106 6 -1.20e-3 15 -10 インプットファイルとは別のファイルを取り込む場合は”infl:”コマンドを使う 　　{}でファイル名，[n1-n2]で取り込む行の範囲を指定する インプットファイルの一行目に“file = インプットファイル名”を書く(重要) “off”したセクションの下にある場合、無視されるので注意 12

ユーザー定義定数と数式の利用 lec03.inp 13 ユーザー定義定数の利用： 同じ値を使う場合は，予め定義しておくと便利！ 数式の利用： file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] set: c1[5] [ S o u r c e ] totfact = 1.0 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 PHITSでは定数を定義して、インプットファイル内でその定数を参照することが可能。 “set: ci[x]” 名前はc1からc99まで使用可能 書いた箇所より下の部分で有効となる。また、何度も同じ名前を定義できる “off”で飛ばしたセクションの下にある場合、無視されるので注意 数式の利用： PHITS入力ファイルでは，FORTRAN形式の数式が利用可能 円周率（pi）は定義しなくても使える べき乗は「**」なので注意（例えば、p(c1)2はpi*c1**2） cosやexpなどの関数も使える（例えばcos(pi/2)。次頁に関数一覧） 13

利用できる関数一覧 14 関数 説明 sin(x) 正弦 cos(x) 余弦 tan(x) 正接 asin(x) 逆正弦 acos(x) 逆余弦 atan(x) 逆正接 sinh(x) 双曲線正弦 cosh(x) 双曲線余弦 tanh(x) 双曲線正接 関数 説明 int(x) 整数部分 nint(x) 四捨五入した整数値 abs(x) 絶対値 exp(x) eを底とする指数関数 log(x) 自然対数 log10(x) 常用対数 sqrt(x) 平方根 max(x,y,…) 最大値 min(x,y,…) 最小値 mod(x,y) xをyで割った余り 14

課題1 玉ねぎ体系全体を覆うように、ビームの幅を拡げてみましょう。 ユーザー定義定数c1の値を大きくする lec03.inp 15 角柱形状の線源 (s-type=2) set: c1[5] [ S o u r c e ] totfact = 1.0 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 x軸とy軸に関して1辺10cmの正方形を発生領域とし、+z軸をビームの方向とする線源 track_xz.eps 半径25cmの球を覆うことができるように、ビームの幅を拡げてみよう 15

課題1の答え合わせ 玉ねぎ体系全体を覆うように、ビームの幅を拡げてみましょう。 lec03.inp 16 角柱形状の線源 (s-type=2) set: c1[25] [ S o u r c e ] totfact = 1.0 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 track_xz.eps 発生領域となる正方形の1辺の長さを50cmをとすることで玉ねぎ体系を覆うことが可能となった 16

規格化について lec03.inp 17 PHITSのタリー結果は，発生する線源数*（試行 回数[source]）あたりで出力される 　　（*厳密には，線源粒子のウェイトあたり） 測定値などと比較するには，Bqや/cm2/sで表され る線源の発生頻度で規格化する必要がある タリー結果を定数倍（totfactで定義）することによ り，直接， /cm2/sやmGy/hなどの単位に換算した 値を出力できるようになる set: c1[25] [ S o u r c e ] totfact = 1.0 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 17

課題2 タリーで得られる値を1,000倍にしてみましょう。 計算結果に掛けたい値をtotfactに設定する lec03.inp 18 set: c1[25] [ S o u r c e ] totfact = 1.0 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 10-3となっている目盛りを1に増やすにはどうすればよいだろう？ track_xz.eps 18

課題2の答え合わせ タリーで得られる値を1,000倍にしてみましょう。 lec03.inp 19 set: c1[25] [ S o u r c e ] totfact = 1000 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 タリー結果全体が1000倍され、目盛りが100(=1)となっている track_xz.eps 19

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 20

体積、面積計算の必要性 reg=meshとして各領域毎の物理量を計算するタリーの場合、その領域の体積や表面積を計算しておく必要があります*。 指定した領域の 単位体積あたりの熱量は？ 指定した表面の 単位面積あたりの粒子フルエンスは？ 指定した領域が球や円柱などであれば解析的に計算可能。 しかし、複雑な形状の場合はどうすれば良いだろうか？ →　PHITSを用いてモンテカルロ積分を行う ＊各領域の体積は[volume]セクションで、表面積は[t-cross]の”area”で与えることができます。 21

モンテカルロ積分 乱数を用いて定積分の近似解を求める方法。積分範囲の内か外かを判定することにより、積分値を確率的に概算する。計算精度は試行回数に依存する。 a cm 1辺a cmの正方形の的を用意し、その的の内部にランダムに点を打って、色の付いた部分に当たれば+1点とする。 左の例の場合、1/2の確率で+1点されるので、総得点を試行回数で割った値は，試行回数を増やせば1/2に近づいていく。これに的の面積a2を掛けると、色の付いた部分の面積（a2/2）を概算することができる。 a cm 総得点を試行回数で割り、その値に的の面積を掛けることで任意の面積を概算することができます。また、この考え方は3次元にも拡張でき、同様に体積も概算できます 22

（線源領域に関する要素はtotfactを使って結果に反映させる） PHITSによる体積計算の原理 PHITSでは[t-track]によって得られるtrack lengthを用いることで体積の概算が可能となる。 求める体積 [cm3] ＝線の長さの合計（総得点）[cm] 　÷線の数(試行回数) 　×線源領域（的の面積）[cm2] ＝ 平均track length 　 ×線源領域（的の面積）[cm2] [t-track]の結果 （線源領域に関する要素はtotfactを使って結果に反映させる） 線源領域 的 的の面積と線源領域が同じとなる条件で計算を実行する 23

課題3 玉ねぎ体系の各層の体積を計算してみよう。 線源領域（的の面積）をc1を使ってtotfactに設定する file=volume.outとなっている[t-track]を有効にする lec03.inp set: c1[25] [ S o u r c e ] totfact = 1000 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 [ T - T r a c k ] off mesh = reg reg = 101 102 103 104 105 ・ ・ ・ ・ ・ ・ file = volume.out 各層の体積はどの位になるだろう？ 24

課題3の答え合わせ 玉ねぎ体系の各層の体積を計算してみよう。 lec03.inp 25 体系の体積や面積を計算する場合: icntl=5, set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 [ T - T r a c k ] mesh = reg reg = 101 102 103 104 105 ・ ・ ・ ・ ・ ・ file = volume.out 体系の体積や面積を計算する場合: icntl=5, 線源領域を計算したい体系が全部入る大きさにする volume.outを開いて計算結果を確認してみよう！ 25

体積計算の結果 26 モンテカルロ積分の要領で体系の体積を求めることが可能 volume.out 各層の体積は？ ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 6.3277E+02 0.1630 2 102 1.0000E+00 3.6913E+03 0.0897 3 103 1.0000E+00 9.3260E+03 0.0561 4 104 1.0000E+00 1.8670E+04 0.0361 5 105 1.0000E+00 3.2614E+04 0.0211 各層の体積は？ 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 内側の結果ほど一致していない →　統計量が十分ではない →　統計量を増やすには？ 26

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 27

試行回数（ヒストリ数）の設定 モンテカルロ計算の統計精度は、シミュレーションの試行回数に依存しています。したがって、信頼できる計算結果を得るためには、十分な数の試行回数を設定する必要があります。 maxcas (D=10) 1バッチのヒストリー数 maxbch バッチ数 rseed (D=0.0) rseed < 0 rseed = 0 rseed > 0 初期乱数オプション 計算開始時間から初期乱数を決定 デフォルトの初期乱数 rseedの値を初期乱数とする irskip (D=0) irskip > 0 irskip < 0 乱数のコントロール irskip回ヒストリーをスキップする irskip回乱数をスキップする maxcas ×maxbch ＝全試行回数 初期設定では，常に同じ乱数を使う 毎回違う結果を得たい場合は，rseed < 0とする 28

統計処理の時間が気にならない程度に，maxcasとmaxbchを調整する ヒストリーとバッチの関係 Q. ヒストリー数とは？ 　　[source]セクションで指定した線源を発生させる回数→ぶどうの粒の数 Q. バッチ数とは？ 　　ある一定のヒストリー数（maxcas）を束（バッチ）にして，その束を実行する回数（maxbch） 　　→　maxcas: 1房当たりのぶどうの粒の数　　maxbch：ぶどうの房の数 Q. 各バッチの終了時には何をするの？ 　　タリーの結果を取りまとめて平均値や統計誤差を導出し，途中経過を出力する (itall=1) 　　→収穫したぶどうの房の味見をする Q. なぜ束（バッチ）に分けるの？ 　　一気に全てのヒストリー数を実行すると，何か間違っていたときに再計算が大変！ 　　　→味見もせずに全てのぶどうを収穫すると，収穫時期を間違えたときに大損害！ 　　各ヒストリー終了時に統計処理をすると，計算時間が掛かりすぎてしまう 　　　→１粒１粒，味見をしながら収穫すると，収穫に時間が掛かりすぎてしまう 統計処理の時間が気にならない程度に，maxcasとmaxbchを調整する 例）1バッチ当たりの計算時間を2～3分にする 29

課題4 統計精度を高めるために、試行回数を増やしてみましょう。 lec03.inp maxcasを1,000や10,000に増やす 統計誤差が小さくなっていることを確認する [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 統計誤差 （相対誤差） volume.out ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 6.3277E+02 0.1630 2 102 1.0000E+00 3.6913E+03 0.0897 3 103 1.0000E+00 9.3260E+03 0.0561 4 104 1.0000E+00 1.8670E+04 0.0361 5 105 1.0000E+00 3.2614E+04 0.0211 試行回数を10倍にすると誤差の値はどの程度小さくなるか？ 30

課題4の答え合わせ 統計精度を高めるために、試行回数を増やしてみましょう。 lec03.inp 31 試行回数と統計誤差の関係 volume.out ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 6.3277E+02 0.1630 2 102 1.0000E+00 3.6913E+03 0.0897 3 103 1.0000E+00 9.3260E+03 0.0561 4 104 1.0000E+00 1.8670E+04 0.0361 5 105 1.0000E+00 3.2614E+04 0.0211 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.2623E+02 0.0187 2 102 1.0000E+00 3.6038E+03 0.0089 3 103 1.0000E+00 9.8977E+03 0.0054 4 104 1.0000E+00 1.9304E+04 0.0035 5 105 1.0000E+00 3.2052E+04 0.0021 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.2586E+02 0.0588 2 102 1.0000E+00 3.6515E+03 0.0279 3 103 1.0000E+00 9.7779E+03 0.0172 4 104 1.0000E+00 1.9166E+04 0.0111 5 105 1.0000E+00 3.2095E+04 0.0067 [ P a r a m e t e r s ] icntl = 5 maxcas = 10000 maxbch = 10 file(6) = phits.out set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 [ P a r a m e t e r s ] icntl = 5 maxcas = 1000 maxbch = 10 file(6) = phits.out set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 正解 （解析解） 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 試行回数と統計誤差の関係 試行回数を増やすと正解に近づく 試行回数を10倍する毎に、誤差はおよそ1/√10になる 31

再開始計算 lec03.inp 32 統計量が足りない場合に、過去のタリー結果を読み込んで、続き計算を行うことが可能です。 istdev=-1を加えて、再開始計算を実行してみましょう！ [ P a r a m e t e r s ] icntl = 5 maxcas = 10000 maxbch = 10 file(6) = phits.out istdev = -1 set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 再開始計算が実行された場合、コンソール画面にメッセージが表示されます。 計算が終了した各バッチの情報も出力されます。 32

バッチ計算の活用 PHITSでは、タリーの結果を各バッチが終了する度に出力させることができ、それを見ながらモンテカルロ計算を中断することができます。 itall (D=0) = 0 = 1 タリーの出力をバッチ毎に行うオプション 出力なし 同じファイルに出力 batch.out 計算の途中で1を0に書き換え保存すると その次のバッチで計算をやめる。 1 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- 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 33

課題5 batch.outを使って、計算を中断させてみましょう lec03.inp 34 istdevの行は$マークでコメントアウト（再開始計算は行わない） itallを使って途中結果を出力させる 　　（結果が刻々と変化するのを確認する） maxbchを100に増やして長時間の計算を実行 batch.outを開いて中断 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 itall = maxcas = 10000 maxbch = 10 file(6) = phits.out istdev = -1 34

課題5の答え合わせ batch.outを使って、計算を中断させてみましょう lec03.inp 35 PHITSを実行する track_xz.epsを更新し、結果が刻々と変化するのを確認する batch.outを開いて中断 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 itall = 1 maxcas = 10000 maxbch = 100 file(6) = phits.out $ istdev = -1 batch.out 0 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- 1 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- track_xz.eps 途中のバッチで終わっているかphits.outを開いて確認してみましょう。 35

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 36

切断エネルギーとは？ 注目していない、或いは重要でないと思われる粒子の輸送や核反応の発生をストップし、数値計算の量を減らすことが可能です。 初期設定では輸送しない粒子（1MeV未満の低エネルギー中性子，光子，電子など）を輸送するように変更できます。 粒子 パラメータ 切断エネルギー 陽子 emin(1) 1 MeV 中性子 emin(2) p+, p0, p- emin(3), emin(4), emin(5) m+, m- emin(6), emin(7) K+, K0, K- emin(8), emin(9), emin(10) 電子, 陽電子, 光子 emin(12), emin(13), emin(14) 109 MeV d, t, 3He, a emin(15), emin(16), emin(17), emin(18) 1 MeV/u 原子核 emin(19) 37

通常の輸送計算を実行 lec03.inp 38 物質を考慮した輸送計算を実行し、玉ねぎ構造の各層を通過する粒子の振る舞いを調べましょう [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 file(6) = phits.out $ istdev = -1 修正して実行 [ 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 ・ ・ ・ ・ ・ ・ offを消して [t-cross]を 有効にする タリー面として、どの領域の境界面が指定されているか確認しよう！ [T-cross]とは、 指定した面を通過する粒子を計算するタリー mesh=regの場合、r-inの領域からr-outの領域に入射した粒子をタリーする 38

通常の輸送計算を実行 lec03.inp 39 物質を考慮した輸送計算を実行し、玉ねぎ構造の各層を通過する粒子の振る舞いを調べましょう [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 file(6) = phits.out $ istdev = -1 cross.eps（3ページ目） 3ページ目は、領域104から領域103に入射した粒子のエネルギー分布。陽子は、約85MeV以下の領域で連続的な分布をとっている。 1,2ページ目はどうなっているだろう？ track_xz.eps 39

課題6 中性子の切断エネルギーを設定し、50MeV以下の中性子を輸送させない計算を実行してみましょう lec03.inp 40 [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = file(6) = phits.out $ istdev = -1 emin(2)を追加で設定する（中性子の指定する番号は2） 切断エネルギーの設定により、50MeV以下の中性子はどうなるだろうか？ cross.eps（3ページ目） 40

課題6の答え合わせ 中性子の切断エネルギーを設定し、50MeV以下の中性子を輸送させない計算を実行してみましょう lec03.inp 41 [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 50 file(6) = phits.out $ istdev = -1 cross.eps（3ページ目） 50MeV未満の中性子はカットされている →輸送されず各層を通過していない 41

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 42

断面積データライブラリとは？ 光子・電子・陽電子並びに低エネルギー中性子（20MeV以下）の断面積は，ターゲットやエネルギーに複雑に依存するため，画一的なモデルでは表現できない → 断面積データライブラリが必要 断面積ライブラリを使った計算は，計算時間が掛かる場合があるため，デフォルトでは，それらの粒子の切断エネルギーは高く設定されている 正しい輸送計算結果を得るためには，断面積データライブラリを使うように[parameters]セクションでeminなどを設定する必要がある 113Cdの中性子反応断面積（JENDL-4.0より） 43

データライブラリの設定 データライブラリファイルの確認 アドレス（xsdir）ファイルの確認 アドレスファイルの１行目の確認 c:/phits/XS/neu アドレス（xsdir）ファイルの確認 c:/phits/data/xsdir.jnd アドレスファイルの１行目の確認 datapath=c:/phits/XS [Parameters]セクションにおいて“emin(i)” “dmax(i)” ，“file(7)”パラメータを設定 中性子用核データライブラリ アドレスファイル データライブラリを含むフォルダ名 以上はWindowsの場合です。Macの場合は、/Users/ユーザ名/phits/***のようになりますので、***の部分をそれぞれの内容に置き換えてください。このユーザ名がわからない場合は、xsdirの1行目にあるdatapathを参照してください。 44

課題7 中性子用核データライブラリを利用して、 0.1meV以上の中性子を考慮した計算を実行しましょう lec03.inp 45 [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 50 dmax(2) = file(6) = phits.out file(7) = $ istdev = -1 emin(2)を1.0e-10とする（10-10[MeV]=0.1[meV]） dmax(2)を追加で設定し20MeVとする（付属している中性子のデータライブラリの上限は20MeV） file(7)=c:/phits/data/xsdir.jndを追加する（Windowsの場合。環境に応じたアドレスファイル名を指定する。） file=track_eng.outとなっている[t-track]のoffを消す [ T - T r a c k ] off mesh = reg reg = (101 102 103 104 105) e-type = 2 ・ ・ ・ ・ ・ ・ axis = eng unit = 1 part = proton neutron photon alpha file = track_eng.out epsout = 1 玉ねぎ構造5層の領域の平均値 （）で領域をまとめて、そのまとめた領域の平均値を求める 45

課題7の答え合わせ1 中性子用核データライブラリを利用して、 0.1meV以上の中性子を考慮した計算を実行しましょう lec03.inp [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 1.0e-10 dmax(2) = 20 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd $ istdev = -1 [ T - T r a c k ] mesh = reg reg = (101 102 103 104 105) e-type = 2 ・ ・ ・ ・ ・ ・ axis = eng unit = 1 part = proton neutron photon alpha file = track_eng.out epsout = 1 track_eng.eps 1MeV以下の中性子も輸送されている！ 46

課題7の答え合わせ2 中性子用核データライブラリを利用して、 0.1meV以上の中性子を考慮した計算を実行しましょう lec03.inp [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 1.0e-10 dmax(2) = 20 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd $ istdev = -1 phits.out（の最後の部分。520行目辺り） --------------------------------------------------------------- CPU time and number of event called in PHITS ・ ・ ・ ・ ・ ・ dklos = 0. hydro = 288. n-data = 42958. h-data = 0. p-data = 0. e-data = 0. p-egs5 = 0. e-egs5 = 0. [ T - T r a c k ] mesh = reg reg = (101 102 103 104 105) e-type = 2 ・ ・ ・ ・ ・ ・ axis = eng unit = 1 part = proton neutron photon alpha file = track_eng.out epsout = 1 phits.outで、核データを使ったことが確認できる 47

課題8 光子・電子用データライブラリを利用して、光子と電子を考慮した輸送計算を実行しましょう lec03.inp 48 [ P a r a m e t e r s ] icntl = 0 　itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 1.0e-10 dmax(2) = 20 emin(12) = emin(13) = emin(14) = dmax(12) = dmax(13) = dmax(14) = file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd $ istdev = -1 電子・陽電子・光子の切断エネルギーemin(12),emin(13),emin(14)を1MeVに設定する 電子・陽電子・光子に関するライブラリの上限エネルギーdmax(12),dmax(13),dmax(14)を1GeVとする 注）中性子の設定emin(2)&dmax(2)は，そのまま残しておいて下さい PHITSを実行し、track_eng.epsの図で光子のスペクトルを確認しましょう！また、phits.outの最後の欄でp-dataとe-dataの値を確認しましょう！ 追加演習 [t-track](track_xz.out)のpartを変更し、中性子・電子・陽電子・光子の飛跡を確認して見ましょう！ 48

課題8の答え合わせ 光子・電子用データライブラリを利用して、光子と電子を考慮した輸送計算を実行しましょう lec03.inp 49 [ P a r a m e t e r s ] 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 dmax(14) = 1000.0 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd $ istdev = -1 track_eng.eps 10MeV以下の光子(photon)が輸送されている！ 49

課題9 光子・電子の輸送にEGS5を利用しましょう lec03.inp 50 negs=1を追加する [ P a r a m e t e r s ] 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 dmax(14) = 1000.0 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd file(20) = negs = $ istdev = -1 negs=1を追加する file(20)=c:/phits/XS/egsを追加する（Windowsの場合。Macは各環境に応じて変更してください。） PHITSを実行し、phits.outの最後の欄でp-egs5とe-egs5の値を確認しましょう！ 50

課題9の答え合わせ 光子・電子の輸送にEGS5を利用しましょう lec03.inp 51 phits.out（の最後の部分。540行目辺り） [ P a r a m e t e r s ] 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 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 phits.out（の最後の部分。540行目辺り） --------------------------------------------------------------- CPU time and number of event called in PHITS ・ ・ ・ ・ ・ ・ dklos = 118. hydro = 289. n-data = 43949. h-data = 0. p-data = 0. e-data = 0. p-egs5 = 1477. e-egs5 = 42529. 51

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 52

γ脱励起オプション igamma: 核反応で生成される放射性核種からのγ線放出を考慮するためのオプション。デフォルトはγ線放出を考慮しないので注意が必要 igamma (D=0) = 0 = 1 = 2 = 3 残留核のγ崩壊オプション γ崩壊を考慮しない γ崩壊を考慮する γ崩壊をEBITEMモデルを用いて考慮する γ崩壊とアイソマー生成をEBITEMモデルを用いて考慮する Version 2.64以降の奨励値は igamma = 2 53

荷電粒子のビームライン設計オプション nspred & nedisp: 荷電粒子の物質中での角度・エネルギー分散を考慮するためのオプション。加速器ビームのビームライン設計の際，不可欠となる nspred (D=0) = 0 = 1 = 2 = 10 荷電粒子のクーロン散乱（角度分散）オプション クーロン散乱を考慮しない クーロン散乱をNMTCモデルを用いて考慮する クーロン散乱をLynchの式を用いて考慮する クーロン散乱をATIMAの式を用いて考慮する nedisp 荷電粒子のエネルギー分散オプション エネルギー分散を考慮しない エネルギー分散をLandau Vavilovの式を用いて考慮する エネルギー分散をATIMAの式を用いて考慮する ビームライン設計時の奨励値は nspred = 2 & nedisp = 1 54

核反応モデルの変更 PHITSには、Bertini, JAM, JQMD, INCL, INC-ELFという核反応モデルがあり、状況に応じてユーザーが使い分けることができます。 inclg (D=1) INCLを使用する場合の切り替えオプション ejamnu (D=20.) JAMへの切り替えエネルギー(MeV) ejampi パイオン入射反応のJAMへの切り替えエネルギー (MeV) eqmdnu 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粒子のライブラリー利用の上限エネルギー 55

核反応モデルの切替エネルギー Nucleon 核データ INCL (inclg=1) JAM Pion INCL (inclg=1) JAM (1MeV) emin(i) (=emin) dmax(i) (3.0GeV) einclmax Nucleon 核データ INCL (inclg=1) JAM (1MeV) emin(i) (3.0GeV) einclmax Pion INCL (inclg=1) JAM (10MeV/u) eqmdmin (3.5GeV/u) ejamqmd Nucleus JQMD JAMQMD (d, t, 3He, α) INCL (inclg=1) Kaon, Hyperon JAM 56

イベントジェネレーターモード 核データ＋統計崩壊モデルの利用により、残留核からの放出粒子まで考慮し、かつエネルギーと運動量の保存則を満たした核反応イベントを生成するモード。 イベントジェネレータモードを使った方がよい場合 低エネルギー中性子が放出する陽子やα粒子スペクトルが見たい 残留核の情報（反跳エネルギーなど）が知りたい イベント毎の情報が知りたい（検出器応答関数の計算など） イベントジェネレータモードを使わない方がよい場合 中性子と光子の情報だけ知りたい（遮へい計算など） 中性子の透過率を計算したい 熱中性子の正確な挙動が見たい 57

イベントジェネレーターモードの設定 使用する方法は、核データを使用する設定を施した上で、[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 58

課題10 まず、イベントジェネレータモードを使用せずに、20MeV中性子を線源とする計算を実行してみましょう lec03.inp 59 set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = proton e0 = 150.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 線源粒子を中性子に変更する 線源エネルギーを20MeVとする PHITSを実行し、track_eng.epsの図で陽子やアルファ粒子が検出されるかどうかを確認しましょう！ 59

課題10の答え合わせ まず、イベントジェネレータモードを使用せずに、20MeV中性子を線源とする計算を実行してみましょう lec03.inp set: c1[25] [ S o u r c e ] totfact = (c1*2)**2 s-type = 2 proj = neutron e0 = 20.0 x0 = -c1 x1 = c1 y0 = -c1 y1 = c1 z0 = -30.0 z1 = -30.0 dir = 1.0 track_eng.eps 陽子やアルファ粒子は検出されていない！ →イベントジェネレータモードを使用しない場合、20MeV中性子の反応で荷電粒子は生成されない 60

課題11 イベントジェネレータモードを使用してみましょう lec03.inp 61 [parameters]セクションでe-modeを2として設定する [ P a r a m e t e r s ] 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 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 e-mode = イベントジェネレータモードの有無で、陽子やアルファ粒子が検出されるかどうかを確認しましょう！ 61

課題11の答え合わせ イベントジェネレータモードを使用してみましょう lec03.inp 62 陽子やアルファ粒子が生成されている！ [ P a r a m e t e r s ] 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 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 e-mode = 2 track_eng.eps 陽子やアルファ粒子が生成されている！ 62

入出力ファイルの設定 63 file(6) (D=phits.out) サマリーの出力ファイル名。指定しない時は、標準出力 file(7) (D=c:/phits/data/xsdir.jnd) 断面積のアドレスファイル名 file(15) (D=dumpall.dat) dumpall=1を指定した時の出力ファイル名 file(18) (D=voxel.bin) ivoxel=1,2を指定した際のバイナリデータ用ファイル名 file(20) (D=c:/phits/XS/egs/) EGS用断面積データを格納したフォルダ名 （negs=1の際に必要） file(21) (D=c:/phits/dchain-sp/data/) DCHAIN-SP用データを格納したフォルダ名 file(22) (D=batch.out) 現在のバッチ情報を出力するファイル名 file(23) (D=pegs5) PEGS5用出力ファイル名 63

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 64

課題12 100MeVの陽子ビームを玉ねぎ球に一様に照射する場合を考えます。ビームの流束が1cm2・1秒あたりに1個である場合、ビームを1時間照射した場合の、各層に付与されるエネルギーを計算しましょう PHITSで出力される線源粒子あたりのエネルギーGy/sourceを、Gy単位の物理量に変換するための規格化定数の問題です。 100MeV陽子はこれまでと同じ1辺50cmの正方形の領域から発生すると考える（線源領域は変更しない） unit=0となっている[t-deposit]を有効にする（unit=0は、Gy/sourceを単位とする結果を出力する） [source]セクションにおいて、線源粒子とエネルギーを変更し、Gy/sourceをGyに変換するための係数をtotfactに設定する 100MeVの陽子が1cm2・1秒あたり1個の割合で照射される （ヒント） 1cm2・1秒あたり1個のビームが50cm四方の領域から発生し、それを1時間照射する場合を考えると、発生する陽子数は合計何個[source]になるでしょうか？ deposit_reg.outを開いて計算結果を確認してみよう！ 65

課題12の答え合わせ 100MeVの陽子ビームを玉ねぎ球に一様に照射する場合を考えます。ビームの流束が1cm2・1秒あたりに1個である場合、ビームを1時間照射した場合の、各層に付与されるエネルギーを計算しましょう totfactの計算： 1cm2あたり1個なので、50cm四方の領域から発生する1秒あたりの陽子数は50*50個 よって、このビームを1時間照射すると50*50*3600個[source]の陽子が発生する [t-deposit]でunit=0とした場合は[Gy/source]で出力されるので、totfactを用いて50*50*3600[source]を掛けることで、タリー結果が[Gy]を意味するようになる lec03.inp [ S o u r c e ] totfact = (c1*2)**2*3600 s-type = 2 proj = proton e0 = 100.0 x0 = -c1 ・ ・ ・ ・ ・ ・ [ T-deposit ] title = Deposition energy ・ ・ ・ mesh = reg reg = 101 102 103 104 105 unit = 0 offを消す deposit_reg.out ・ ・ ・ ・ ・ ・ # num reg volume dose r.err 1 101 5.2400E+02 2.8336E-10 0.4636 2 102 3.6650E+03 3.2418E-09 0.2985 3 103 9.9480E+03 9.0401E-10 0.1473 4 104 1.9373E+04 1.0944E-06 0.0147 5 105 3.1940E+04 3.0658E-06 0.0063 66

実習内容 計算モードの選択 便利な入力補助機能 統計的に信頼できる結果を得るための設定 物理的に信頼できる結果を得るための設定 モンテカルロ積分を使った体積・面積計算 試行回数と統計誤差 物理的に信頼できる結果を得るための設定 切断エネルギーの設定 断面積データライブラリの利用 物理モデルの選択 規格化に関する問題 まとめ 67

最適な設定が分からない場合は，奨励設定ファイル(recommendation)をご参照ください まとめ PHITSでは、[Parameters]セクションの設定を変更することで、 様々な計算モードや物理モデルの選択を行うことができる。 各計算モード（icntl）を使い分けることにより、ジオメトリの確 認，ソースのチェック、粒子の輸送計算を行う。 統計精度と試行回数(maxcas, maxbch)は密接に関連して おり、適切な数を設定する必要がある。 低エネルギーの中性子や光子，電子，陽電子等は，切断エ ネルギー（emin）や核データの上限エネルギー（dmax）を変 更することで精度の高い計算結果を得ることができる。その 際，file(7)でアドレスファイルを指定する必要がある 状況に応じてイベントジェネレーターモード(e-mode)など物理 モデルの設定を変更する必要がある。 《休憩はさむ》 まとめ 最適な設定が分からない場合は，奨励設定ファイル(recommendation)をご参照ください 68

宿題 宿題ファイルで，熱中性子まで考慮し，20MeV以下の中性子に対して核データライブラリを使う イベントジェネレータモードを使ってみる 線量-深さ分布の相対標準偏差２％以内を目指す（maxcas, istdev, batch.outなどを活用する） 宿題（3次元体系） 69

宿題（解答例） 陽子（上）・中性子（下）フルエンス 円柱中心（上）と端側（下）における 付与エネルギーの深さ分布 70