geant4 灵敏探测器

geant4 灵敏探测器

---

 如何模拟探测器的响应？

1. SteppingAction 统计每一个 step 在探测器中沉积的能量。
2. 指定探测器为灵敏探测器，统计探测器内每一个 step 的能量沉积。

如何设置一个灵敏探测器？

1. 指定某逻辑体为灵敏探测器
2. 将灵敏探测器纳入 SD 管理器

　　代码举例，取自 example B2a ---- 设置灵敏探测器，将灵敏探测器纳入 SD 管理器，指定某逻辑体为灵敏探测器。

void B2aDetectorConstruction::ConstructSDandField()
{
  // Sensitive detectors

  G4String trackerChamberSDname = "B2/TrackerChamberSD";
  B2TrackerSD* aTrackerSD = new B2TrackerSD(trackerChamberSDname,
                                            "TrackerHitsCollection");
  G4SDManager::GetSDMpointer()->AddNewDetector(aTrackerSD);
  // Setting aTrackerSD to all logical volumes with the same name 
  // of "Chamber_LV".
  SetSensitiveDetector("Chamber_LV", aTrackerSD, true);

  // Create global magnetic field messenger.
  // Uniform magnetic field is then created automatically if
  // the field value is not zero.
  G4ThreeVector fieldValue = G4ThreeVector();
  fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
  fMagFieldMessenger->SetVerboseLevel(1);
  
  // Register the field messenger for deleting
  G4AutoDelete::Register(fMagFieldMessenger);
}

• 第 6 - 9 行，设置灵敏探测器， 并将其纳入灵敏探测器管理器。其中 “TrackerHitsCollection” 是设置的 hitcollectionName，便于在 G4HCofThisEvent 容器中索引。
• 第 11 行，使用 方法 SetSensitiveDetector()设置逻辑体为灵敏探测器，该方法的构造函数如下，这里是将所有名字为 Chamber_LV 的逻辑体统统设置为灵敏探测器。

  void SetSensitiveDetector(const G4String& logVolName,
                  G4VSensitiveDetector* aSD,G4bool multi=false);
  void SetSensitiveDetector(G4LogicalVolume* logVol,
                  G4VSensitiveDetector* aSD);

如何从灵敏探测器中提取信息？

Hit、HitCollection、Event

• Hit：记录 SD 中每个 G4Step 中的信息(能量，时间等)，因此一个 SD 中将产生大量Hit对象。
• HitCollection：一个 SD 中产生大量的 Hit，可以被记录在一个或者多个容器中， 该容器即为 HitCollection 对象。
• G4HCofThisEvent：G4Event 的一个成员对象，为一容器；多个 SD 的多个 HitCollection 对象储存在该容器中。

一个 Hit 类的举例，代码取自 example B2a

/// Tracker hit class
///
/// It defines data members to store the trackID, chamberNb, energy deposit,
/// and position of charged particles in a selected volume:
/// - fTrackID, fChamberNB, fEdep, fPos

class B2TrackerHit : public G4VHit
{
  public:
    B2TrackerHit();
    B2TrackerHit(const B2TrackerHit&);
    virtual ~B2TrackerHit();

    // operators
    const B2TrackerHit& operator=(const B2TrackerHit&);
    G4int operator==(const B2TrackerHit&) const;

    inline void* operator new(size_t);
    inline void  operator delete(void*);

    // methods from base class
    virtual void Draw();
    virtual void Print();

    // Set methods
    void SetTrackID  (G4int track)      { fTrackID = track; };
    void SetChamberNb(G4int chamb)      { fChamberNb = chamb; };
    void SetEdep     (G4double de)      { fEdep = de; };
    void SetPos      (G4ThreeVector xyz){ fPos = xyz; };

    // Get methods
    G4int GetTrackID() const     { return fTrackID; };
    G4int GetChamberNb() const   { return fChamberNb; };
    G4double GetEdep() const     { return fEdep; };
    G4ThreeVector GetPos() const { return fPos; };

  private:

      G4int         fTrackID;
      G4int         fChamberNb;
      G4double      fEdep;
      G4ThreeVector fPos;
};

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

typedef G4THitsCollection<B2TrackerHit> B2TrackerHitsCollection;

