Geant4(その2) 世界一わかりやすい?コード例

ソフトウェア

今回の話題

前回Geant4をインストールしました.今回はユーザコード例(C++)を紹介します.Geant4のexamplesのコード結構難しい…という方向けに,最小限の骨格を明確にしたいと思います.

下図が出来上がりの画面イメージです.リニアックの電子線ターゲット風に,真空中に銅とタングステンの板を重ねて置き,上から電子線を当てます.電子(赤線)が散乱したり制動X線(緑線)が出たりします.

Geant4サンプルアプリ 実行画面イメージ
実行画面イメージ

ファイル構成

今回のコード例のファイル構成とその役割は以下のとおりです.

  • myapp.cc … メインプログラム
  • include/ … ユーザ定義クラスのI/F
    • MyDetectorConstruction.hh … 物体モデル
    • MyPrimaryGeneratorAction.hh … 線源モデル
    • MySteppingAction.hh … イベント処理
  • src/ … ユーザ定義クラスの実装
    • MyDetectorConstruction.cc … 物体モデル(形状,材質,配置)
    • MyPrimaryGeneratorAction.cc … 線源モデル(粒子のサンプリング)
    • MySteppingAction.cc … イベント処理(データの集計・出力処理)

ユーザ定義の各クラス(My○○)は,いわゆるGeant4カーネルとうまく連携できるように,それぞれ役割に応じたインタフェースクラス(G4V○○という名前,Vはvirtual classの意)を基底とし必要なI/Fを実装します.以下,順に見ていきましょう.

メインプログラム

ポイントは,Run managerで計算体系を構築する(36-39行目),UI managerでユーザレベルのコマンドを実行する(50行目),の2点です.

#include "MyDetectorConstruction.hh"
#include "MyPrimaryGeneratorAction.hh"
#include "MySteppingAction.hh"
#include "G4PhysListFactory.hh"
#include "G4RunManager.hh"
#include "G4UImanager.hh"
#include "G4VisExecutive.hh"
#include "G4UIExecutive.hh"
#include "G4SystemOfUnits.hh"  // 単位換算係数(mm,MeVなど)
#include "CLHEP/Random/Randomize.h"
#include <map>
#include <string>
#include <cstdlib>
#include <iostream>

int main(int argc,char* argv[]){
  // デフォルト設定
  std::map<std::string,double> cfg = { {"E",6*MeV}, {"R",0.1*mm}, {"D0",20*mm}, {"D1",2*mm}, {"D2",1*mm} };  // 設定項目=>値
  std::string macro="/dev/stdin", phys="FTFP_BERT";  // マクロ, 物理リスト
  bool visual=false;  // 可視化ON/OFF

  // オプション設定                                                                                                                                                                    
  for(int i=1; i<argc; i++){
    std::string arg = std::string(argv[i]);
    if (arg[0]!='-') macro=arg;
    else if (arg=="-d") cfg["D2"]= atof(argv[++i])*mm;   // W層厚[mm]
    else if (arg=="-e") cfg["E"] = atof(argv[++i])*MeV;  // ビームエネルギー[MeV]
    else if (arg=="-r") cfg["R"] = atof(argv[++i])*mm;   // ビーム径[mm]
    else if (arg=="-s") CLHEP::HepRandom::setTheSeed(atoi(argv[++i]));  // 乱数初期化
    else if (arg=="-v") visual = true;  // 可視化ON
    else { std::cerr<<"Usage: "<<argv[0]<<" [-d MM] [-e MEV] [-r MM] [-s SEED] [-v] [MACRO]"<<std::endl; return 0; }
  }

  // 計算体系構築 - 物理リスト, 物体モデル, 線源モデル, イベント処理                                                                             
  G4RunManager runManager;
  runManager.SetUserInitialization(G4PhysListFactory().GetReferencePhysList(phys));
  runManager.SetUserInitialization(new MyDetectorConstruction(cfg));
  runManager.SetUserAction(new MyPrimaryGeneratorAction(cfg));
  runManager.SetUserAction(new MySteppingAction(cfg));
  runManager.Initialize();

  // シミュレーション実行 - 外部マクロ実行など                                                                      
  if ( visual==true ){  // visual & interactive
    G4VisExecutive visManager;
    visManager.Initialize();  // 表示系準備
    G4UImanager::GetUIpointer()->ApplyCommand("/control/execute "+macro);
    G4UIExecutive(argc,argv).SessionStart();  // 対話セッション
  }
  else {  // batch
    G4UImanager::GetUIpointer()->ApplyCommand("/control/execute "+macro);
  }
}

