GATE10の使い方 ２ ~物理モデルの比較~

GATEはバックグラウンドでGeant4を動かしています。

Geant4には様々な物理モデルが搭載されていおり、目的に応じて最適なものを選択する必要があります。これは、すべての種類の放射線、広いエネルギー範囲の物理現象を説明する一意のモデルが存在しないためです。

Geant4で使用できる物理モデルはここに示されています。

今回はそれを順に実行して結果を比較していきます。

まず最初にコードのフローチャートを示します。親プロセスが最初に呼ばれ、その中で子プロセスが何度も呼ばれる（物理モデルの数だけ繰り返す）という流れです。

次にコードの全体像を示します。コメントやコードの整形はChatGPTにお願いしました。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

"""
Theme01:
- 指定したエネルギーのガンマ線を水ファントムに照射する
- 統計情報（時間・イベント数・トラック数・ステップ数など）を取得する
- 物理リスト（特にEMオプション）の違いによる計算結果・速度の違いを比較する

このスクリプトは2つの動き方があります。

(1) 親プロセス（まとめ役）
    - 物理リストの候補を順番に「別プロセス」で実行する
    - それぞれの結果を JSONL ファイルに追記していく
    - 最後に一覧表を表示する

(2) 子プロセス（1つの物理リストを実行する役）
    - 環境変数 GATE_PHYSICS_LIST に指定された物理リストでシミュレーションを実行
    - 結果を JSONL に1行追記する（親プロセスが回収する）
"""

import os
import sys
import json
import subprocess
from pathlib import Path

import opengate as gate
from opengate.utility import g4_units

# =============================================================================
# 単位の定義（Geant4の単位系）
# =============================================================================

m = g4_units.m
cm = g4_units.cm
mm = g4_units.mm
keV = g4_units.keV
Bq = g4_units.Bq
sec = g4_units.second  # ※今のコードでは使っていないが、読みやすさのため残してある

# =============================================================================
# 設定（ユーザーが変えやすい場所）
# =============================================================================

CONFIG = {
    # 線源
    "source_particle": "gamma",
    "source_energy_keV": 511.0,
    "source_n": 200000,

    # ファントム（水箱）
    "phantom_material": "G4_WATER",
    "phantom_size_cm": [40.0, 40.0, 40.0],          # [X, Y, Z] cm
    "phantom_translation_cm": [0.0, 0.0, 25.0],     # [X, Y, Z] cm

    # 比較する物理リスト
    "physics_lists": [
        "QGSP_BERT_EMV",
        "QGSP_BERT",
        "QGSP_BERT_EMZ",
        "QGSP_BIC",
        "QGSP_BIC_EMY",
        "QGSP_INCLXX",
        "FTFP_BERT",
        "FTFP_BERT_EMV",
        "FTFP_BERT_EMZ",
    ],
}

# =============================================================================
# 1) シミュレーションを「作る」関数
# =============================================================================

def create_simulation(physics_list_name: str):
    """
    指定された物理リスト名で Simulation を組み立てる。
    ここではまだ実行しない（runしない）。
    """

    # Simulation オブジェクトを作成
    sim = gate.Simulation()

    # -------------------------------------------------------------------------
    # World（最上位ボリューム）設定
    # -------------------------------------------------------------------------
    world = sim.world
    world.size = 
    world.material = "G4_AIR"

    # -------------------------------------------------------------------------
    # 水ファントム（箱）を作成して World の中に置く
    # -------------------------------------------------------------------------
    waterbox = sim.add_volume("Box", "Waterbox")
    waterbox.size = [value_cm * cm for value_cm in CONFIG["phantom_size_cm"]]
    waterbox.translation = [value_cm * cm for value_cm in CONFIG["phantom_translation_cm"]]
    waterbox.material = CONFIG["phantom_material"]
    waterbox.color = [0, 0, 1, 0.3]  # 青・半透明（可視化する場合に見やすい）

    # -------------------------------------------------------------------------
    # 線源（ガンマ線）を作成
    # -------------------------------------------------------------------------
    source = sim.add_source("GenericSource", "gamma_source")
    source.particle = CONFIG["source_particle"]
    source.energy.mono = CONFIG["source_energy_keV"] * keV

    # 位置：点線源、ファントム前方10 cm（Z方向マイナス側）
    source.position.type = "point"
    source.position.translation = [0, 0, -10 * cm]

    # 向き：Z軸正方向へ打つ
    source.direction.type = "momentum"
    source.direction.momentum = [0, 0, 1]

    # 発生させる粒子数
    source.n = CONFIG["source_n"]

    # -------------------------------------------------------------------------
    # 物理リストを設定
    # -------------------------------------------------------------------------
    sim.physics_manager.physics_list_name = physics_list_name

    # Production cut（生成カット）設定：例として0.1 mm
    sim.physics_manager.global_production_cuts.gamma = 0.1 * mm
    sim.physics_manager.global_production_cuts.electron = 0.1 * mm
    sim.physics_manager.global_production_cuts.positron = 0.1 * mm
    sim.physics_manager.global_production_cuts.proton = 0.1 * mm

    # -------------------------------------------------------------------------
    # 統計情報アクターを追加
    # -------------------------------------------------------------------------
    stats_actor = sim.add_actor("SimulationStatisticsActor", "Stats")
    stats_actor.track_types_flag = True  # 粒子種別ごとのトラック数も記録する

    return sim, waterbox, source, stats_actor