extern G4ThreadLocal G4Allocator<B2TrackerHit>* B2TrackerHitAllocator;

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

inline void* B2TrackerHit::operator new(size_t)
{
  if(!B2TrackerHitAllocator)
      B2TrackerHitAllocator = new G4Allocator<B2TrackerHit>;
  return (void *) B2TrackerHitAllocator->MallocSingle();
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

inline void B2TrackerHit::operator delete(void *hit)
{
  B2TrackerHitAllocator->FreeSingle((B2TrackerHit*) hit);
}

• 在上述例子中，提取了 tarckID，Chamber number，沉积能量，位置信息。
• 选择在 hit 中记录哪些量，是在 hit 类中明确的。
• 第 47 行，将存放 B2TrackerHit 的容器类型重命名为 B2TrackerHitsCollection，方便接下来定义使用。

一个 SD 类的举例，代码取自 example B2a

#ifndef B2TrackerSD_h
#define B2TrackerSD_h 1

#include "G4VSensitiveDetector.hh"

#include "B2TrackerHit.hh"

#include <vector>

class G4Step;
class G4HCofThisEvent;

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

/// B2Tracker sensitive detector class
///
/// The hits are accounted in hits in ProcessHits() function which is called
/// by Geant4 kernel at each step. A hit is created with each step with non zero 
/// energy deposit.

class B2TrackerSD : public G4VSensitiveDetector
{
  public:
    B2TrackerSD(const G4String& name, 
                const G4String& hitsCollectionName);
    virtual ~B2TrackerSD();
  
    // methods from base class
    virtual void   Initialize(G4HCofThisEvent* hitCollection);
    virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
    virtual void   EndOfEvent(G4HCofThisEvent* hitCollection);

  private:
    B2TrackerHitsCollection* fHitsCollection;
};

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#endif

#include "B2TrackerSD.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4SDManager.hh"
#include "G4ios.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B2TrackerSD::B2TrackerSD(const G4String& name,
                         const G4String& hitsCollectionName) 
 : G4VSensitiveDetector(name),
   fHitsCollection(NULL)
{
  collectionName.insert(hitsCollectionName);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B2TrackerSD::~B2TrackerSD() 
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B2TrackerSD::Initialize(G4HCofThisEvent* hce)
{
  // Create hits collection

  fHitsCollection 
    = new B2TrackerHitsCollection(SensitiveDetectorName, collectionName[0]); 

  // Add this collection in hce

  G4int hcID 
    = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
  hce->AddHitsCollection( hcID, fHitsCollection ); 
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

G4bool B2TrackerSD::ProcessHits(G4Step* aStep, 
                                     G4TouchableHistory*)
{  
  // energy deposit
  G4double edep = aStep->GetTotalEnergyDeposit();

  if (edep==0.) return false;

  B2TrackerHit* newHit = new B2TrackerHit();

  newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
  newHit->SetChamberNb(aStep->GetPreStepPoint()->GetTouchableHandle()
                                               ->GetCopyNumber());
  newHit->SetEdep(edep);
  newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());

  fHitsCollection->insert( newHit );

  //newHit->Print();

  return true;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B2TrackerSD::EndOfEvent(G4HCofThisEvent*)
{
  if ( verboseLevel>1 ) { 
     G4int nofHits = fHitsCollection->entries();
     G4cout << G4endl
            << "-------->Hits Collection: in this event they are " << nofHits 
            << " hits in the tracker chambers: " << G4endl;
     for ( G4int i=0; i<nofHits; i++ ) (*fHitsCollection)[i]->Print();
  }
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

• 当一个 step 的能量沉积不是 0 的时候，就会创建一个 hit，言外之意是，不是所有的 step 都有能量沉积。
• 灵敏探测器类继承自 G4VSensitiveDetector 类，其中必须实现的方法是 ProcessHits(G4Step*aStep,G4TouchableHistory*ROhist) ，该方法是从 step 中提取信息，写入 hit 中，并且将 hit 放入到 HitCollection 的容器中。
• 第 15 行，必须要设置 "collectionName" string vector，原文：In the derived class constructor, name(s) of hits collection(s) which are made by the sensitive detector must be set to "collectionName" string vector. 
• EndOfEvent(G4HCofThisEvent*) 不是一定要实现的方法，在该例中，在 EndOfEvent() 方法中提取了每一个 Event 共产生了多少个hit.

一个 EventAction 类的举例，用于提取灵敏探测器的信息，代码取自 example B5（第一段代码取自 DetectorConstruction.cc）

void B5DetectorConstruction::ConstructSDandField()
{
  // sensitive detectors -----------------------------------------------------
  auto sdManager = G4SDManager::GetSDMpointer();
  G4String SDname;
  
  auto hodoscope1 = new B5HodoscopeSD(SDname="/hodoscope1");
  sdManager->AddNewDetector(hodoscope1);
  fHodoscope1Logical->SetSensitiveDetector(hodoscope1);

  auto hodoscope2 = new B5HodoscopeSD(SDname="/hodoscope2");
  sdManager->AddNewDetector(hodoscope2);
  fHodoscope2Logical->SetSensitiveDetector(hodoscope2);
  
  auto chamber1 = new B5DriftChamberSD(SDname="/chamber1");
  sdManager->AddNewDetector(chamber1);
  fWirePlane1Logical->SetSensitiveDetector(chamber1);

  auto chamber2 = new B5DriftChamberSD(SDname="/chamber2");
  sdManager->AddNewDetector(chamber2);
  fWirePlane2Logical->SetSensitiveDetector(chamber2);
  
  auto emCalorimeter = new B5EmCalorimeterSD(SDname="/EMcalorimeter");
  sdManager->AddNewDetector(emCalorimeter);
  fCellLogical->SetSensitiveDetector(emCalorimeter);
  
  auto hadCalorimeter = new B5HadCalorimeterSD(SDname="/HadCalorimeter");
  sdManager->AddNewDetector(hadCalorimeter);
  fHadCalScintiLogical->SetSensitiveDetector(hadCalorimeter);

  // magnetic field ----------------------------------------------------------
  fMagneticField = new B5MagneticField();
  fFieldMgr = new G4FieldManager();
  fFieldMgr->SetDetectorField(fMagneticField);
  fFieldMgr->CreateChordFinder(fMagneticField);
  G4bool forceToAllDaughters = true;
  fMagneticLogical->SetFieldManager(fFieldMgr, forceToAllDaughters);

  // Register the field and its manager for deleting
  G4AutoDelete::Register(fMagneticField);
  G4AutoDelete::Register(fFieldMgr);
}    

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

• 注意到，在这里初始化各个灵敏探测器时，没有设置其对应的 HitCollectionName，这是因为在对应的灵敏探测器类中，已经设置了该名称。

  B5EmCalorimeterSD::B5EmCalorimeterSD(G4String name)
  : G4VSensitiveDetector(name), 
    fHitsCollection(nullptr), fHCID(-1)
  {
    collectionName.insert("EMcalorimeterColl");
  }

#ifndef B5EventAction_h
#define B5EventAction_h 1


#include "G4UserEventAction.hh"
#include "globals.hh"

#include <vector>

/// Event action

class B5EventAction : public G4UserEventAction
{
public:
    B5EventAction();
    virtual ~B5EventAction();
    
    virtual void BeginOfEventAction(const G4Event*);
    virtual void EndOfEventAction(const G4Event*);

    std::vector<G4double>& GetEmCalEdep() { return fEmCalEdep; }
    std::vector<G4double>& GetHadCalEdep() { return fEmCalEdep; }
    
private:
    G4int fHodHC1ID;
    G4int fHodHC2ID;
    G4int fDriftHC1ID;
    G4int fDriftHC2ID;
    G4int fEmCalHCID;
    G4int fHadCalHCID;
    std::vector<G4double> fEmCalEdep;
    std::vector<G4double> fHadCalEdep;
};

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#endif

#include "B5EventAction.hh"
#include "B5HodoscopeHit.hh"
#include "B5DriftChamberHit.hh"
#include "B5EmCalorimeterHit.hh"
#include "B5HadCalorimeterHit.hh"
#include "B5Constants.hh"
#include "B5Analysis.hh"

#include "G4Event.hh"
#include "G4RunManager.hh"
#include "G4EventManager.hh"
#include "G4HCofThisEvent.hh"
#include "G4VHitsCollection.hh"
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4ios.hh"


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B5EventAction::B5EventAction()
: G4UserEventAction(), 
  fHodHC1ID(-1),
  fHodHC2ID(-1),
  fDriftHC1ID(-1),
  fDriftHC2ID(-1),
  fEmCalHCID(-1),
  fHadCalHCID(-1),
  fEmCalEdep(kNofEmCells, 0.), 
  fHadCalEdep(kNofHadCells, 0.)
{
  // set printing per each event
  G4RunManager::GetRunManager()->SetPrintProgress(1);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

B5EventAction::~B5EventAction()
{}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B5EventAction::BeginOfEventAction(const G4Event*)
{
  if (fHodHC1ID==-1) {
    auto sdManager = G4SDManager::GetSDMpointer();
    fHodHC1ID = sdManager->GetCollectionID("hodoscope1/hodoscopeColl");
    fHodHC2ID = sdManager->GetCollectionID("hodoscope2/hodoscopeColl");
    fDriftHC1ID = sdManager->GetCollectionID("chamber1/driftChamberColl");
    fDriftHC2ID = sdManager->GetCollectionID("chamber2/driftChamberColl");
    fEmCalHCID = sdManager->GetCollectionID("EMcalorimeter/EMcalorimeterColl");
    fHadCalHCID = sdManager->GetCollectionID("HadCalorimeter/HadCalorimeterColl");
  }
}     

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void B5EventAction::EndOfEventAction(const G4Event* event)
{
  auto hce = event->GetHCofThisEvent();
  if (!hce) {
      G4ExceptionDescription msg;
      msg << "No hits collection of this event found." << G4endl; 
      G4Exception("B5EventAction::EndOfEventAction()",
                  "B5Code001", JustWarning, msg);
      return;
  }

  // Get hits collections 
  auto hHC1 
    = static_cast<B5HodoscopeHitsCollection*>(hce->GetHC(fHodHC1ID));
    
  auto hHC2 
    = static_cast<B5HodoscopeHitsCollection*>(hce->GetHC(fHodHC2ID));
    
  auto dHC1 
    = static_cast<B5DriftChamberHitsCollection*>(hce->GetHC(fDriftHC1ID));
    
  auto dHC2 
    = static_cast<B5DriftChamberHitsCollection*>(hce->GetHC(fDriftHC2ID));
    
  auto ecHC 
    = static_cast<B5EmCalorimeterHitsCollection*>(hce->GetHC(fEmCalHCID));
    
  auto hcHC 
    = static_cast<B5HadCalorimeterHitsCollection*>(hce->GetHC(fHadCalHCID));
    
  if ( (!hHC1) || (!hHC2) || (!dHC1) || (!dHC2) || (!ecHC) || (!hcHC) ) {
      G4ExceptionDescription msg;
      msg << "Some of hits collections of this event not found." << G4endl; 
      G4Exception("B5EventAction::EndOfEventAction()",
                  "B5Code001", JustWarning, msg);
      return;
  }   
  
  //
  // Fill histograms & ntuple
  // 
  
  // Get analysis manager
  auto analysisManager = G4AnalysisManager::Instance();
 
  // Fill histograms
 
  auto nhit  = dHC1->entries();
  analysisManager->FillH1(0, nhit );

  for (auto i=0;i<nhit ;i++) {
    auto hit = (*dHC1)[i];
    auto localPos = hit->GetLocalPos();
    analysisManager->FillH2(0, localPos.x(), localPos.y());
  }
 
  nhit  = dHC2->entries();
  analysisManager->FillH1(1, nhit );

  for (auto i=0;i<nhit ;i++) {
    auto hit = (*dHC2)[i];
    auto localPos = hit->GetLocalPos();
    analysisManager->FillH2(1, localPos.x(), localPos.y());
  }
      
  // Fill ntuple
  
  // Dc1Hits
  analysisManager->FillNtupleIColumn(0, dHC1->entries());
  // Dc2Hits
  analysisManager->FillNtupleIColumn(1, dHC1->entries());
  
  // ECEnergy
  G4int totalEmHit = 0;
  G4double totalEmE = 0.;
  for (auto i=0;i<kNofEmCells;i++) {
    auto hit = (*ecHC)[i];
    auto edep = hit->GetEdep();
    if (edep>0.) {
      totalEmHit++;
      totalEmE += edep;
    }
    fEmCalEdep[i] = edep;
  }
  analysisManager->FillNtupleDColumn(2, totalEmE);

  // HCEnergy
  G4int totalHadHit = 0;
  G4double totalHadE = 0.;
  for (auto i=0;i<kNofHadCells;i++) {
    auto hit = (*hcHC)[i];
    auto edep = hit->GetEdep();
    if (edep>0.) {
        totalHadHit++;
        totalHadE += edep;
    }
    fHadCalEdep[i] = edep;
  }
  analysisManager->FillNtupleDColumn(3, totalHadE);

  // Time 1
  for (auto i=0;i<hHC1->entries();i++) {
    analysisManager->FillNtupleDColumn(4,(*hHC1)[i]->GetTime());
  }
    
  // Time 2
  for (auto i=0;i<hHC2->entries();i++) {
    analysisManager->FillNtupleDColumn(5,(*hHC2)[i]->GetTime());
  }
  
  analysisManager->AddNtupleRow();  
  
  //
  // Print diagnostics
  // 
  
  auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
  if ( printModulo==0 || event->GetEventID() % printModulo != 0) return;
  
  auto primary = event->GetPrimaryVertex(0)->GetPrimary(0);
  G4cout 
    << G4endl
    << ">>> Event " << event->GetEventID() << " >>> Simulation truth : "
    << primary->GetG4code()->GetParticleName()
    << " " << primary->GetMomentum() << G4endl;
  
  // Hodoscope 1
  nhit  = hHC1->entries();
  G4cout << "Hodoscope 1 has " << nhit  << " hits." << G4endl;
  for (auto i=0;i<nhit ;i++) {
    auto hit = (*hHC1)[i];
    hit->Print();
  }

  // Hodoscope 2
  nhit  = hHC2->entries();
  G4cout << "Hodoscope 2 has " << nhit  << " hits." << G4endl;
  for (auto i=0;i<nhit ;i++) {
    auto hit = (*hHC2)[i];
    hit->Print();
  }

  // Drift chamber 1
  nhit  = dHC1->entries();
  G4cout << "Drift Chamber 1 has " << nhit  << " hits." << G4endl;
  for (auto layer=0;layer<kNofChambers;layer++) {
    for (auto i=0;i<nhit ;i++) {
      auto hit = (*dHC1)[i];
      if (hit->GetLayerID()==layer) hit->Print();
    }
  }

  // Drift chamber 2
  nhit  = dHC2->entries();
  G4cout << "Drift Chamber 2 has " << nhit  << " hits." << G4endl;
  for (auto layer=0;layer<kNofChambers;layer++) {
    for (auto i=0;i<nhit ;i++) {
      auto hit = (*dHC2)[i];
      if (hit->GetLayerID()==layer) hit->Print();
    }
  }

  // EM calorimeter
  G4cout << "EM Calorimeter has " << totalEmHit << " hits. Total Edep is "
    << totalEmE/MeV << " (MeV)" << G4endl;

  // Had calorimeter
  G4cout << "Hadron Calorimeter has " << totalHadHit << " hits. Total Edep is "
    << totalHadE/MeV << " (MeV)" << G4endl;
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

• 第 47 - 52 行，是为了得到不同的灵敏探测器所产生的 hitCollection 在 G4HCofThisEvent 容器中的位置，使用灵敏探测器管理器类的成员函数  GetCollectionID(G4String colName) 返回对应的 hitCollection 的 ID，其参数为 SDname/hitCollectionName，中间需要加"/"，入示例中所示。
• 在  EndOfEventAction() 中，首先需要定义一个  GetHCofThisEvent 容器对象，通过 hitCollectionID 索引得到 hitCollection 的容器，在文档中的第 69 - 94 行。
• 第 130 - 156 行，这里用了累加，将每一个 hitCollection 容器中 hit 对象的能量相加，得到一个 Event 中，在一个灵敏探测器中沉积的总的能量。