Geant4 からの情報の取り出し方

Geant4 は放射線シミュレーションツールキットです。シミュレーションコードの検証には，結果をその場で可視化するのがとても有用です。

Geant4 ツールキットを使ってアプリケーションを開発する際，シミュレーションの情報を取り出すには G4AnalysisManager クラスを使います。

用意するメソッド

G4AnalysisManager クラスの各メソッドを適宜呼び出すためのメソッドを用意しましょう。ここでは ResultLogger クラスとしてまとめておきます。

class ResultLogger
{
public:
  ResultLogger();
  ~ResultLogger();

  // 最初: 記録しておくデータ構造とファイルの確保
  void Book();
  // 最後: データの書き出し
  void Save();

  // イベントの開始時: 生成した光子のIDを記録する
  void LogPhotonID(int photon_id);
  // イベントの終了時: 光子が散乱した座標を記録する
  void LogScatterPosition(int photon_id, const G4ThreeVector &position);
  // イベントの終了時: 光子が検出器に入射した時のエネルギーを記録する
  void LogPhotonEnergy(int photon_id, double energy);
}

実際に取り出したい情報を記録するメソッドは最後の LogPhotonID, LogScatterePosition, LogPhotonEnergy の3つです。これらをどこで呼び出せばよいかを説明していきます。

各メソッドの実装はこんな感じです。

void ResultLogger::LogPhotonID(int photon_id) {
  fill_tuple(index_tuple_photon_id, photon_id);
}

void ResultLogger::fill_tuple(G4int tuple_id, const int photon_id) {
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->FillNtupleIColumn(tuple_id, 0, photon_id);
  analysisManager->AddNtupleRow(tuple_id);
}

Book と Save

Book と Save に関しては Geant4 公式サンプルどおりです。

void ResultLogger::Book() {
  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->SetVerboseLevel(1);

  // Create directories
  analysisManager->SetNtupleDirectoryName("ntuple");

  // Create tuple
  analysisManager->SetFirstNtupleId(0);
  index_tuple_photon_id = analysisManager->CreateNtuple("photon_id", "photon_id");
  analysisManager->CreateNtupleIColumn(index_tuple_photon_id, "photon_id");
  analysisManager->FinishNtuple(index_tuple_photon_id);

  fFactoryOn = true;

  // 他の必要なタプルなどを作成する...

  // Open an output file
  G4bool fileOpen = analysisManager->OpenFile((OutputResult::path() / "analysis").string());
  if (!fileOpen) {
    G4ExceptionDescription msg;
    msg << "Could not open file " << analysisManager->GetFileName() << ".\n";
    G4Exception("ResultLogger::Book()", "io_error", FatalException, msg);
  }
}

タプル作成は型ごとに使うメソッドが異なりますが，メソッド名に型が入っています。

void ResultLogger::Save() {
  if (!fFactoryOn) return;

  auto analysisManager = G4AnalysisManager::Instance();
  analysisManager->Write();
  analysisManager->CloseFile();

  delete G4AnalysisManager::Instance();
  fFactoryOn = false;
}

Save も公式サンプルコードどおりです。

イベント開始時

G4UserEventAction を継承した EventAction クラスを作り，仮想関数 G4UserEventAction::BeginOfEventAction をオーバーライドして EventAction::BeginOfEventAction を実装しましょう。これはイベント開始時に呼び出されます。したがって，このメソッドでは粒子を生成したときの情報が取り出せます。

void EventAction::BeginOfEventAction(const G4Event *event) {
  // initialisation per event
  auto id = event->GetEventID();
  result_logger->LogPhotonID(id);
}

僕のコードでは ResultLogger のインスタンスは EventAction のコンストラクタで渡していますが，設計次第ではもっといい方法があると思います。

イベントの途中

ここでは散乱ですが，粒子となにかの相互作用はイベントの途中で起こりますね。こういう情報は G4UserSteppingAction の子クラスとして SteppingAction クラスを用意しておき，その SteppingAction::UserSteppingAction 内で取得します。取得した情報は EventAction で一旦保持しておき，イベント終了時に ResultLogger に渡すことにします。（ResultLogger に直接渡しても問題はないと思います。）

void EventAction::BeginOfEventAction(const G4Event *event) {
  // initialisation per event
  auto id = event->GetEventID();
  result_logger->LogPhotonID(id);
-}
+ // 散乱位置記録の準備
+ scatter_position.clear();
+ scatter_position.shrink_to_fit();
+}

void EventAction::EndOfEventAction(const G4Event *event) {
  G4int eventID = event->GetEventID();
  for (const auto &position: scatter_position) {
    result_logger->LogScatterPosition(eventID, position);
  }
}

void EventAction::PushScatterPosition(const G4ThreeVector &position) {
  scatter_position.push_back(position);
}

散乱位置を記録するための std::vector<G4ThreeVector> scatter_position の初期化を BeginOfEventAction でしています。

void SteppingAction::UserSteppingAction(const G4Step *step) {

  // 物理過程の name が "your_process" であるときのみ，相互作用した直後の座標を記録する
  if (step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "your_process") {
    fEventAction->PushScatterPosition(step->GetPostStepPoint()->GetPosition());
  };
}

コード規約が全然守られていないのは今は無視してください。

イベント終了時

ここでは検出器に入ったら終わりにします。正確には EventAction::EndOfEventAction が呼び出される時点，すなわち粒子の飛跡が 0 になるタイミングが終了です。

void SteppingAction::UserSteppingAction(const G4Step *step) {

  // 物理過程の name が "your_process" であるときのみ，相互作用した直後の座標を記録する
  if (step->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "your_process") {
    fEventAction->PushScatterPosition(step->GetPostStepPoint()->GetPosition());
  };
-}
+ // 粒子が検出器ボリュームに入ったら光子のエネルギーを記録する
+ G4VPhysicalVolume *volume = step->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
+ if (volume == fDetConstruction->GetDetectorPV()) {
+   const auto track = step->GetTrack();
+   DoDetectorTasks(*track);
+ }
+}

void SteppingAction::DoDetectorTasks(const G4Track &track) {
  fEventAction->RegisterPhotonEnergy(track.GetKineticEnergy() / keV);
}

void EventAction::EndOfEventAction(const G4Event *event) {
  G4int eventID = event->GetEventID();

  for (const auto &position: scatter_position) {
    result_logger->LogScatterPosition(eventID, position);
  }
-}
+ result_logger->LogPhotonEnergy(eventID, photon_energy);
+}

このようにして，節目節目で様々な情報を取り出すことができます。SteppingAction 以上に細かい時点での情報の取得はパフォーマンスに影響を及ぼすため推奨されていませんが，イベント開始と終了，相互作用の発生点で情報が取れればあらかたの用途には間に合うのではないでしょうか。