モンテカルロ法による粒子輸送計算
平山 英夫、波戸 芳仁
KEK, 高エネルギー加速器研究機構
2022.8.8

モンテカルロ法• 乱数を使用して問題を解く手法を “モンテカルロ法”という。
– “モンテカルロ”という名称は、J. von Neumann と S. M.
Ulamがつけたものである。

• 従って、乱数は、モンテカルロ法にとって最も重要なものである。

乱数の発生方法• サイコロやルーレット等を使用する– 非常に遅い• 乱数表を使用する– 統計的な性質が良く調べられている– 乱数表全体を計算機のデータとして持つ必要がある– 現時点では、非常に速い乱数発生法ではない• 放射性同位元素の崩壊の様な物理現象を利用した乱数を使用する– 数値化が容易でない、安定性と再現性の点で問題がある。
– https://www.kek.jp/ja/press/202207071401/

高エネルギー加速器研究機構（KEK）素粒子原子核研究所・理論センター量子基礎論グループの筒井泉特別准教授は、三菱電機株式会社情報技術総合研究所の鶴丸豊広主席技師長、東京大学大学院工学系研究科物理工学専攻の佐々木寿彦講師との共同研究において、パリティ（３次元空間の反転に対する）対称性を持つ原子核の崩壊現象から生成される量子乱数のランダム性と秘匿性に、厳密な証明を与えることに世界で初めて成功しました。

疑似乱数
••
••

最初にたね乱数, R0,を適当に選び、 Rn+1=f(Rn)の形の漸化式(普通は合同式)によって順次乱数を作り出す方法。
疑似乱数は、 mを法とする剰余である。
最大でもm個の異なった整数しか存在しない。従って、疑似乱数は有限の周期がある。
良い疑似乱数の条件––––•

速く発生できる
周期が長い
再現性がある
統計的特性が良い発生した疑似乱数をmで割ることにより、 0 から 1の疑似乱数を作ることができる

線形合同法による疑似乱数• D. H. Lehmer が提唱した線形合同法が最も広く使われている。 Rn+1=mod(aRn+b,m)– mod(aRn+b,m)は、 aRn+bをmで割ったときの余り• a, b 及び m は、正の整数で m は、コンパイラーで使用可能な整数の最大値である。
NameRANDU
SLAC RAN1SLAC RAN6

a65539b0m23169069
663608491

00231231

電卓を使った疑似乱数の発生• R0=3, a=5, b=0 , m=16として、10個の乱数を発生する。
• ある段階で同じ乱数が発生する事を確かめる。
• 上記の場合の周期は、いくつになるか ?
• R0を変えた場合、どうなるか調べる。

n

Rn0
1

R0=3

Rn*5
15

Rn+1=mod(Rn*a,m)
15/16=0・・・余り15 R1=15

23mod(Rn*a,m)： Rn*aをmで割ったときの余り

nRn

Rn*5

Rn+1=mod(Rn*a,m)0 R0=31 R1=151515/16=0・・・余り15 R1=152 R2=113 R3=74 R4=35555/16=3・・・余り73535/16=2・・・余り3 R4=375 75/16=4・・・余り11 R2=11R3=7

乱数を使ったp の計算• Table 1 (SLAC RAN6を使って作成したもの)の任意の場所から順番(左から、右に)に2つの乱数(x,h)を選ぶ。
• 次の条件を満たすペアの数を数える。
R = x 2 + h 2 £ 1.0
AREA= p /4
p= 4 x (条件を満たすペアの割合)

表１の乱数の使い方

Trial No123456
78910
A/10=0.8ξ

ηR

0.896
0.759
0.251
0.371
0.492
0.789
0.397
0.5760.517
0.909

0.618
0.690
0.094
0.148
0.519
0.567
0.179
0.341
0.583
0.380

1.088
1.026
0.268
0.399
0.715
0.972
0.435
0669
0.779
0.985

(A/10)*4=3.2

R≦1○○
○
○
○○
○
○
A=8

