geant4 灵敏探测器

如何模拟探测器的响应？

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

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

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

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

 1 void B2aDetectorConstruction::ConstructSDandField()
 2 {
 3   // Sensitive detectors
 4 
 5   G4String trackerChamberSDname = "B2/TrackerChamberSD";
 6   B2TrackerSD* aTrackerSD = new B2TrackerSD(trackerChamberSDname,
 7                                             "TrackerHitsCollection");
 8   G4SDManager::GetSDMpointer()->AddNewDetector(aTrackerSD);
 9   // Setting aTrackerSD to all logical volumes with the same name 
10   // of "Chamber_LV".
11   SetSensitiveDetector("Chamber_LV", aTrackerSD, true);
12 
13   // Create global magnetic field messenger.
14   // Uniform magnetic field is then created automatically if
15   // the field value is not zero.
16   G4ThreeVector fieldValue = G4ThreeVector();
17   fMagFieldMessenger = new G4GlobalMagFieldMessenger(fieldValue);
18   fMagFieldMessenger->SetVerboseLevel(1);
19   
20   // Register the field messenger for deleting
21   G4AutoDelete::Register(fMagFieldMessenger);
22 }

第 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

 1 /// Tracker hit class
 2 ///
 3 /// It defines data members to store the trackID, chamberNb, energy deposit,
 4 /// and position of charged particles in a selected volume:
 5 /// - fTrackID, fChamberNB, fEdep, fPos
 6 
 7 class B2TrackerHit : public G4VHit
 8 {
 9   public:
10     B2TrackerHit();
11     B2TrackerHit(const B2TrackerHit&);
12     virtual ~B2TrackerHit();
13 
14     // operators
15     const B2TrackerHit& operator=(const B2TrackerHit&);
16     G4int operator==(const B2TrackerHit&) const;
17 
18     inline void* operator new(size_t);
19     inline void  operator delete(void*);
20 
21     // methods from base class
22     virtual void Draw();
23     virtual void Print();
24 
25     // Set methods
26     void SetTrackID  (G4int track)      { fTrackID = track; };
27     void SetChamberNb(G4int chamb)      { fChamberNb = chamb; };
28     void SetEdep     (G4double de)      { fEdep = de; };
29     void SetPos      (G4ThreeVector xyz){ fPos = xyz; };
30 
31     // Get methods
32     G4int GetTrackID() const     { return fTrackID; };
33     G4int GetChamberNb() const   { return fChamberNb; };
34     G4double GetEdep() const     { return fEdep; };
35     G4ThreeVector GetPos() const { return fPos; };
36 
37   private:
38 
39       G4int         fTrackID;
40       G4int         fChamberNb;
41       G4double      fEdep;
42       G4ThreeVector fPos;
43 };
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 
47 typedef G4THitsCollection<B2TrackerHit> B2TrackerHitsCollection;
48 
49 extern G4ThreadLocal G4Allocator<B2TrackerHit>* B2TrackerHitAllocator;
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52 
53 inline void* B2TrackerHit::operator new(size_t)
54 {
55   if(!B2TrackerHitAllocator)
56       B2TrackerHitAllocator = new G4Allocator<B2TrackerHit>;
57   return (void *) B2TrackerHitAllocator->MallocSingle();
58 }
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
62 inline void B2TrackerHit::operator delete(void *hit)
63 {
64   B2TrackerHitAllocator->FreeSingle((B2TrackerHit*) hit);
65 }

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

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

 1 #ifndef B2TrackerSD_h
 2 #define B2TrackerSD_h 1
 3 
 4 #include "G4VSensitiveDetector.hh"
 5 
 6 #include "B2TrackerHit.hh"
 7 
 8 #include <vector>
 9 