計算体系構築は,Geant4カーネルに下記4点を渡すことでおこないます.I/Fとしては,カッコ内のクラスを基底にもつオブジェクトとしてRun managerに渡します.

  • 線源モデル (G4VUserPrimaryGeneratorAction)
  • 物体モデル (G4VUserDetectorConstruction)
  • 各種イベント処理 (G4○○Action)
  • 物理リスト (G4VUserPhysicsList)

「これこれこういう線源粒子を,こういう物体に照射します.こういうイベント(散乱とか物体境界超えとか)が生じたらこう処理(集計や出力など)します.各種物理過程の計算には指定した物理リストのアルゴリズムを使います」というわけです.これで計算の準備が整いました.

計算実行するには,”/run/beamOn 100″ のような「コマンド」をUI managerに直接渡すか,50行目のように「マクロ」から読込みます(コマンド,マクロは後日紹介).するとあとはGeant4カーネルが線源粒子を生成し,物体内で輸送計算し,指定通り集計・出力処理してくれます.

なお物理リスト(physics list)について詳しくはこちら.Geant4は各種物理過程について適用条件が異なる複数の計算モデルを持っているので,どれを使うかユーザが指定します.今回は定番セットのFTFP_BERTを選択.

線源モデル

線源粒子を確率的に生成します.Geant4カーネル内で線源粒子生成時に呼ばれるイベントハンドラです(ヘッダ11行目).

#pragma once
#include "G4VUserPrimaryGeneratorAction.hh"  // 線源モデル基底クラス
#include "G4ParticleGun.hh"
#include <map>
#include <string>

class MyPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction {
  public:
    MyPrimaryGeneratorAction(const std::map<std::string,double>& _cfg);
    virtual ~MyPrimaryGeneratorAction();
    virtual void GeneratePrimaries(G4Event*);  // 粒子生成イベントハンドラ

  private:
    G4ParticleGun gun;
    std::map<std::string,double> cfg;
};
#include "MyPrimaryGeneratorAction.hh"
#include "G4Electron.hh"
#include "CLHEP/Random/RandGauss.h"

inline double nrand(){ return CLHEP::RandGauss::shoot(); }  // 標準正規分布

MyPrimaryGeneratorAction::MyPrimaryGeneratorAction(const std::map<std::string,double>& _cfg)
  : G4VUserPrimaryGeneratorAction(), gun(), cfg(_cfg) {}

MyPrimaryGeneratorAction::~MyPrimaryGeneratorAction() {}

void MyPrimaryGeneratorAction::GeneratePrimaries(G4Event* ev) {
  // 粒子属性を所望の確率分布からサンプリング
  double d0=cfg["D0"], r=cfg["R"], e=cfg["E"];
  gun.SetParticleDefinition(G4Electron::Definition());
  gun.SetParticleEnergy(e);
  gun.SetParticlePosition({r*nrand(), r*nrand(), d0/2});
  gun.SetParticleMomentumDirection({0,0,-1});
  gun.GeneratePrimaryVertex(ev);
}

処理概要:

  • 粒子属性を所望の確率分布からサンプリングします(15-18行目).今回は位置以外固定.
  • 粒子を生成し,Geant4カーネルに渡します(19行目).

Geant4カーネルとの連携の仕組:

  • Geant4カーネルが線源モデルとして扱えるよう,G4VUserPrimaryGeneratorAction を基底クラスとします(ヘッダ7行目).
  • Geant4カーネルに線源モデルを登録します(myapp.cc 38行目).
  • Geant4カーネル内の計算で,線源粒子生成時に GeneratePrimaries() が呼ばれます.

線源粒子は,粒子の種類,運動エネルギー(=質量含まず),位置,進行方向で規定します.固定値の属性はコンストラクタで1度設定するだけでもOKです.

物体モデル

ビーム照射を受ける各物体の形状,材質,配置を規定します.これもGeant4カーネル内で物体生成/更新時に呼ばれるイベントハンドラと言って良いでしょう(ヘッダ10行目).

#pragma once
#include "G4VUserDetectorConstruction.hh"  // 物体モデル基底クラス
#include <map>
#include <string>

class MyDetectorConstruction : public G4VUserDetectorConstruction {
  public:
    MyDetectorConstruction(const std::map<std::string,double>& _cfg);
    virtual ~MyDetectorConstruction();
    virtual G4VPhysicalVolume* Construct();  // 物体モデル生成