参考になるWebページ等• 乱数について– http://www.nikonet.or.jp/spring/sanae/MathTop
ic/montecalro/montecalro.htm– さいころの数字、乱数を使ったπの計算• 疑似乱数とモンテカルロ法– http://www.sm.rim.or.jp/~shishido/pie.html– 繰り返し数を増やすことによるπの変化• Newton 2009年8月号29ページ

他の方式の乱数• Marasaglia-Zaman 乱数– G. Masaglia and A. Zaman, “A New Class of Random NumberGenerator”, Annals of Applied Probability 1(1991)462-480.
– 周期が非常に長い –2144 ~1043– 乱数の制御が若干複雑– 32-bit の計算機ならば、どの計算機にも適用できる• RANLUX 乱数– F. James, “A Fortran implementation of the high-quality
pseudorandom number generators”, Comp. Phys Comm 79 (1994)
111-114.
– 周期は10 171– 1-231のシードにより、オーバーラップすることのない独立した乱数を発生することができる– egs5の乱数発生方法として採用

離散型確率過程のサンプリング• x1, x2,......,xnを確率p1, p2,......., pnを持つn個の独立背反物理事象とする。(例えば、光子の反応における、光電吸収、コンプトン散乱、電子対生成とそれぞれの発生確率)nF ( xn ) = å pi = 1

累積分布関数i =1• hを0から1の間の一様な乱数とすると、以下の条件の時、事象xiが起きることになる。
i -1

ij =1j =1

F ( xi -1 ) = å p j £ h < F ( xi ) = å p j

離散型サンプリングの導入（１）
例題） 光電効果：３０％、コンプトン散乱：５０％、対生成：２０％を乱数をもちいてサンプリングせよ。
0.60.50.4
0.30.20.10
光電効果 コンプトン散乱 対生成

離散型サンプリングの導入（２）
1.2
10.8

乱数
h0.6

コンプトン
散乱の起こる確率0.40.20

光電効果 コンプトン散乱

対生成「累積分布関数」 または「積み上げ計算」

離散型確率過程の例
光子の反応の種類の決定• 光子の反応における光電吸収、コンプトン散乱、電子対生成とそれぞれの発生確率をPphoto, PCompt , Ppairとする。 .
Pphoto +PCompt + Ppair =1•h £ Pphoto , の時は、光電吸収• Pphoto £ h < Pphoto + PCompt , の時は、コンプトン散乱•Pphoto + PCompt £ h , の時は、電子対生成

サンプリング法 (連続型確率過程)•

ある物理現象がxとx+dxで発生する確率をf(x)dxとする。[a£x£b]f(x)を確率密度関数(PDF)という。
b( ò f ( x)dx = 1)a•

累積分布関数(CDF:F(x))
xF ( x) = ò f ( xi )dxia•

hを0から1の間の一様な乱数とすると、
xh = F ( x) = ò f ( xi )dxiaと関係づける事ができる。上式からxはx = F -1 (h )l

この式を解析的に解くことが可能な場合には、xを直接計算して決めが可能であり、“直接サンプリング法”と呼ぶ。

F(x)1h0axbx

直接サンプリング法の例-飛行距離の決定• 1個の入射粒子が単位距離当たりに衝突する確率をStとする時、lと l+dlの間で最初に衝突が起きる確率p(l )dl = e - St l St dllh = P(l ) = ò p(l1 )dl1 = 1 - e -S lt0

1
l = - ln(1 - h ) = -l ln(1 - h )Stl: 飛行距離l: 平均自由行程(mean free path)h: 乱数 (1-h とhは等価)

l = -l ln(h )無限体系

光子
初期条件：エネルギー、位置、方向e0, x0, y0, z0, u0, v0, w0

反応点までの距離 l の決定

l=-ln(d)/µ移動後の座標x=x0+u0l, y=y0+v0l, z=z0+w0ll

初期条件：エネルギー、位置、方向

e0, x0, y0, z0, u0, v0, w0