10 class G4Step;
11 class G4HCofThisEvent;
12 
13 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
14 
15 /// B2Tracker sensitive detector class
16 ///
17 /// The hits are accounted in hits in ProcessHits() function which is called
18 /// by Geant4 kernel at each step. A hit is created with each step with non zero 
19 /// energy deposit.
20 
21 class B2TrackerSD : public G4VSensitiveDetector
22 {
23   public:
24     B2TrackerSD(const G4String& name, 
25                 const G4String& hitsCollectionName);
26     virtual ~B2TrackerSD();
27   
28     // methods from base class
29     virtual void   Initialize(G4HCofThisEvent* hitCollection);
30     virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* history);
31     virtual void   EndOfEvent(G4HCofThisEvent* hitCollection);
32 
33   private:
34     B2TrackerHitsCollection* fHitsCollection;
35 };
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
38 
39 #endif

 1 #include "B2TrackerSD.hh"
 2 #include "G4HCofThisEvent.hh"
 3 #include "G4Step.hh"
 4 #include "G4ThreeVector.hh"
 5 #include "G4SDManager.hh"
 6 #include "G4ios.hh"
 7 
 8 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 9 
10 B2TrackerSD::B2TrackerSD(const G4String& name,
11                          const G4String& hitsCollectionName) 
12  : G4VSensitiveDetector(name),
13    fHitsCollection(NULL)
14 {
15   collectionName.insert(hitsCollectionName);
16 }
17 
18 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
19 
20 B2TrackerSD::~B2TrackerSD() 
21 {}
22 
23 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
24 
25 void B2TrackerSD::Initialize(G4HCofThisEvent* hce)
26 {
27   // Create hits collection
28 
29   fHitsCollection 
30     = new B2TrackerHitsCollection(SensitiveDetectorName, collectionName[0]); 
31 
32   // Add this collection in hce
33 
34   G4int hcID 
35     = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
36   hce->AddHitsCollection( hcID, fHitsCollection ); 
37 }
38 
39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
40 
41 G4bool B2TrackerSD::ProcessHits(G4Step* aStep, 
42                                      G4TouchableHistory*)
43 {  
44   // energy deposit
45   G4double edep = aStep->GetTotalEnergyDeposit();
46 
47   if (edep==0.) return false;
48 
49   B2TrackerHit* newHit = new B2TrackerHit();
50 
51   newHit->SetTrackID  (aStep->GetTrack()->GetTrackID());
52   newHit->SetChamberNb(aStep->GetPreStepPoint()->GetTouchableHandle()
53                                                ->GetCopyNumber());
54   newHit->SetEdep(edep);
55   newHit->SetPos (aStep->GetPostStepPoint()->GetPosition());
56 
57   fHitsCollection->insert( newHit );
58 
59   //newHit->Print();
60 
61   return true;
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
66 void B2TrackerSD::EndOfEvent(G4HCofThisEvent*)
67 {
68   if ( verboseLevel>1 ) { 
69      G4int nofHits = fHitsCollection->entries();
70      G4cout << G4endl
71             << "-------->Hits Collection: in this event they are " << nofHits 
72             << " hits in the tracker chambers: " << G4endl;
73      for ( G4int i=0; i<nofHits; i++ ) (*fHitsCollection)[i]->Print();
74   }
75 }
76 
77 //....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）

 1 void B5DetectorConstruction::ConstructSDandField()
 2 {
 3   // sensitive detectors -----------------------------------------------------
 4   auto sdManager = G4SDManager::GetSDMpointer();
 5   G4String SDname;
 6   
 7   auto hodoscope1 = new B5HodoscopeSD(SDname="/hodoscope1");
 8   sdManager->AddNewDetector(hodoscope1);
 9   fHodoscope1Logical->SetSensitiveDetector(hodoscope1);
10 
11   auto hodoscope2 = new B5HodoscopeSD(SDname="/hodoscope2");
12   sdManager->AddNewDetector(hodoscope2);
13   fHodoscope2Logical->SetSensitiveDetector(hodoscope2);
14   
15   auto chamber1 = new B5DriftChamberSD(SDname="/chamber1");
16   sdManager->AddNewDetector(chamber1);
17   fWirePlane1Logical->SetSensitiveDetector(chamber1);
18 
19   auto chamber2 = new B5DriftChamberSD(SDname="/chamber2");
20   sdManager->AddNewDetector(chamber2);
21   fWirePlane2Logical->SetSensitiveDetector(chamber2);
22   
23   auto emCalorimeter = new B5EmCalorimeterSD(SDname="/EMcalorimeter");
24   sdManager->AddNewDetector(emCalorimeter);
25   fCellLogical->SetSensitiveDetector(emCalorimeter);
26   
27   auto hadCalorimeter = new B5HadCalorimeterSD(SDname="/HadCalorimeter");
28   sdManager->AddNewDetector(hadCalorimeter);
29   fHadCalScintiLogical->SetSensitiveDetector(hadCalorimeter);
30 
31   // magnetic field ----------------------------------------------------------
32   fMagneticField = new B5MagneticField();
33   fFieldMgr = new G4FieldManager();
34   fFieldMgr->SetDetectorField(fMagneticField);
35   fFieldMgr->CreateChordFinder(fMagneticField);
36   G4bool forceToAllDaughters = true;
37   fMagneticLogical->SetFieldManager(fFieldMgr, forceToAllDaughters);
38 
39   // Register the field and its manager for deleting
40   G4AutoDelete::Register(fMagneticField);
41   G4AutoDelete::Register(fFieldMgr);
42 }    
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

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

 1 #ifndef B5EventAction_h
 2 #define B5EventAction_h 1
 3 
 4 
 5 #include "G4UserEventAction.hh"
 6 #include "globals.hh"
 7 
 8 #include <vector>
 9 
10 /// Event action
11 
12 class B5EventAction : public G4UserEventAction
13 {
14 public:
15     B5EventAction();
16     virtual ~B5EventAction();
17     
18     virtual void BeginOfEventAction(const G4Event*);
19     virtual void EndOfEventAction(const G4Event*);
20 
21     std::vector<G4double>& GetEmCalEdep() { return fEmCalEdep; }
22     std::vector<G4double>& GetHadCalEdep() { return fEmCalEdep; }
23     
24 private:
25     G4int fHodHC1ID;
26     G4int fHodHC2ID;
27     G4int fDriftHC1ID;
28     G4int fDriftHC2ID;
29     G4int fEmCalHCID;
30     G4int fHadCalHCID;
31     std::vector<G4double> fEmCalEdep;
32     std::vector<G4double> fHadCalEdep;
33 };
34 
35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
36 
37 #endif

  1 #include "B5EventAction.hh"
  2 #include "B5HodoscopeHit.hh"
  3 #include "B5DriftChamberHit.hh"
  4 #include "B5EmCalorimeterHit.hh"
  5 #include "B5HadCalorimeterHit.hh"
  6 #include "B5Constants.hh"
  7 #include "B5Analysis.hh"
  8 
  9 #include "G4Event.hh"
 10 #include "G4RunManager.hh"
 11 #include "G4EventManager.hh"
 12 #include "G4HCofThisEvent.hh"
 13 #include "G4VHitsCollection.hh"
 14 #include "G4SDManager.hh"
 15 #include "G4SystemOfUnits.hh"
 16 #include "G4ios.hh"
 17 
 18 
 19 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 20 
 21 B5EventAction::B5EventAction()
 22 : G4UserEventAction(), 
 23   fHodHC1ID(-1),
 24   fHodHC2ID(-1),
 25   fDriftHC1ID(-1),
 26   fDriftHC2ID(-1),
 27   fEmCalHCID(-1),
 28   fHadCalHCID(-1),
 29   fEmCalEdep(kNofEmCells, 0.), 
 30   fHadCalEdep(kNofHadCells, 0.)
 31 {
 32   // set printing per each event
 33   G4RunManager::GetRunManager()->SetPrintProgress(1);
 34 }
 35 
 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 37 
 38 B5EventAction::~B5EventAction()
 39 {}
 40 
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 42 
 43 void B5EventAction::BeginOfEventAction(const G4Event*)
 44 {
 45   if (fHodHC1ID==-1) {
 46     auto sdManager = G4SDManager::GetSDMpointer();
 47     fHodHC1ID = sdManager->GetCollectionID("hodoscope1/hodoscopeColl");
 48     fHodHC2ID = sdManager->GetCollectionID("hodoscope2/hodoscopeColl");
 49     fDriftHC1ID = sdManager->GetCollectionID("chamber1/driftChamberColl");
 50     fDriftHC2ID = sdManager->GetCollectionID("chamber2/driftChamberColl");
 51     fEmCalHCID = sdManager->GetCollectionID("EMcalorimeter/EMcalorimeterColl");
 52     fHadCalHCID = sdManager->GetCollectionID("HadCalorimeter/HadCalorimeterColl");
 53   }
 54 }     
 55 
 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 57 
 58 void B5EventAction::EndOfEventAction(const G4Event* event)
 59 {
 60   auto hce = event->GetHCofThisEvent();
 61   if (!hce) {
 62       G4ExceptionDescription msg;
 63       msg << "No hits collection of this event found." << G4endl; 
 64       G4Exception("B5EventAction::EndOfEventAction()",
 65                   "B5Code001", JustWarning, msg);
 66       return;
 67   }
 68 
 69   // Get hits collections 
 70   auto hHC1 
 71     = static_cast<B5HodoscopeHitsCollection*>(hce->GetHC(fHodHC1ID));
 72     
 73   auto hHC2 
 74     = static_cast<B5HodoscopeHitsCollection*>(hce->GetHC(fHodHC2ID));
 75     
 76   auto dHC1 
 77     = static_cast<B5DriftChamberHitsCollection*>(hce->GetHC(fDriftHC1ID));
 78     
 79   auto dHC2 
 80     = static_cast<B5DriftChamberHitsCollection*>(hce->GetHC(fDriftHC2ID));
 81     
 82   auto ecHC 
 83     = static_cast<B5EmCalorimeterHitsCollection*>(hce->GetHC(fEmCalHCID));
 84     
 85   auto hcHC 
 86     = static_cast<B5HadCalorimeterHitsCollection*>(hce->GetHC(fHadCalHCID));
 87     
 88   if ( (!hHC1) || (!hHC2) || (!dHC1) || (!dHC2) || (!ecHC) || (!hcHC) ) {
 89       G4ExceptionDescription msg;
 90       msg << "Some of hits collections of this event not found." << G4endl; 
 91       G4Exception("B5EventAction::EndOfEventAction()",
 92                   "B5Code001", JustWarning, msg);
 93       return;
 94   }   
 95   
 96   //
 97   // Fill histograms & ntuple
 98   // 
 99   
100   // Get analysis manager
101   auto analysisManager = G4AnalysisManager::Instance();
102  
103   // Fill histograms
104  
105   auto nhit  = dHC1->entries();
106   analysisManager->FillH1(0, nhit );
107 
108   for (auto i=0;i<nhit ;i++) {
109     auto hit = (*dHC1)[i];
110     auto localPos = hit->GetLocalPos();
111     analysisManager->FillH2(0, localPos.x(), localPos.y());
112   }
113  
114   nhit  = dHC2->entries();
115   analysisManager->FillH1(1, nhit );
116 
117   for (auto i=0;i<nhit ;i++) {
118     auto hit = (*dHC2)[i];
119     auto localPos = hit->GetLocalPos();
120     analysisManager->FillH2(1, localPos.x(), localPos.y());
121   }
122       
123   // Fill ntuple
124   
125   // Dc1Hits
126   analysisManager->FillNtupleIColumn(0, dHC1->entries());
127   // Dc2Hits
128   analysisManager->FillNtupleIColumn(1, dHC1->entries());
129   
130   // ECEnergy
131   G4int totalEmHit = 0;
132   G4double totalEmE = 0.;
133   for (auto i=0;i<kNofEmCells;i++) {
134     auto hit = (*ecHC)[i];
135     auto edep = hit->GetEdep();
136     if (edep>0.) {
137       totalEmHit++;
138       totalEmE += edep;
139     }
140     fEmCalEdep[i] = edep;
141   }
142   analysisManager->FillNtupleDColumn(2, totalEmE);
143 
144   // HCEnergy
145   G4int totalHadHit = 0;
146   G4double totalHadE = 0.;
147   for (auto i=0;i<kNofHadCells;i++) {
148     auto hit = (*hcHC)[i];
149     auto edep = hit->GetEdep();
150     if (edep>0.) {
151         totalHadHit++;
152         totalHadE += edep;
153     }
154     fHadCalEdep[i] = edep;
155   }
156   analysisManager->FillNtupleDColumn(3, totalHadE);
157 
158   // Time 1
159   for (auto i=0;i<hHC1->entries();i++) {
160     analysisManager->FillNtupleDColumn(4,(*hHC1)[i]->GetTime());
161   }
162     
163   // Time 2
164   for (auto i=0;i<hHC2->entries();i++) {
165     analysisManager->FillNtupleDColumn(5,(*hHC2)[i]->GetTime());
166   }
167   
168   analysisManager->AddNtupleRow();  
169   
170   //
171   // Print diagnostics
172   // 
173   
174   auto printModulo = G4RunManager::GetRunManager()->GetPrintProgress();
175   if ( printModulo==0 || event->GetEventID() % printModulo != 0) return;
176   
177   auto primary = event->GetPrimaryVertex(0)->GetPrimary(0);
178   G4cout 
179     << G4endl
180     << ">>> Event " << event->GetEventID() << " >>> Simulation truth : "
181     << primary->GetG4code()->GetParticleName()
182     << " " << primary->GetMomentum() << G4endl;
183   
184   // Hodoscope 1
185   nhit  = hHC1->entries();
186   G4cout << "Hodoscope 1 has " << nhit  << " hits." << G4endl;
187   for (auto i=0;i<nhit ;i++) {
188     auto hit = (*hHC1)[i];
189     hit->Print();
190   }
191 
192   // Hodoscope 2
193   nhit  = hHC2->entries();
194   G4cout << "Hodoscope 2 has " << nhit  << " hits." << G4endl;
195   for (auto i=0;i<nhit ;i++) {
196     auto hit = (*hHC2)[i];
197     hit->Print();
198   }
199 
200   // Drift chamber 1
201   nhit  = dHC1->entries();
202   G4cout << "Drift Chamber 1 has " << nhit  << " hits." << G4endl;
203   for (auto layer=0;layer<kNofChambers;layer++) {
204     for (auto i=0;i<nhit ;i++) {
205       auto hit = (*dHC1)[i];
206       if (hit->GetLayerID()==layer) hit->Print();
207     }
208   }
209 
210   // Drift chamber 2
211   nhit  = dHC2->entries();
212   G4cout << "Drift Chamber 2 has " << nhit  << " hits." << G4endl;
213   for (auto layer=0;layer<kNofChambers;layer++) {
214     for (auto i=0;i<nhit ;i++) {
215       auto hit = (*dHC2)[i];
216       if (hit->GetLayerID()==layer) hit->Print();
217     }
218   }
219 
220   // EM calorimeter
221   G4cout << "EM Calorimeter has " << totalEmHit << " hits. Total Edep is "
222     << totalEmE/MeV << " (MeV)" << G4endl;
223 
224   // Had calorimeter
225   G4cout << "Hadron Calorimeter has " << totalHadHit << " hits. Total Edep is "
226     << totalHadE/MeV << " (MeV)" << G4endl;
227 }
228 
229 //....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 中，在一个灵敏探测器中沉积的总的能量。