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

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

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

## 用意するメソッド

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

```cpp
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つです。これらをどこで呼び出せばよいかを説明していきます。

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

```cpp
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 公式サンプルどおりです。

```cpp
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);
  }
}
```

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

```cpp
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` を実装しましょう。これはイベント開始時に呼び出されます。したがって，このメソッドでは粒子を生成したときの情報が取り出せます。

```cpp
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` に直接渡しても問題はないと思います。）

```diff
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();
+}
```

```cpp
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` でしています。

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

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

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

## イベント終了時

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

```diff
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);
+ }
+}
```

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

```diff
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 以上に細かい時点での情報の取得はパフォーマンスに影響を及ぼすため推奨されていませんが，イベント開始と終了，相互作用の発生点で情報が取れればあらかたの用途には間に合うのではないでしょうか。