反応の種類の決定
光電吸収：a, コンプトン散乱:b, 電子対生成：ｃd<=a/(a+b+c):光電吸収a/(a+b+c)< d<=(a+b)/(a+b+c):コンプトン散乱d>(a+b)/(a+b+c):電子対生成

x=x0+u0l, y=y0+v0l, z=z0+w0ll

初期条件：エネルギー、位置、方向

e0, x0, y0, z0, u0, v0, w0

電子生成粒子の追跡

個々の粒子のエネルギー、
方向の決定

光子
l

x=x0+u0l, y=y0+v0l, z=z0+w0l

初期条件：エネルギー、位置、方向

e0, x0, y0, z0, u0, v0, w0

ｄ>l：反応場所まで移動

d<=l：距離 d だけ移動
同じ物質の場合：反応場所までの距離=l-d
リージョン境界物質が異なる場合：改めて反応場所を決定

反応場所

x=x0+u0l, y=y0+v0l, z=z0+w0lld
境界までの直線距離(d)を計算
初期条件：エネルギー、位置、方向

e0, x0, y0, z0, u0, v0, w0

情報の採取粒子が移動：エネルギー付与
飛程長
リージョン境界

境界を横切る光電吸収：エネルギー付与
カットオフエネルギー以下
追跡中止(例：体系外)

境界までの直線距離(d)を計算
ld

初期条件：エネルギー、位置、方向

e0, x0, y0, z0, u0, v0, w0

電子や陽電子は、物質中で多数回の弾性散乱をするので、この弾性散乱を光子と同じように扱うことは難しい。

Condensed History Technique

反応点までの距離の決定
l=-ln(d)/S

反応場所までの距離を多くの細かいステップに分割し、各ステップにおける多数回の弾性散乱による実際の飛程、方向や位置の変化を“多重散乱モデルを使って評価する

llx 5x11xl2xl4l3

電子あるいは陽電子
初期条件：エネルギー、位置、方向e0, x0, y0, z0, u0, v0, w0xxl6

l7xl8

荷電粒子は、移動に伴い電離あるいは励起により、そのエネルギーの一部を失う
各ステップでのエネルギー付与
真の飛程 x 阻止能(dE/dx)

手計算による光子の輸送計算• Table 1の乱数（SLAC RAN6で求めたもの）– 任意の場所から始めて良いが、その後は順番に(左から右。右端に来たら、次の行へ)表の乱数を使用する– 計算結果の有効数字は３桁

表１の乱数の使い方

手計算による光子の輸送計算1 (Fig.1)•

図1のように、厚さ50 cmの物質Aがあるとする。
––––•0.5 MeVの光子が物質Aの左側から垂直に入射すると仮定する。
平均自由行程は20 cmとする。
光電吸収とコンプトン散乱の比は、 1:1とする。
コンプトン散乱後、光子はエネルギーも方向も変わらないとする。

例1– 最初の乱数:0.234 -- l=-20.0 x ln(0.234)=29.0– 29.0(cm)<50.0(cm)– 次の乱数:0.208 (<0.5) -- 光電吸収（終了）•

例2–––––

次の乱数:0.906 -- l=-20.0 x ln(0.906)=1.97
1.97(cm)<50.0(cm)次の乱数:0.716 (>0.5) -- コンプトン散乱
次の乱数：0.996 -- l=-20.0 x ln(0.996)=0.0802
0.0802(cm)<50.0-1.97(cm)Single layerNo.

d(cm) Random
numberl(cm)d>lExp.150.0Exp.2d£l0.23429.0*0.20850.00.906

1.97*0.716*48.03

0.996 0.0802*0.600

*

47.95

0.183

34.00.868*

13.950.351

20.9

*

Random Photo. Compt.
number*手計算による光子の輸送計算2 (Fig.2)• 図2のように、厚さ30cmの物質Aの後ろに厚さ20cmの物質Bがあるとする。
– 0.5MeVの光子が物質Aの左側から垂直に入射すると仮定する。
– 物質Aの平均自由行程、光電吸収とコンプトン比は先と同じとする。
– 物質Bの平均自由行程は3cmとする。
– 物質Bの光電吸収とコンプトン散乱の比は、3:1とする。
– 先と同様、コンプトン散乱後、光子はエネルギーも方向も変わらないとする。

