Geant4でピクセル検出器を作って画像認識したい（Sensitive Volume編）

はじめに

途中で切れている記事の下書きがあったので、公開だけしておきます。
このプロジェクトから離れて久しいため、補足は多分ないです…

物理的背景

γ線はコンプトン散乱と光電吸収を起こす。一方でβ線は制動放射をしてγ線を放出するか、そのまま吸収される。
なので、とりあえずはγ線以外の粒子を選んで取れば、現実に即した検出器になりそう。

情報の格納

粒子の発射一発毎にピクセル上のエネルギー分布の画像が欲しい。
SensitiveVolume自体は粒子のステップ毎に呼ばれることから、イベント毎に考えることには不適切。ラン（指定した回数粒子を発射する）の間保有されるHitCollectionに情報を格納する。

どの物理量が欲しいか

必須のデータ

• TrackID(生成された粒子の識別番号)
• ParticleID（粒子の種類のID）
• Row & Col（どのピクセルにいるか）
• Energy Deposit(落としたエネルギーの量)

あると便利かも or 後に使うかものデータ

• copyNumber（複数の検出器を使う時のための識別番号）
• nStep（そのトラックでのステップ番号）
• nhit（ヒットの全体での通し番号）

SensitiveVolumeの組み立て

情報の初期化

void SensitiveVolume::Initialize(G4HCofThisEvent* HCTE)
{
  no_Step = 0;
  // create hit collection(s)
  hitsCollection= new PixelHitsCollection(SensitiveDetectorName,
					  collectionName[0]); 
  
  // push H.C. to "Hit Collection of This Event"
  G4int hcid= GetCollectionID(0);
  HCTE-> AddHitsCollection(hcid, hitsCollection);

}

粒子の情報の格納先になるHitCollectionを作っておく。

粒子の情報を得る

G4bool SensitiveVolume::ProcessHits(G4Step* aStep, G4TouchableHistory*)
{
  // energy deposit in this step and its accumulation over steps
  G4double edep = aStep->GetTotalEnergyDeposit();
  
  // Retrieve information from the track object
  G4int nTrack = aStep->GetTrack()->GetTrackID();
  G4int nStep  = aStep->GetTrack()->GetCurrentStepNumber();
  G4String particleName = aStep->GetTrack()->GetDefinition()->GetParticleName();
  G4int pid  = aStep->GetTrack()->GetDefinition()->GetPDGEncoding();

現在のステップから落としたエネルギー、トラックID、ステップ番号、粒子の種類とその名前が得られる。
前のステップとの差からエネルギーが計算されている。

  // Get PreStepPoint and TouchableHandle objects
  G4StepPoint*      preStepPoint = aStep->GetPreStepPoint();
  G4TouchableHandle theTouchable = preStepPoint->GetTouchableHandle();
  G4int copyNum = theTouchable->GetCopyNumber();
 
  G4int mid = 0; //mother ID
  G4int nhit = hitsCollection->entries();

  // Touchable information: Volume name and copy# at the current step
  G4String volumeName = theTouchable->GetVolume()->GetName();  

論理物体の境界では必ずステップが起こるので、その時は直前のステップの座標を元にするようにする。
theTouchableは自分で配置した物体の総称で、そこから今いる論理物体のIDや名前が得られる。
midは今は使わないが、複数の検出器を使う時は、
theTouchable->GetcopyNumber(3)
で、その検出器のIDを持ってくる。

//LogVol_PixElmt
  G4int row = theTouchable->GetCopyNumber();  // was copyNo

  G4String motherVolumeName;
  G4int col; // was motherCopyNo;
  if (volumeName != "LogVol_PixElmt" ) {
  motherVolumeName = theTouchable->GetVolume(1)->GetName();
  col = theTouchable->GetCopyNumber(1);
  } 
  else {
  motherVolumeName = "None";
  col = 0;
  }

粒子がどのピクセルにいるのかを
今のlogVol_ElmtのlogVol_EnvLの中でのコピーナンバー(row)
そのlogVol_EnvLのlogVol_EnvGの中でのコピーナンバー(col)
を得ることで特定する。

  // a primary track with non-zero energy deposition is registerd
  if ( edep==0. ) return false;
  if (particleName=="gamma") return false;

  PixelHit* hit = new PixelHit( nhit, edep, no_Step, pid, nTrack, mid, copyNum, row, col);
  hitsCollection->insert(hit);
  return true;
}