  private:
    std::map<std::string,double> cfg;
};
#include "MyDetectorConstruction.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4Exception.hh"

// 物質取得ユーティリティ
G4Material* get_material(const std::string& name) {
  // NIST物質DBから検索
  auto* nist = G4NistManager::Instance();
  auto* mat = nist->FindOrBuildMaterial(name);
  if (mat!=0) return mat;  // DB内蔵orユーザ定義済み物質

  // ユーザ定義物質を追加
  std::vector<int>    natom;  // 組成(原子数)
  std::vector<double> wfrac;  // 組成(重量比)
  if (name=="Vacuum") mat = nist->ConstructNewMaterial(name, {"N","O"}, wfrac={0.7,0.3}, 1e-8*gram/cm3);
  if (name=="Water" ) mat = nist->ConstructNewMaterial(name, {"H","O"}, natom={2,1}, 1*gram/cm3);
  if (mat==0) G4Exception("get_material()", ("unknown material="+name).c_str(), FatalException, "");
  return mat;
}

MyDetectorConstruction::MyDetectorConstruction(const std::map<std::string,double>& _cfg)
  : G4VUserDetectorConstruction(), cfg(_cfg) {}

MyDetectorConstruction::~MyDetectorConstruction() {}

G4VPhysicalVolume* MyDetectorConstruction::Construct() {
  double d0=cfg["D0"], d1=cfg["D1"], d2=cfg["D2"];  // d0:領域0の1辺,d1,d2:領域1,2の厚さ

  // 全体領域0(真空) - 形状, 材質, 配置を指定
  auto* bx0 = new G4Box("box0", d0/2, d0/2, d0/2);
  auto* lv0 = new G4LogicalVolume(bx0, get_material("Vacuum"), "World");
  auto* pv0 = new G4PVPlacement(0, {0,0,0}, lv0, "pvol0", 0, false, 0, true);
  // (回転*,移動,論理領域*,物理領域名,親領域*,コピー使用,コピー番号,重複確認)

  // 子領域1(Cu層)
  auto* bx1 = new G4Box("box1", d0/2, d0/2, d1/2);
  auto* lv1 = new G4LogicalVolume(bx1, get_material("G4_Cu"), "Region1");
  /*pv1 =*/ new G4PVPlacement(0, {0,0,0}, lv1, "pvol1", lv0, false, 0, true);

  // 孫領域2(W層)
  auto* bx2 = new G4Box("box2", d0/2, d0/2, d2/2);
  auto* lv2 = new G4LogicalVolume(bx2, get_material("G4_W"), "Region2");
  /*pv2 =*/ new G4PVPlacement(0, {0,0,(d1-d2)/2}, lv2, "pvol2", lv1, false, 0, true);

  return pv0;  // 全体を返す
}

処理概要:

  • 各材質を生成します(12-13, 19-20行目).NIST物質DBから選ぶか,組成と密度で規定.
  • 各領域を生成します(34-36,40-42,45-47行目).形状,材質,配置で規定.
  • 全体領域を返します(49行目).Geant4カーネルとの連携上の決まり.

Geant4カーネルとの連携の仕組:

  • Geant4カーネルが物体モデルとして扱えるよう,G4VUserDetectorConstruction を基底クラスとします(ヘッダ6行目).
  • Geant4カーネルに物体モデルを登録します(myapp.cc 37行目).
  • Geant4カーネル内の計算で,物体モデル生成時に Construct() が呼ばれます.

物体モデルのデータ構造は次のようになっていて,コード例ではこれを設定しています.大雑把に言って,物理領域=論理領域+配置,論理領域=形状+材質,材質=組成+密度です.各物理領域は親子関係によりツリー構造をなしており,全体領域(親がnullのルートノード)からすべての子孫領域の情報を取得できます.

  • 物理領域 (G4VPhysicalVolume*)
    • 論理領域 (G4LogicalVolume*)
      • 形状 (G4VSolid*)
      • 材質 (G4Material*)
        • 組成
        • 密度
      • 子の物理領域 (G4VPhysicalVolume*[])
    • 親に対する回転量 (G4RotationMatrix*)
    • 親に対する移動量 (G4ThreeVector)
    • 親の論理領域 (G4VLogicalVolume*)

形状を表す具象クラス(全てG4VSolidから派生)はいろいろあります.たとえば,

  • 直方体(G4Box),球面(G4Sphere),管(G4Tubs),円錐台(G4Cons)などの単純形状.
  • 既存形状の論理演算形状(G4IntersectionSolid, G4SubtractionSolid, G4UnionSolid).