# =============================================================================
# 2) シミュレーションを「実行して結果を表示・保存する」関数
# =============================================================================

def run_simulation_and_print_report(physics_list_name: str):
    """
    create_simulation() で組み立てたシミュレーションを実行し、
    統計情報を表示し、必要なら JSONL ファイルに結果を追記する。
    """

    # まずシミュレーションを構築（ここで失敗することがあるのでtry）
    try:
        sim, waterbox, source, _stats_actor = create_simulation(physics_list_name)
    except Exception as exc:
        print("=" * 80)
        print("テーマ1: 物理リスト比較")
        print("=" * 80)
        print(f"物理リスト: {physics_list_name}")
        print(f"物理リストの初期化に失敗しました: {exc}")
        print("=" * 80 + "\n")
        return

    # -------------------------------------------------------------------------
    # 実行前の情報表示（初心者が条件を追えるように）
    # -------------------------------------------------------------------------
    print("=" * 80)
    print("テーマ1: 物理リスト比較")
    print("=" * 80)
    print(f"物理リスト: {physics_list_name}")
    print(f"粒子種類: {source.particle}")
    print(f"エネルギー: {source.energy.mono / keV:.1f} keV")
    print(f"粒子数: {source.n:,}")
    print(f"ファントム材質: {waterbox.material}")
    print(
        f"ファントムサイズ: {waterbox.size[0] / cm:.1f} x "
        f"{waterbox.size / cm:.1f} x {waterbox.size / cm:.1f} cm"
    )
    print("=" * 80)
    print("\nシミュレーション開始...")

    # -------------------------------------------------------------------------
    # 実行
    # -------------------------------------------------------------------------
    sim.run()

    print("\nシミュレーション完了！")
    print("=" * 80)

    # -------------------------------------------------------------------------
    # 統計情報の取得と表示
    # -------------------------------------------------------------------------
    output = sim.get_actor("Stats")  # "Stats" という名前のアクターの結果を取り出す

    print("\n【統計情報】")
    print(output)

    # duration は ns 単位なので秒に変換（ns -> s）
    duration_s = output.counts.duration * 1e-9

    events = output.counts.events
    tracks = output.counts.tracks
    steps = output.counts.steps

    # 1秒あたりイベント数（PPS）
    if duration_s > 0:
        pps = events / duration_s
    else:
        pps = 0.0

    print("\n【詳細統計】")
    print(f"総実行時間: {duration_s:.2f} 秒")
    print(f"総イベント数: {events:,}")
    print(f"総トラック数: {tracks:,}")
    print(f"総ステップ数: {steps:,}")
    print(f"1秒あたりの粒子数 (PPS): {pps:.2e}")

    # 粒子種別トラック数（存在するときだけ表示）
    if output.counts.track_types:
        print("\n【粒子種別トラック数】")
        for particle, count in sorted(output.counts.track_types.items()):
            print(f"  {particle:15s}: {count:,}")

    print("\n" + "=" * 80)

    # -------------------------------------------------------------------------
    # 結果を JSONL 形式で追記（親プロセスが後でまとめる）
    # -------------------------------------------------------------------------
    result = {
        "physics_list": physics_list_name,
        "duration_s": duration_s,
        "events": events,
        "tracks": tracks,
        "steps": steps,
        "pps": pps,
    }

    results_file = os.environ.get("GATE_RESULTS_FILE")
    if results_file:
        with open(results_file, "a", encoding="utf-8") as f:
            f.write(json.dumps(result) + "\n")

# =============================================================================
# 3) 物理リストごとに「別プロセス」で実行する関数
# =============================================================================

def run_one_physics_list_in_subprocess(physics_list_name: str):
    """
    物理リスト1つを「別プロセス」で実行する。
    - 子プロセス側は環境変数 GATE_PHYSICS_LIST を見て単発実行する
    - こうすることで、物理リスト初期化が互いに影響しにくくなる
    """

    env = os.environ.copy()
    env["GATE_PHYSICS_LIST"] = physics_list_name

    # 結果ファイル名も子プロセスに引き継ぐ
    if "GATE_RESULTS_FILE" in os.environ:
        env["GATE_RESULTS_FILE"] = os.environ["GATE_RESULTS_FILE"]

    # このスクリプト自身を「pythonで」実行する
    cmd = [sys.executable, os.path.abspath(__file__)]
    subprocess.run(cmd, env=env, check=False)