エネルギーを落とさない過程、γ線の過程はHitCollectionに格納しないようにする。
PixelHitは後述するコードの中で定義される、HitCollectionの中に格納できる形にする関数。

PixelHit

ステップ一つ一つをHitCollectionに格納できるようにする。

PixelHit.hh

#ifndef CAL_HIT_H
#define CAL_HIT_H
 
#include "G4VHit.hh"
#include "G4THitsCollection.hh"
#include "G4Allocator.hh"

class PixelHit : public G4VHit {
private:
  G4int id;
  G4double edep;
  G4int nstep;
  G4int pid;
  G4int tid;
  G4int mid;
  G4int copy;
  G4int row;
  G4int col;
  //

  //
public:
  PixelHit();
  PixelHit(G4int aid, G4double aedep, G4int anstep, 
	   G4int apid, G4int atid, G4int amid, G4int acopy,
	   G4int arow, G4int acol);
  virtual ~PixelHit();

  // copy constructor & assignment operator
  PixelHit(const PixelHit& right);
  const PixelHit& operator=(const PixelHit& right);
  G4int operator==(const PixelHit& right) const;
  
  // new/delete operators
  void* operator new(size_t);
  void operator delete(void* aHit);
  
  // set/get functions
  void SetID(G4int aid) { id = aid; }
  void SetEdep(G4double aedep){ edep = aedep; }
  void SetNstep(G4int anstep ){ nstep = anstep; }
  void SetPID(G4int apid) { pid = apid; }
  void SetTID(G4int atid) { tid = atid; }
  void SetMID(G4int amid) { mid = amid; }
  void SetNcopy(G4int acopy) { copy = acopy; }
  void SetRow(G4int arow) {row = arow;}
  void SetCol(G4int acol) {col = acol;}
  
  G4int    GetID() const { return id; }
  G4double GetEdep() const { return edep; }
  G4int    GetNstep() const { return nstep; }
  G4int    GetPID() const { return pid; }
  G4int    GetTID() const { return tid; }
  G4int    GetMID() const { return mid; }
  G4int    GetNcopy() const { return copy; }
  G4int    GetRow() const { return row; }
  G4int    GetCol() const { return col; }

  // methods
  virtual void Draw();
  virtual void Print();
};

// ====================================================================
// inline functions
// ====================================================================
inline PixelHit::PixelHit(const PixelHit& right)
  : G4VHit()
{
  id= right.id;
  edep= right.edep;
  nstep= right.nstep;
  pid= right.pid;
  tid= right.tid;
  mid= right.mid;
  copy= right.copy;
  row = right.row;
  col = right.col;
}

inline const PixelHit& PixelHit::operator=(const PixelHit& right)
{
  id= right.id;
  edep= right.edep;
  nstep= right.nstep;
  pid= right.pid;
  tid= right.tid;
  mid= right.mid;
  copy= right.copy;
  row = right.row;
  col = right.col;

  return *this;
}

inline G4int PixelHit::operator==(const PixelHit& right) const 
{
   return (this==&right) ? 1 : 0; 
}

typedef G4THitsCollection<PixelHit> PixelHitsCollection;
extern G4Allocator<PixelHit> PixelHitAllocator; 

inline void* PixelHit::operator new(size_t)
{
  void* aHit= (void*)PixelHitAllocator.MallocSingle();
  return aHit;
}

inline void PixelHit::operator delete(void* aHit)
{
  PixelHitAllocator.FreeSingle((PixelHit*) aHit);
}

#endif