| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 #include "Par04SensitiveDetector.hh"
27
28 #include "Par04EventInformation.hh" // for Par04EventInformation
29 #include "Par04Hit.hh" // for Par04Hit, Par04HitsCollection
30
31 #include "G4Event.hh" // for G4Event
32 #include "G4EventManager.hh" // for G4EventManager
33 #include "G4HCofThisEvent.hh" // for G4HCofThisEvent
34 #include "G4SDManager.hh" // for G4SDManager
35 #include "G4Step.hh" // for G4Step
36 #include "G4Track.hh" // for G4Track
37
38 #include <CLHEP/Vector/Rotation.h> // for HepRotation
39 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
40 #include <G4CollectionNameVector.hh> // for G4CollectionNameVector
41 #include <G4FastHit.hh> // for G4FastHit
42 #include <G4FastTrack.hh> // for G4FastTrack
43 #include <G4RotationMatrix.hh> // for G4RotationMatrix
44 #include <G4StepPoint.hh> // for G4StepPoint
45 #include <G4THitsCollection.hh> // for G4THitsCollection
46 #include <G4ThreeVector.hh> // for G4ThreeVector
47 #include <G4VSensitiveDetector.hh> // for G4VSensitiveDetector
48 #include <G4VUserEventInformation.hh> // for G4VUserEventInformation
49 #include <cmath> // for floor
50 #include <cstddef> // for size_t
51 #include <vector> // for vector
52
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54
55 Par04SensitiveDetector::Par04SensitiveDetector(G4String aName) : G4VSensitiveDetector(aName)
56 {
57 collectionName.insert("hits");
58 }
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61 Par04SensitiveDetector::Par04SensitiveDetector(G4String aName, G4ThreeVector aNb,
62 G4ThreeVector aSize)
63 : G4VSensitiveDetector(aName), fMeshNbOfCells(aNb), fMeshSizeOfCells(aSize)
64 {
65 collectionName.insert("hits");
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70 Par04SensitiveDetector::~Par04SensitiveDetector() = default;
71
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74 void Par04SensitiveDetector::Initialize(G4HCofThisEvent* aHCE)
75 {
76 fHitsCollection = new Par04HitsCollection(SensitiveDetectorName, collectionName[0]);
77 if (fHitCollectionID < 0) {
78 fHitCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID(fHitsCollection);
79 }
80 aHCE->AddHitsCollection(fHitCollectionID, fHitsCollection);
81
82 // reset entrance position
83 fEntrancePosition.set(-1, -1, -1);
84 fEntranceDirection.set(-1, -1, -1);
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88
89 G4bool Par04SensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory*)
90 {
91 G4double edep = aStep->GetTotalEnergyDeposit();
92 if (edep == 0.) return true;
93
94 auto hit = RetrieveAndSetupHit(aStep->GetPostStepPoint()->GetPosition());
95 if (hit == nullptr) return true;
96
97 // Add energy deposit from G4Step
98 hit->AddEdep(edep);
99 // Increment the counter
100 hit->AddNdep(1);
101
102 // Fill time information from G4Step
103 // If it's already filled, choose hit with earliest global time
104 if (hit->GetTime() == -1 || hit->GetTime() > aStep->GetTrack()->GetGlobalTime())
105 hit->SetTime(aStep->GetTrack()->GetGlobalTime());
106
107 // Set hit type to full simulation (only if hit is not already marked as fast
108 // sim)
109 if (hit->GetType() != 1) hit->SetType(0);
110
111 return true;
112 }
113
114 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115
116 G4bool Par04SensitiveDetector::ProcessHits(const G4FastHit* aHit, const G4FastTrack* aTrack,
117 G4TouchableHistory*)
118 {
119 G4double edep = aHit->GetEnergy();
120 if (edep == 0.) return true;
121
122 auto hit = RetrieveAndSetupHit(aHit->GetPosition());
123 if (hit == nullptr) return true;
124
125 // Add energy deposit from G4FastHit
126 hit->AddEdep(edep);
127 // Increment the counter
128 hit->AddNdep(1);
129
130 // Fill time information from G4FastTrack
131 // If it's already filled, choose hit with earliest global time
132 if (hit->GetTime() == -1 || hit->GetTime() > aTrack->GetPrimaryTrack()->GetGlobalTime()) {
133 hit->SetTime(aTrack->GetPrimaryTrack()->GetGlobalTime());
134 }
135
136 // Set hit type to fast simulation (even if hit was already marked as full
137 // sim, overwrite it)
138 hit->SetType(1);
139
140 return true;
141 }
142
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144
145 Par04Hit* Par04SensitiveDetector::RetrieveAndSetupHit(G4ThreeVector aGlobalPosition)
146 {
147 if (fEntrancePosition.x() == -1) {
148 auto info = dynamic_cast<Par04EventInformation*>(
149 G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetUserInformation());
150 if (info == nullptr) return nullptr;
151 fEntrancePosition = info->GetPosition();
152 fEntranceDirection = info->GetDirection();
153 }
154
155 auto delta = aGlobalPosition - fEntrancePosition;
156
157 // Calculate rotation matrix along the particle momentum direction
158 // It will rotate the shower axes to match the incoming particle direction
159 G4RotationMatrix rotMatrix = G4RotationMatrix();
160 double particleTheta = fEntranceDirection.theta();
161 double particlePhi = fEntranceDirection.phi();
162 rotMatrix.rotateZ(-particlePhi);
163 rotMatrix.rotateY(-particleTheta);
164 G4RotationMatrix rotMatrixInv = CLHEP::inverseOf(rotMatrix);
165
166 delta = rotMatrix * delta;
167
168 G4int rhoNo = std::floor(delta.perp() / fMeshSizeOfCells.x());
169 G4int phiNo = std::floor((CLHEP::pi + delta.phi()) / fMeshSizeOfCells.y());
170 G4int zNo = std::floor(delta.z() / fMeshSizeOfCells.z());
171
172 std::size_t hitID =
173 fMeshNbOfCells.x() * fMeshNbOfCells.z() * phiNo + fMeshNbOfCells.z() * rhoNo + zNo;
174
175 if (zNo >= fMeshNbOfCells.z() || rhoNo >= fMeshNbOfCells.x() || zNo < 0) {
176 return nullptr;
177 }
178
179 auto hit = fHitsMap[hitID].get();
180 if (hit == nullptr) {
181 fHitsMap[hitID] = std::unique_ptr<Par04Hit>(new Par04Hit());
182 hit = fHitsMap[hitID].get();
183 hit->SetPhiId(phiNo);
184 hit->SetRhoId(rhoNo);
185 hit->SetZid(zNo);
186 hit->SetRot(rotMatrixInv);
187 hit->SetPos(fEntrancePosition
188 + rotMatrixInv * G4ThreeVector(0, 0, (zNo + 0.5) * fMeshSizeOfCells.z()));
189 }
190 return hit;
191 }
192
193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194
195 void Par04SensitiveDetector::EndOfEvent(G4HCofThisEvent*)
196 {
197 for (const auto& hits : fHitsMap) {
198 fHitsCollection->insert(new Par04Hit(*hits.second.get()));
199 }
200 fHitsMap.clear();
201 }
202