# =============================================================================
# 4) 結果ファイル（JSONL）を読み込む関数
# =============================================================================

def load_results_jsonl(path: Path):
    """
    JSONL（1行1JSON）を読み込んで、Pythonのリストにする。
    """
    results = []
    if not path.exists():
        return results

    with open(path, "r", encoding="utf-8") as f:
        for line in f:
            line = line.strip()
            if line:
                results.append(json.loads(line))
    return results

# =============================================================================
# 5) 結果を表形式で表示する関数
# =============================================================================

def print_results_table(results):
    """
    results（辞書のリスト）を、見やすい表で表示する。
    """
    if not results:
        return

    print("\n" + "=" * 80)
    print("物理リスト比較結果（一覧）")
    print("=" * 80)

    header = (
        f"{'PhysicsList':<18} "
        f"{'Time[s]':>9} "
        f"{'Events':>10} "
        f"{'Tracks':>10} "
        f"{'Steps':>10} "
        f"{'PPS':>12}"
    )
    print(header)
    print("-" * len(header))

    for r in results:
        print(
            f"{r['physics_list']:<18} "
            f"{r['duration_s']:>9.2f} "
            f"{r['events']:>10,} "
            f"{r['tracks']:>10,} "
            f"{r['steps']:>10,} "
            f"{r['pps']:>12.2e}"
        )

# =============================================================================
# メイン処理
# =============================================================================

if __name__ == "__main__":

    # ------------------------------------------------------------
    # 子プロセスとして呼ばれた場合：
    #   環境変数 GATE_PHYSICS_LIST がある → その物理リストで単発実行
    # ------------------------------------------------------------
    physics_list_from_env = os.environ.get("GATE_PHYSICS_LIST")
    if physics_list_from_env:
        run_simulation_and_print_report(physics_list_from_env)
        sys.exit(0)

    # ------------------------------------------------------------
    # 親プロセスとして呼ばれた場合：
    #   全物理リストを順番に「別プロセス」で実行して結果を集める
    # ------------------------------------------------------------
    results_path = Path("theme01_physics_list_results.jsonl")

    # 以前の結果があれば削除して、まっさらから始める
    if results_path.exists():
        results_path.unlink()

    # 子プロセスが結果を書き込めるように、結果ファイルの場所を環境変数で共有
    os.environ["GATE_RESULTS_FILE"] = str(results_path.resolve())

    # 物理リスト候補を順番に実行（それぞれ別プロセス）
    for physics_list in CONFIG["physics_lists"]:
        run_one_physics_list_in_subprocess(physics_list)

    # すべて終わったら結果ファイルを読み込んで一覧表示
    results = load_results_jsonl(results_path)
    print_results_table(results)

実行しますと、それぞれの物理モデルでの結果が表示されます。

また、最後にすべてのモデルの比較結果が表形式で表示されます。

物理リスト比較結果（一覧）

使用した物理モデルは（１）ハドロンモデルと（２）電磁モデルに分けて考えます。

ハドロン系の違い（QGSP_BERT、QGSP_BIC、QGSP_INCLXX、FTFP_BERT など）よりも、末尾に付いている電磁相互作用（EM）オプションの違い（EMV、EMY、EMZ）の方が、計算時間とステップ数に大きく影響します。まぁ入射エネルギーが511keVだからハドロン反応が関係しないってのがあると思います。

EMオプションを付けると、デフォルトの電磁物理が置き換わり、EMV（Option 1）は「精度は低いが速い」、EMY（Option 3）は「ガンマ線と荷電粒子輸送をより詳細に扱うため標準より遅くなる」、EMZ（Option 4）は「低エネルギーと標準パッケージから最良のモデルを選ぶ最も高精度の構成で、標準より遅い」と説明されています。

今回の結果hその通りになっているのが分かります。自分が行うシミュレーションにおいて、どの程度の細かさが必要なのかによってEMVを使うのか、EMYを使うのかを選択すると良いと思います。

上記のコードではMaterialやProduction cutについても指定していますので、あわせて説明していきます。

GATE 10は、Geant4のNISTマテリアルデータベースと、GATEカスタムマテリアルデータベースをサポートしています。自分で物質を定義することもできます。世の中に存在しないものまで作れます。

また、上記のコードではカットオフ（プロダクションカット）も指定しています。他のモンテカルロコードがエネルギーで指定するカットオフですが、Geant4ではプロダクションカットは、二次粒子の飛程を指定します。しきい値を下回った場合、二次粒子は生成されず、その場でエネルギーが付与されます。production cut は「粒子を止めるルール」ではないと説明されています。そのため、エネルギーが0になるまで追跡されます。

ちょっと理解が難しいところですね。