計算例• 添付の乱数を使用した例– 最初の乱数:0.329 -- l=-20.0 x ln(0.329)=22.2– 22.2(cm)<30.0(cm)– 次の乱数0.612 (>0.5) -- コンプトン散乱– 次に乱数:0.234 --l=-20.0 x ln(0.234)=29.0– 29.0(cm)>30.0-22.2(cm)– AとBの境界まで移動 (30.0cm)– 次の乱数:0.281 --l=-3.0 x ln(0.281)=3.80– 3.80(cm)<20.0(cm)Medium Ad(cm)Randomnumber

l(cm)

d>l

30.0

0.329

22.2

*

7.8

0.234

29.0

No.
Exp.1

d£l Random Photo. Compt
number0.612**

Medium Bd(cm)Randomnumberl(cm)d>ld£lRandom Photo. Compt
number

20.00

0.281

3.80

*

0.906

*

16.200.716

1.00

*

0.996

*

15.200.600

1.53

*

0.183*

複雑だがより実際に近い光子の追跡• 10cmのアルミ板について、以下の仮定の下で光子の飛跡を追跡する。
• 入射光子のエネルギーは、0.5MeVとする。
• 光子のエネルギーに関係なくコンプトン散乱において、光子の散乱角は90度単位とし、散乱の確率は同じとする。
• 散乱後の光子のエネルギーは、次式で計算する。
E=E0
æ E ö
1 + ç 0 ÷(1 - cosq )è 0.511 ø

複雑だがより実際に近い光子の追跡• 散乱の方位角は、0度と180度が1:1の確率で起きるとする。(コンプトン散乱は、X-Z平面で生じる事になる。粒子の進行方向に対して、左側を0度とする。)• mfp 及び反応の分岐比(branching ratio)は、図4と5から読みとる。
• 光子のカットオフエネルギーは、0.05MeVとする。

二次粒子生成を伴わない電子の飛跡• 7図の様な厚さ 1 mm にアルミニウム平板を考える。
• 1.0 MeV の電子が左側から垂直に入射する。
• 多重散乱による移動距離の補正を無視する。
• d-線や制動輻射等二次粒子生成反応は起きないとする。
• 電子のステップ長は、電子のエネルギーに関係な
く 0.01 cm とする。

二次粒子生成を伴わない電子の飛跡• 電子のエネルギーに関係なく、多重散乱による方向の変化は、90°間隔で、同じ確率で生じる。
– 0°-- 1/3, 90°-- 1/3, 180°-- 1/3• 多重散乱後の方位角は、 0º 又は180º とし、発生確率は同じとする。0º を進行方向の左側、180ºを右側とする。
• 電子のエネルギーに関係なく、非弾性散乱によるエネルギー損失は、0.01cmあたり0.04 MeVとする。
• 電子のカットオフエネルギーを0.01 MeVとする。

ランダム・ウォークの例• Newton 2009年8月号42-43ページ– 不規則な動きを再現する「ランダム・ウォーク」
として、1次元、2次元、3次元の例が示されている。
– 移動距離は、常に同じで一歩とし、次元を変えるとどの様な動きになるかが例示されている。

改定記録• 乱数の発生及びπ計算結果の例追加 2007.8• 下記を小変更 2007.11.29– Page7,8に「15/16=0・・・余り15 」を追加
– Page 11のRを小数点以下5桁から3桁に変更– Page 31の光電効果の後に「終了」を追加• 下記小変更 2007.12.6– Page 35 一番下の行 6.20 cm→3.80 cm• 2重層の最初の厚さを30cmにしたことによる修正(34,36,37頁） 2008.2.8• 例(31,35頁)で、lがλになっていたので修正 2008.2.8• Newton 2009年8月号の記事を紹介 2009.7.3• 24ページの(a+c)を(a+b)に修正 2009.8.8• 14ページの不等式の添え字を1だけ減らす* 2013.1.17