計算効率を上げるためのvariance reduction(粒子のウェイト)の利用 A

Slides:



Advertisements
Similar presentations
P HI T S 陽子ビームで雪だるまを溶かそう Multi-Purpose Particle and Heavy Ion Transport code System title 年 3 月改訂.
Advertisements

宇宙線ミューオンの測定 久野研究室 4回生 卒業研究 荒木 慎也 宮本 紀之 室井 章. 目次 実験内容 測定方法・結果 ・検出装置とセットアップ 解析 ・バックグラウンド除去 ・検出効率 ・立体角 ・文献 値との比較 まとめ.
放射線の計算や測定における統計誤 差 「平均の誤差」とその応用( 1H) 2 項分布、ポアソン分布、ガウス分布 ( 1H ) 最小二乗法( 1H )
P HI T S α 線, β 線, γ 線,中性子線を止めるに は? Multi-Purpose Particle and Heavy Ion Transport code System title 年 5 月改訂.
P HI T S PSFC4PHITS の使い方 Multi-Purpose Particle and Heavy Ion Transport code System Title 年 2 月改訂.
計算効率を上げるためのvariance reduction(粒子のウェイト)の利用
PHITS 講習会 基礎実習(III): 計算条件の設定
Multi-Purpose Particle and Heavy Ion Transport code System
Multi-Purpose Particle and Heavy Ion Transport code System
PHITS講習会 基礎実習(I): 体系及び線源の定義
PHITS講習会 基礎実習(II): Tally(タリー)の定義
平成20年度 核融合科学研究所共同研究 研究会 「負イオン生成および負イオンビーム加速とその応用」 プロセスプラズマのPIC計算のモデリング
強度変調回転放射線治療(VMAT)の シミュレーション
Multi-Purpose Particle and Heavy Ion Transport code System
Multi-Purpose Particle and Heavy Ion Transport code System
Multi-Purpose Particle and Heavy Ion Transport code System
物理化学(メニュー) 0-1. 有効数字 0-2. 物理量と単位 0-3. 原子と原子量 0-4. 元素の周期表 0-5.
Multi-Purpose Particle and Heavy Ion Transport code System
放射能計算の基礎と DCHAIN-SPの応用実習
PHITSによるX線治療シミュレーション 基本操作復習編
ParaViewを用いたPHITS 計算結果の3次元表示
放射線の計算や測定における統計誤差 「平均の誤差」とその応用(1H) 2項分布、ポアソン分布、ガウス分布(1H) 最小二乗法(1H)
Determination of the number of light neutrino species
相対論的重イオン衝突実験PHENIX におけるシミュレーションによる charm粒子測定の可能性を探る
中性子過剰核での N = 8 魔法数の破れと一粒子描像
埼玉大学大学院理工学研究科 物理機能系専攻 物理学コース 06MP111 吉竹 利織
PHITS 講習会 基礎実習(III): 計算条件の設定
natMg+86Krの反応による 生成核からのβ線の測定と GEANTによるシミュレーションとの比較
Multi-Purpose Particle and Heavy Ion Transport code System
誘導放射能の評価方法: PHITSとDCHAIN-SPの接続計算
α線,β線,γ線,中性子線を止めるには?
京都大学大学院医学研究科 画像応用治療学・放射線腫瘍学 石原 佳知
Multi-Purpose Particle and Heavy Ion Transport code System
For the PHENIX collaboration
応用実習用資料 Environmental radioactivity
鉄核媒質中の閾値近傍における 中性π中間子生成実験
Multi-Purpose Particle and Heavy Ion Transport code System
MeV internal meeting Oct. 2, 2015
計算効率を上げるためのvariance reduction(粒子のウェイト)の利用 A
論文講読 Measurement of Neutrino Oscillations with the MINOS Detectors in the NuMI Beam 2009/11/17 Zenmei Suzuki.
PHITS 講習会 基礎実習(III): 計算条件の設定
応用実習用資料 Neutron target
IAEA phase space fileを用いた X線治療シミュレーション
理研RIBFにおける 中性子過剰Ne同位体の核半径に関する研究
Multi-Purpose Particle and Heavy Ion Transport code System
Multi-Purpose Particle and Heavy Ion Transport code System
治療用フィルムによる線量分布測定の 基礎的検討Ⅱ
Multi-Purpose Particle and Heavy Ion Transport code System
EGSに対応した粒子軌跡と 計算体系の3次元表示ソフト - CGVIEW -
Charmonium Production in Pb-Pb Interactions at 158 GeV/c per Nucleon
Multi-Purpose Particle and Heavy Ion Transport code System
B物理ゼミ Particle Detectors:Claus Grupen, Boris Shwartz (Particle id
電子後方散乱の モンテカルロ計算と実験の比較 総研大 桐原 陽一 KEK 波戸 芳仁、平山 英夫、岩瀬 広.
Dark Matter Search with μTPC(powerd by μPIC)
Multi-Purpose Particle and Heavy Ion Transport code System
Multi-Purpose Particle and Heavy Ion Transport code System
NaIシンチレーターを使った 放射線検出システムの開発
SciFi を用いたΣ+p散乱実験での (ほろ苦い)思い出
天体核実験用の窓無しガス標的と ガス循環系の開発
卒業論文発表 中性子ハロー核14Beの分解反応 物理学科4年 中村研究室所属   小原雅子.
京大理 身内賢太朗 平成22年度東京大学宇宙線研究所 共同利用研究成果発表会
ユーザーコードに記述する事項の概要 2010年7月21日 KEK 波戸.
α線,β線,γ線,中性子線を止めるには?
計算と実測値の比較 高エネルギー加速器研究機構 平山 英夫.
Multi-Purpose Particle and Heavy Ion Transport code System
α decay of nucleus and Gamow penetration factor ~原子核のα崩壊とGamowの透過因子~
大型ヘリカル装置における実座標を用いた 粒子軌道追跡モンテカルロコードの開発
計算効率を上げるためのvariance reduction(粒子のウェイト)の利用 B
Recoil catcher法による質量数90領域の
荷電粒子の物質中でのエネルギー損失と飛程
Presentation transcript:

計算効率を上げるためのvariance reduction(粒子のウェイト)の利用 A PHITS Multi-Purpose Particle and Heavy Ion Transport code System 計算効率を上げるためのvariance reduction(粒子のウェイト)の利用 A PHITS講習会 入門実習 2018年8月改訂

実習内容 1.はじめに 2.中性子深層透過計算 [importance]  Geometry splitting and Russian Roulette (セルインポータンス法) 3.薄膜からの生成粒子計算 [forced collisions]  強制衝突法

中性子深層透過計算 厚い遮蔽体の中性子輸送を計算し、深さ毎の実効線量率分布を計算 まずはimp.inpを走らせて,計算時間を計ってみよう [ P a r a m e t e r s ] icntl = 0 maxcas = 8000 maxbch = 1 ... [ T - T r a c k ] (3つ目) mesh = xyz … axis = xz z-txt = Effective dose [pSv/source] file = imp-dose-xz.out part = neutron epsout = 1 gshow = 1 multiplier = all mat mset1 all ( 1.0 -202 ) ヒストリ数は8,000 計算時間はphits.outの最後の方(total CPU time)に出力される 半径50cm、深さ180cmのコンクリート円柱 空気 コンクリート 14MeV 中性子 実効線量は、[t-track]で計算した中性子フルエンスに線量換算係数を乗じて導出

実効線量の計算 [t-track]で計算したフルエンスは、格納もしくは[multiplier]セクションで定義した複数の関数と掛け合わせることができる /phits/data/multiplier/フォルダに 格納されている線量換算係数データ一覧の抜粋 [ T - T r a c k ] (2,3つ目) … multiplier = all mat mset1 all ( 1.0 -202 ) ID番号 線量換算係数データ pSv・cm2 -200 H*(10) -201 ICRP60定義に基づく実効線量(前方照射) -202 ICRP103定義に基づく実効線量(前方照射) -203 ICRP103定義に基づく実効線量(等方照射) (C k): C 規格化定数, k ID番号 線源強度 線量換算係数 フルエンスを線源強度で規格化したとき、(pSv)単位で線量を導出 今回は規格化定数Cに線源強度を代入しないため、 フルエンス(1/cm2/source) x 線量換算係数(pSv・cm2) = 線源1個あたりの線量(pSv/source) [multiplier]セクションについてはマニュアルを参照 4

中性子深層透過計算 シングルコア( Intel Core i5-4590 CPU @ 3.30GHz)の計算 効率良く計算したい。 ヒストリー数 8,000 total cpu time = 6.39 sec imp-dose-xz.eps ヒストリー数 500,000 total cpu time = 215.36 sec imp-dose-xz.eps 効率良く計算したい。 分散低減法 ウェイトの概念を利用

モンテカルロ計算におけるWeightの概念 例:track length タリー   Li: i番目の粒子の軌跡長さ Wi: i番目の粒子のウェイト n0: 総ヒストリー数 ウェイト: モンテカルロ計算における粒子の重要度、通常の計算は1*。 ウェイトを操作して、統計上発生しにくいイベントを人為的に発生させる。 ただし、ウェイトを操作した場合、ヒストリー毎の頻度分布([t-deposit], output = depositなど)の計算ができなくなるので注意が必要。 *ただし、核データを使用した中性子輸送計算については、e-modeを使用しない場合にウェイトが変化する可能性があります。

実習内容 1.はじめに 2.中性子深層透過計算 [importance]  Geometry splitting and Russian Roulette (セルインポータンス法) 3.薄膜からの生成粒子計算 [forced collisions]  強制衝突法

中性子がセル1とセル2を横切る時、I2/I1を設定 セルインポータンス法 領域(セル)にインポータンスIを設定 中性子がセル1とセル2を横切る時、I2/I1を設定 I2/I1 > 1: Particle Split (粒子の分割) 例 I2/I1 が整数の場合  I2/I1 = 3: 3つの粒子にsplit(分割) I2/I1 が非整数の場合 I1=1 I2=3 I2/I1 = 2.75: 75%の確率で3つの粒子に分割 25%の確率で2つの粒子に分割 W=1/3 W=1/3 W=1/3 W=1 I1=1 I2=1/3 W=3 I2/I1 < 1: Russian Roulette (消滅か生存) 1/3の確率 W=1 消滅 2/3の確率

1ヒストリーのときの中性子飛跡を確認 imp-hist.inp imp-hist.inpを走らせて,中性子の飛跡を確認してみよう 空気 [importance] part = neutron reg imp 10 1 1 1 2 1 imp-hist.inp imp-hist.inpを走らせて,中性子の飛跡を確認してみよう 空気 コンクリート (半径20cm, 深さ8cm×2) 5MeV 中性子 I10=1.0 I1=1.0 I2=1.0 標準では、インポータンスは全て1 imp-trajectory.eps

インポータンスを大きくしたときの中性子挙動を確認 [importance] part = neutron reg imp 10 1 1 2 2 4 imp-hist.inp ここで反応が起きた I10=1 I1=2 I2=4 この時点で既に2つに分かれている imp-trajectory.eps

インポータンスを小さくしたときの中性子挙動を確認 [importance] part = neutron reg imp 10 1 1 1/2 2 1/4 imp-hist.inp 1/2の確率で生存 1/2の確率で消滅 I10=1 I1=1/2 I2=1/4 imp-trajectory.eps

インポータンスを極端に大きくしたときの中性子挙動を確認 [importance] part = neutron reg imp 10 1 1 2 2 100 粒子が分割しすぎて計算時間が無駄になってしまう (統計はあまり良くならない) imp-hist.inp I10=1 I1=2 I2=100 imp-trajectory.eps セル間のインポータンスの比は最大2~3程度が良い “A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380

インポータンスを使用した中性子深層透過計算例 Ii+1/Ii = 2.5 と設定 imp.inpの [importance]セクションを有効にして、実行してみよう [importance] part = neutron reg imp 1 2.5**0 2 2.5**1 3 2.5**2 4 2.5**3 5 2.5**4 6 2.5**5 7 2.5**6 8 2.5**7 9 2.5**8 10 2.5**9 11 2.5**10 12 2.5**11 半径50cm、深さ180cmのコンクリート円柱 セルの厚さは15cm 半径1cmの14MeV中性子をz=0からz軸に沿って入射 14MeV 中性子

インポータンスを使用した中性子深層透過計算例 3.30GHz, single ヒストリー数 8,000 [importance] 未使用 total cpu time = 6.39 sec Dose Relative error imp-dose-xz_err.eps imp-dose-xz.eps [importance] 使用 total cpu time = 48.39 sec 1 相対誤差の確認重要(赤1.0、黄 約0.1、緑 約0.01 )

インポータンスを使用した中性子深層透過計算例 With importance Without importance 元データ:imp-dose-reg.out 領域(15cm厚)毎の線量変化 14

インポータンス利用の注意点 セル間の大きなインポータンスの比は良くない。多くの粒子分割を一度に引き起こす。 良い例 粒子 1 2 4 8 16 32 悪い例 粒子 1 1 8 8 8 32 セル間のインポータンスの比は最大2~3程度が良い 参考文献 “A Sample Problem for Variance Reduction in MCNP” LA-10363-MS DE86 004380

インポータンスの誤った使用例 先ほどの結果と比較して相対誤差が大きい 悪い例: 大きなインポータンスのギャップを作って実行してみよう。 Cell 5と6の間に急激なImportanceギャップ 大きなインポータンスのギャップを作って実行してみよう。 [importance] part = neutron reg imp 1 2.5**0 2 2.5**0 3 2.5**0 4 2.5**0 5 2.5**0 6 2.5**5 7 2.5**6 8 2.5**7 9 2.5**8 10 2.5**9 11 2.5**10 12 2.5**11 imp-dose-xz.eps imp-dose-xz_err.eps 先ほどの結果と比較して相対誤差が大きい

実習内容 1.はじめに 2.中性子深層透過計算 [importance]  Geometry splitting and Russian Roulette (セルインポータンス法) 3.薄膜からの生成粒子計算 [forced collisions]  強制衝突法

強制衝突とは? 放射線を衝突の確率の低いターゲット(薄膜)に強制的に衝突させる。 応用例: 薄膜を用いた断面積測定実験を模擬するときに利用 10μm厚さSi 10μm厚さSi 100MeV陽子 強制衝突なし 陽子は薄膜をそのまま通過。 強制衝突あり 陽子は薄膜と必ず衝突。 応用例: 薄膜を用いた断面積測定実験を模擬するときに利用

薄膜を用いた断面積測定実験を模擬しよう force.inpを実行してみよう [t-product]: ターゲット内で生成するα粒子とSiの生成エネルギー分布 → 真の断面積 [t-cross]:検出器に入射したα粒子とSiのエネルギー分布 → 測定上の断面積 force.inp 強制衝突無し:殆どの入射粒子は衝突しない force.inpを実行してみよう maxcas = 5000 maxbch = 5 product.eps 半径0.5cm, 10μm厚さSi 100MeV陽子 真空 検出器 cross.eps

強制衝突法 薄膜から生成する粒子を解析するときに利用 二つに分割 非衝突粒子 非衝突のウェイト:Wi×exp(-Sd) d:セルを横切る距離 強制衝突領域 S: 巨視的断面積 衝突位置は断面積に従って確率的に決定 強制衝突法はヒストリー分散を減らすが、ヒストリー当たりの計算時間が増える。 セルでは一つ以上の衝突が起こる。 fcl:強制衝突コントロールパラメータ  (通常 fcl = -1 を使用) fcl = -1:強制衝突による生成粒子は通常の衝突をさせる(ウェイトカットオフは考慮しない) fcl = 1:強制衝突による生成粒子は,ウェイトカットオフになるまで強制衝突させる

強制衝突させてみよう 強制衝突の設定を行おう ターゲット内での粒子生成時点の情報をタリー force.inp [Forced Collisions] off part = proton reg fcl 1 -1.0 offを削除して実行 product.eps 強制衝突による生成粒子のウェイトが、ウェイトカットオフwc2を下回る。 ⇒ Russian Rouletteが働くため、粒子輸送されにくい。 α track-xz.eps cross.eps

強制衝突で発生させた2次粒子を輸送するには? 強制衝突で生成する粒子を輸送させるため、 ウェイトカットオフwc2を極力低く設定してみよう。 force.inp [parameters] ………. wc2(1) = 1.0E-12 # weight cutoff of proton wc2(18) = 1.0E-12 # weight cutoff of Alpha wc2(19) = 1.0E-12 # weight cutoff of Nucleus product.eps α Si track-xz.eps cross.eps 多くの生成Siがターゲット内で止まるため,この実験からは断面積が正しく測定できない

まとめ 統計がたまらない計算では、粒子のウェイト(重要度)を調整することにより計算時間を短縮することが可能 深層透過計算には、セルインポータンス法が有効 薄膜における生成粒子情報を計算するにはforced collisionsが有効 タリー領域が小さすぎて統計がたまらない場合はポイントタリーが有効 *ポイントタリー[t-point]の使い方はutility\tpointを参照 *隣り合うセル間のインポータンス値の比は2~3以下が望ましい *複雑な体形でインポータンスの設定が難しい場合は、Weight Window Generatorを使うのが便利(weight Bで説明)