材質(G4Material)は主に組成と密度で規定しますが,例のように,

  • NIST物質DBには,単一元素物質,水,空気,人体組成ほかがあります(12-13行目).
  • Geant4では完全な真空を定義できません.低密度の気体などとします(19行目).

子領域が親領域からハミ出ていたり,親子関係にない領域が重なっていたりすると,計算がおかしくなる(落ちはしない)のでご注意を.

イベント処理

Geant4カーネル内で呼ばれるイベントハンドラです.今回はstep時(輸送計算で粒子が移動したとき)に呼ばれるタイプ(ヘッダ10行目).

#pragma once
#include "G4UserSteppingAction.hh"  // step時イベントハンドラ基底クラス
#include <map>
#include <string>

class MySteppingAction : public G4UserSteppingAction {
  public:
    MySteppingAction(const std::map<std::string,double>& _cfg);
    virtual ~MySteppingAction();
    virtual void UserSteppingAction(const G4Step*);  // step時イベントハンドラ

  private:
    std::map<std::string,double> cfg;
};
#include "MySteppingAction.hh"
#include "G4Step.hh"
#include "G4SystemOfUnits.hh"
#include <iostream>

MySteppingAction::MySteppingAction(const std::map<std::string,double>& _cfg)
  : G4UserSteppingAction(), cfg(_cfg) {}

MySteppingAction::~MySteppingAction() {}

// step時イベントハンドラ - 領域境界で状態出力する例
void MySteppingAction::UserSteppingAction(const G4Step* step) {
  const G4StepPoint* p1 = step->GetPostStepPoint();
  if (p1->GetStepStatus() != G4StepStatus::fGeomBoundary) return;

  const G4ThreeVector& pos = p1->GetPosition();
  const G4ThreeVector& dir = p1->GetMomentumDirection();
  std::cout << p1->GetPhysicalVolume()->GetName() << "\t"
            << p1->GetCharge() << "\t"
            << p1->GetKineticEnergy()/MeV << "\t"
            << pos.x()/mm << "\t" << pos.y()/mm << "\t" << pos.z()/mm << "\t"
            << dir.x() << "\t" << dir.y() << "\t" << dir.z()
            << std::endl;
}

処理概要:

  • 処理対象を絞ります(14行目).ここでは次の領域の境界に到達したケース.
  • 粒子情報などを出力してみます(18行目).実用性は考慮せず…

Geant4カーネルとの連携の仕組:

  • Geant4カーネルがstep時のイベントハンドラとして扱えるよう,G4UserSteppingAction を基底クラスとします(ヘッダ6行目).
  • Geant4カーネルにイベントハンドラを登録します(myapp.cc 39行目).
  • Geant4カーネル内の計算で,step発生時に UserSteppingAction() が呼ばれます.

イベントハンドラとして使えるクラスはほかにも下記などがあります.それぞれ呼ばれるタイミングが違うので(下記は怪しいかも),うまく使い分け・併用してください.5番目のSensitive Detectorについては,登録手続きなどやや複雑ですが,これを駆使してスコアリングなどするとプロっぽい感じに?

  • G4UserRunAction … Run(ビーム照射)の開始時と終了時
    • G4UserEventAction … Event(線源粒子・二次粒子の一部始終)の開始時と終了時
      • G4UserTrackingAction … Tracking(個々の粒子軌跡)の始点と終点
        • G4UserSteppingAction … Stepping(粒子軌跡のある点から次まで)ごと
        • G4VSensitiveDetector …. 指定した論理領域中のStepごと

まとめ

以上で簡単なGeant4アプリが完成しましたが,如何だったでしょうか.

割愛した内容は下記を含めたくさんあります.まだ先は長い…

  • スコアリング関連(Sensitive Detector,Hitsなど)
  • 物体モデル関連(回転,繰返し,CADデータ読込みなど)
  • 並列計算関連(マルチスレッド,MPI)
  • コマンドのカスタマイズ(Messenger)
  • 描画,解析関連

最後に,冒頭述べた「最小限の骨格を」の観点で無理やりまとめておきます.

「ユーザ定義クラスは,みんなGeant4カーネルのイベントハンドラだった!」

次回予告

「Geant4コードのビルドと実行」ということで,今回のコードを動かしてみます.


Accutheraでは加速器・医療システム・機械学習・放射線シミュレーションなどの専門技術を軸に開発やコンサルティングを承っております。お困りの案件がございましたら、こちら からお気軽にお問合せください。

タイトルとURLをコピーしました