等方的ビーム
- 乱数を使って等方的ビームを打つためのコードを覚え書く。途中まで書いたところで、同じようなことを書いているページを見つけてしまった。書いたところは残しておく。さらにうまい方法が Wikipedia に載ってた。
3 次元
等方的ビームを打つためは、単位球面上に一様に分布する (x, y, z) の組を作れば良い。球面座標系では (MathWorld - Spherical Coordinates)
であり (theta と phi をページ内の注意を参考に逆にした)、ヤコビアンは
- である。よって、単位球面上の微小面積要素は
と書ける。これより、
に対しては
の重みを付け、
に対しては重みを付けないで乱数を振れば、単位球面上に一様に分布する (x, y, z) の組を作ることができる。
に従う確率密度分布を作るには、これの積分の逆関数を使って一様乱数を振れば良い (参考: 乱数とモンテカルロ法)。すなわち、
を使って、一様な分布をもつ確率変数
の変換を行う。具体的には、以下のような変換を行えば良い。ここで、
はそれぞれ
の範囲で一様に分布する確率変数。
ここで、
を
に代入すれば、
(Inverse Trigonometric Functions) より
となる。 ROOT のマクロの実装 (isotropic_3d.C) を用いて結果をプロットすると以下のような図になる。
左上の図が、xyz 空間上に (x, y, z) の点をプロットしたもの。ぐりぐり回転すれば分かるが、単位球面上に一様に分布している。真ん中上の図では
の分布が正しく
に従っている事が分かる。左上の図では、
が一様に分布している事が分かる。x,y,z の分布は、中段に表してあり、一様分布となる。なんでこうなるのかは後で考える。下段には、xy, xz, yz の相関を表してあり、動径方向では、外に行くほど密度が高くなるが、角度方向を見れば一様に分布していることが分かる。
4 次元
4 次元の場合も 3 次元の場合とほぼ同じとなる。MathWorld - Hypersphere より、4 次元の場合の球面座標系は
- となる (theta と psi を逆にした)。これより、ヤコビアンは
- となる。よって、単位超球面上の微小面積(体積?)要素は
となる (計算ノート)。これより
に対しては
の重みを付け、
に対しては
の重みを付け、
に対しては重みを付けないで乱数を振れば、単位超球面上に一様に分布する
の組を作ることができる。
に従う確率密度分布を得るには、これを積分し逆関数を得ればよい。ただし、
であるので (参考: I. S. Gradshteyn and M. Ryzhik, Tables of Integrals, Series and Products (Academic Press, New York, 2007), 7th ed., p. 153. Sect. 2.513 5)、逆関数は解析関数にならないが、数値的には解ける。ここで
と置き、この関数の逆関数を
と置けば、単位超球面上に一様に分布する
の組を得る式は以下のようになる。
ここで
は
の範囲で一様に分布する確率変数。これを ROOT で実装すると isotropic_4d.C のようなマクロとなる。ここで、 逆関数
の値を得るために、TF1 クラスの TF1::GetX メソッドを用いた。マクロの結果は以下のようになる。
左上の図が、4 次元の相関を取った図で、x, y, z 軸が x2, x3, x4, に相当し、色が x1に相当する。なんだかよくわからないが、動径方向の内側は赤で 1 に近くなり、外側に行くほど青く、0 に近くなる。よって、超球面上に分布していそうに思える。また、動径方向に垂直な方向に一様に分布している。隣の 4 つの図は、 x1x2x3, x1x2x4, x1x3x4, x2x3x4, の相関を 3 次元でプロットしたもの。 動径方向の外に行くほど密度が濃くなる。また、動径方向に垂直な方向については一様に分布していることが分かる。中段は、各変数の分布。
が期待通りの分布
になっている。x1, x2, x3, x4, の分布はすべて同じ形となり、それぞれの分布では横軸の 0 を中心に左右対称となっているので、正しそう。なぜ円形の分布になるのかは後で考える。下段の図は 2 変数間の相関を表しており、すべて円内に一様に分布する。すべての図で、角度方向に一様に分布しているので、x1, x2, x3, x4, は一様に単位超球面上に分布しているように思える。