Geant4 で新しい物理過程をつくる

Geant4 で新しい物理過程を入れたくて、OSSだしドキュメントも充実しているし論文も講習会資料もあるけど、どこを読めばいいのかわからない、または読んだけどわからないという人向け。

実装済みの物理過程はたくさんあるのですが、それでも未実装な物理過程もあります。そんなニッチな物理過程をシミュレーションしたいときは、自分で実装します。

僕が実装したのが離散過程でしたので、以下は離散過程を中心に述べます。

Geant4 の初歩的な知識と C++多態性の理解を前提とします。

ステップを理解する

ステップという概念があります。ステップの切り方と、切れ目で何が起こるか分かれば実装の大きな助けになります。

まず2017年初心者講習会のレクチャー資料「放射線シミュレーションの概要」の p.12 を見てください。

run --> event --> step --> track

括りの大きさでいうと、こんな階層構造ができています。

プログラム1回の実行 (run) が、粒子の生成から停止までの流れ (event) を任意の回数だけ計算します。Event 中では物理過程*1が何回か起こります。物理過程で区切られた区間をステップ (step) といいます。現在粒子がどういう状態にあるかはトラック (track) が保持しています。

物理過程で区切られた区間をステップとしたので、ステップの終点では必ず何かしら物理過程が起こります。これこそ実装したかった新しい物理過程です。このとき起こる現象を、仮想関数 G4VDiscreteProcess::PostStepDoIt() のオーバーライドとして実装します。

ステップの終点でおめあての物理過程が起こることはわかりましたが、ではステップの終点が起こる場所を決めるのは誰でしょうか。

ステップの長さを決める

ステップの終点、つまり物理過程の発生点までの距離を決めるのは、サイコロです。サイコロの出目の決まり方は「放射線シミュレーションの概要」の p.17 や Appendix を見てください。累積確率がサイコロの出目にまで高まる距離がステップの終点です。ここで重要なのは、これで決まるのは Number of Mean Free Path (NMFP) という無次元量であって、物理的な距離ではないということです。定義として、 NMFP は発生点までの物理的な距離を平均自由行程で割ったものです。平均自由行程数とでもいいましょうか。 *2

NMFP は平均自由行程の何倍かという無次元量なので、物質によりません。仮に粒子がサイコロをふって、 NMFP = 2 だけ飛ぶと決めたとします。この粒子の飛跡が 2 つの異なる物質をまたいだとし、片方の物質で NMFP = 1 だけ、もう一方で NMFP = 1 だけ飛んだとしたら、それぞれの物質中で、その物質における平均自由行程ぶんだけ飛んだことになります。

NMFP で次に起こる相互作用を決める

NMFP を導入するのは計算量の抑制が目的だろうと思います。NMFP なしだと複数の相互作用から次に起こるものを決めるのが大変です。

まず、ある 1 種類だけの相互作用だけが起きると仮定します。 NMFP は一旦忘れます。現実で粒子が飛ぶときには、物質の性質から散乱断面積が決まり、平均自由行程が決まります。粒子が飛んでいく方向に物質が複数あっても、累積確率が十分に高くなるまで逐一平均自由行程を取得していけばステップの長さはいずれわかります。

ここで、 1 種類だけ相互作用が起こるという仮定を外し、任意の種類の相互作用が起こるとします。すると、次に起こる相互作用がどれであるかを決める必要が出てきます。そのために想定しうるすべての相互作用に対してステップの長さを計算し比較する必要があります。

しかし実際に起こる相互作用はステップが最短の 1 つだけで、しかもその相互作用が起こると粒子が飛んでいく方向が変わります。そうなるといま起きたもの以外の相互作用に関するステップ長の計算が無駄です。この無駄は各相互作用について物質とは無関係に NMFP をサイコロで決め、 NMFP が最小であるものを選択することでなくすことができます。

NMFP の比較で次のステップの物理過程が決まることはわかりました。ところで、1 ステップが複数の物質をまたぎそうなときはどうしたらよいでしょうか?この問題に対処するために、物質境界の通過もステップの終了点としましょう。すなわち、候補となる物理過程の NMFP に平均自由行程をかけた実際の空間における距離と、物質の境界までの距離を比較して、境界までの距離のほうが短かったらそこで 1 ステップを終えます。境界をまたいだあとは、また同様の比較を行います。これを繰り返すことで粒子の生成から停止までを追うことができます。

なお、NMFP に平均自由行程をかけて実際の空間における距離に換算する際、物理過程ごとに物質の情報から平均自由行程を計算しますが、この計算は G4VDiscreteProcess::GetMeanFreePath() をオーバーライドしたメソッドとして実装します。

ステップ開始から終了までの流れ

まとめましょう。いま粒子が新たにステップを始めようとしています。

  1. 想定される物理過程すべてについて、 NMFP がサンプリングされます。
  2. NMFP に平均自由行程をかけた実際の空間での距離と、次の物質までの距離のうち最短の過程が次のステップを決めます。
  3. 粒子が輸送されます。このときの輸送距離はすべての候補の残りの距離から差し引かれます。
  4. ステップの輸送が終わると、そのステップを決定した物理過程が起こります。
  5. 次のステップを決めるために、いま起こった過程だけ NMFP が再サンプルされます。すべての過程について、2 に戻って次のステップに行きます。

メソッドが呼び出されるタイミング

ステップが繰り返される流れは以上ですが、各メソッドがいつ呼び出されるかをまとめましょう。

  • PostStepDoIt(): ステップ終了点で物理過程を起こす。
  • GetMeanFreePath(): 次のステップを決めるとき。物質境界をまたいだ直後もこれに含まれる。

これら 2 つは自分で実装する必要があるものです。ついでに次のメソッドを知っておくと理解の助けになると思います。

具体的実装

具体的にどう実装するかについては、既存のプロセスの実装を真似るのが計算量的に安全かと思います。ツールキット開発者向けガイドのハドロンの物理過程のセクションなどが参考になるでしょう。ソースコードを読む際は doxygen は参照箇所がわかって便利ですし、 GitLab の Geant4 リポジトリは検索ができて便利です。(例: 可視光光子のレイリー散乱の実装*3

新しい物理過程をプログラムに組み込む

新しく作った物理過程をプログラムに組み込む方法は2017年初心者講習会の資料や各種ユーザーガイドを参照してください。通常どおり Physics List に追加するだけです。

*1:一般的な言い方では相互作用

*2:平均自由行程を相互作用長 interaction length という分野もあるようです。実際、Geant4 コードでも使われています。 cf. G4VDiscreteProcess.cc

*3:Geant4 は多くの研究者の貢献によって成り立っており、ドキュメントのメンテナンスも研究者らの好意によって成り立っています。それを鑑みると最良のリファレンスはソースコードだろうと思います。