| 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 "Par04EventAction.hh"
27
28 #include "Par04DetectorConstruction.hh" // for Par04DetectorConstruction
29 #include "Par04Hit.hh" // for Par04Hit, Par04HitsCollection
30 #include "Par04ParallelFullWorld.hh"
31
32 #include "G4AnalysisManager.hh" // for G4AnalysisManager
33 #include "G4Event.hh" // for G4Event
34 #include "G4EventManager.hh" // for G4EventManager
35 #include "G4Exception.hh" // for G4Exception, G4ExceptionDesc...
36 #include "G4ExceptionSeverity.hh" // for FatalException
37 #include "G4GenericAnalysisManager.hh" // for G4GenericAnalysisManager
38 #include "G4HCofThisEvent.hh" // for G4HCofThisEvent
39 #include "G4PrimaryParticle.hh" // for G4PrimaryParticle
40 #include "G4PrimaryVertex.hh" // for G4PrimaryVertex
41 #include "G4SDManager.hh" // for G4SDManager
42 #include "G4SystemOfUnits.hh" // for GeV
43 #include "G4THitsCollection.hh" // for G4THitsCollection
44 #include "G4ThreeVector.hh" // for G4ThreeVector
45 #include "G4Timer.hh" // for G4Timer
46 #include "G4UserEventAction.hh" // for G4UserEventAction
47
48 #include <CLHEP/Units/SystemOfUnits.h> // for GeV
49 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
50 #include <algorithm> // for max
51 #include <cmath> // for log10
52 #include <cstddef> // for size_t
53 #include <ostream> // for basic_ostream::operator<<
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57 Par04EventAction::Par04EventAction(Par04DetectorConstruction* aDetector,
58 Par04ParallelFullWorld* aParallel)
59 : G4UserEventAction(),
60 fHitCollectionID(-1),
61 fPhysicalFullHitCollectionID(-1),
62 fPhysicalFastHitCollectionID(-1),
63 fTimer(),
64 fDetector(aDetector),
65 fParallel(aParallel)
66 {
67 fCellNbRho = aDetector->GetMeshNbOfCells().x();
68 fCellNbPhi = aDetector->GetMeshNbOfCells().y();
69 fCellNbZ = aDetector->GetMeshNbOfCells().z();
70 fCalEdep.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
71 fCalRho.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
72 fCalPhi.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
73 fCalZ.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
74 fCalPhysicalEdep.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
75 fCalPhysicalLayer.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
76 fCalPhysicalSlice.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
77 fCalPhysicalRow.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 Par04EventAction::~Par04EventAction() = default;
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85
86 void Par04EventAction::BeginOfEventAction(const G4Event*)
87 {
88 StartTimer();
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92
93 void Par04EventAction::StartTimer()
94 {
95 fTimer.Start();
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99
100 void Par04EventAction::StopTimer()
101 {
102 fTimer.Stop();
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
106
107 void Par04EventAction::EndOfEventAction(const G4Event* aEvent)
108 {
109 G4SDManager::GetSDMpointer()->GetHCtable();
110 StopTimer();
111
112 // Get hits collection ID (only once)
113 if (fHitCollectionID == -1) {
114 fHitCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("hits");
115 }
116 if (fPhysicalFullHitCollectionID == -1) {
117 fPhysicalFullHitCollectionID =
118 G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFullSim");
119 }
120 if (fPhysicalFastHitCollectionID == -1) {
121 fPhysicalFastHitCollectionID =
122 G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFastSim");
123 }
124 // Get hits collection
125 auto hitsCollection =
126 static_cast<Par04HitsCollection*>(aEvent->GetHCofThisEvent()->GetHC(fHitCollectionID));
127 auto physicalFullHitsCollection = static_cast<Par04HitsCollection*>(
128 aEvent->GetHCofThisEvent()->GetHC(fPhysicalFullHitCollectionID));
129 auto physicalFastHitsCollection = static_cast<Par04HitsCollection*>(
130 aEvent->GetHCofThisEvent()->GetHC(fPhysicalFastHitCollectionID));
131
132 if (hitsCollection == nullptr) {
133 G4ExceptionDescription msg;
134 msg << "Cannot access hitsCollection ID " << fHitCollectionID;
135 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
136 }
137 if (physicalFullHitsCollection == nullptr) {
138 G4ExceptionDescription msg;
139 msg << "Cannot access physical full sim hitsCollection ID " << fPhysicalFullHitCollectionID;
140 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
141 }
142 if (physicalFastHitsCollection == nullptr) {
143 G4ExceptionDescription msg;
144 msg << "Cannot access physical fast sim hitsCollection ID " << fPhysicalFastHitCollectionID;
145 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
146 }
147 // Get analysis manager
148 auto analysisManager = G4AnalysisManager::Instance();
149 // Retrieve only once detector dimensions
150 if (fCellSizeZ == 0) {
151 fCellSizeZ = fDetector->GetMeshSizeOfCells().z();
152 fCellSizePhi = fDetector->GetMeshSizeOfCells().y();
153 fCellSizeRho = fDetector->GetMeshSizeOfCells().x();
154 fCellNbRho = fDetector->GetMeshNbOfCells().x();
155 fCellNbPhi = fDetector->GetMeshNbOfCells().y();
156 fCellNbZ = fDetector->GetMeshNbOfCells().z();
157 }
158 if (fPhysicalNbLayers == 0) {
159 fPhysicalNbLayers = fParallel->GetNbOfLayers();
160 fPhysicalNbSlices = fParallel->GetNbOfSlices();
161 fPhysicalNbRows = fParallel->GetNbOfRows();
162 }
163
164 // Retrieve information from primary vertex and primary particle
165 // To calculate shower axis and entry point to the detector
166 auto primaryVertex =
167 G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetPrimaryVertex();
168 auto primaryParticle = primaryVertex->GetPrimary(0);
169 G4double primaryEnergy = primaryParticle->GetTotalEnergy();
170 // Estimate from vertex and particle direction the entry point to the detector
171 // Calculate entrance point to the detector located at z = 0
172 auto primaryDirection = primaryParticle->GetMomentumDirection();
173 auto primaryEntrance =
174 primaryVertex->GetPosition() - primaryVertex->GetPosition().z() * primaryDirection;
175
176 // Resize back to initial mesh size
177 fCalEdep.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
178 fCalRho.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
179 fCalPhi.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
180 fCalZ.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
181 fCalPhysicalEdep.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
182 fCalPhysicalLayer.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
183 fCalPhysicalSlice.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
184 fCalPhysicalRow.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
185
186 // Fill histograms
187 Par04Hit* hit = nullptr;
188 G4double hitEn = 0;
189 G4double totalEnergy = 0;
190 G4int hitNum = 0;
191 G4int totalNum = 0;
192 G4int hitZ = -1;
193 G4int hitRho = -1;
194 G4int hitPhi = -1;
195 G4int hitType = -1;
196 G4int numNonZeroThresholdCells = 0;
197 G4double tDistance = 0., rDistance = 0., phiDistance = 0.;
198 G4double tFirstMoment = 0., tSecondMoment = 0.;
199 G4double rFirstMoment = 0., rSecondMoment = 0.;
200 G4double phiMean = 0.;
201 for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
202 hit = static_cast<Par04Hit*>(hitsCollection->GetHit(iHit));
203 hitZ = hit->GetZid();
204 hitRho = hit->GetRhoId();
205 hitPhi = hit->GetPhiId();
206 hitEn = hit->GetEdep();
207 hitNum = hit->GetNdep();
208 hitType = hit->GetType();
209 if (hitEn > 0) {
210 totalEnergy += hitEn;
211 totalNum += hitNum;
212 tDistance = hitZ * fCellSizeZ;
213 rDistance = hitRho * fCellSizeRho;
214 phiDistance = hitPhi * fCellSizePhi;
215 tFirstMoment += hitEn * tDistance;
216 rFirstMoment += hitEn * rDistance;
217 phiMean += hitEn * phiDistance;
218 analysisManager->FillH1(4, tDistance, hitEn);
219 analysisManager->FillH1(5, rDistance, hitEn);
220 analysisManager->FillH1(10, hitType);
221 if (hitEn > 0.0005) { // e > 0.5 keV
222 fCalEdep[numNonZeroThresholdCells] = hitEn;
223 fCalRho[numNonZeroThresholdCells] = hitRho;
224 fCalPhi[numNonZeroThresholdCells] = hitPhi;
225 fCalZ[numNonZeroThresholdCells] = hitZ;
226 numNonZeroThresholdCells++;
227 analysisManager->FillH1(13, std::log10(hitEn));
228 analysisManager->FillH1(15, hitNum);
229 }
230 }
231 }
232 tFirstMoment /= totalEnergy;
233 rFirstMoment /= totalEnergy;
234 phiMean /= totalEnergy;
235 analysisManager->FillH1(0, primaryEnergy / GeV);
236 analysisManager->FillH1(1, totalEnergy / GeV);
237 analysisManager->FillH1(2, totalEnergy / primaryEnergy);
238 analysisManager->FillH1(3, fTimer.GetRealElapsed());
239 analysisManager->FillH1(6, tFirstMoment);
240 analysisManager->FillH1(7, rFirstMoment);
241 analysisManager->FillH1(12, numNonZeroThresholdCells);
242 analysisManager->FillH1(14, totalNum);
243 // Resize to store only energy hits above threshold
244 fCalEdep.resize(numNonZeroThresholdCells);
245 fCalRho.resize(numNonZeroThresholdCells);
246 fCalPhi.resize(numNonZeroThresholdCells);
247 fCalZ.resize(numNonZeroThresholdCells);
248 analysisManager->FillNtupleDColumn(0, 0, primaryEnergy);
249 analysisManager->FillNtupleDColumn(0, 1, fTimer.GetRealElapsed());
250 // Second loop over hits to calculate second moments
251 for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
252 hit = static_cast<Par04Hit*>(hitsCollection->GetHit(iHit));
253 hitEn = hit->GetEdep();
254 hitZ = hit->GetZid();
255 hitRho = hit->GetRhoId();
256 hitPhi = hit->GetPhiId();
257 if (hitEn > 0) {
258 tDistance = hitZ * fCellSizeZ;
259 rDistance = hitRho * fCellSizeRho;
260 phiDistance = hitPhi * fCellSizePhi;
261 tSecondMoment += hitEn * std::pow(tDistance - tFirstMoment, 2);
262 rSecondMoment += hitEn * std::pow(rDistance - rFirstMoment, 2);
263 analysisManager->FillH1(11, phiDistance - phiMean, hitEn);
264 }
265 }
266 tSecondMoment /= totalEnergy;
267 rSecondMoment /= totalEnergy;
268 analysisManager->FillH1(8, tSecondMoment);
269 analysisManager->FillH1(9, rSecondMoment);
270
271 // Fill ntuple with physical readout data
272 G4double totalPhysicalEnergy = 0;
273 totalNum = 0;
274 hitEn = 0;
275 hitNum = 0;
276 G4int hitLayer = -1;
277 G4int hitRow = -1;
278 G4int hitSlice = -1;
279 numNonZeroThresholdCells = 0;
280 for (size_t iHit = 0; iHit < physicalFullHitsCollection->entries(); iHit++) {
281 hit = static_cast<Par04Hit*>(physicalFullHitsCollection->GetHit(iHit));
282 hitLayer = hit->GetRhoId();
283 hitRow = hit->GetZid();
284 hitSlice = hit->GetPhiId();
285 hitEn = hit->GetEdep();
286 hitNum = hit->GetNdep();
287 if (hitEn > 0) {
288 totalPhysicalEnergy += hitEn;
289 totalNum += hitNum;
290 if (hitEn > 0.0005) { // e > 0.5 keV
291 fCalPhysicalEdep[numNonZeroThresholdCells] = hitEn;
292 fCalPhysicalLayer[numNonZeroThresholdCells] = hitLayer;
293 fCalPhysicalRow[numNonZeroThresholdCells] = hitRow;
294 fCalPhysicalSlice[numNonZeroThresholdCells] = hitSlice;
295 numNonZeroThresholdCells++;
296 analysisManager->FillH1(19, std::log10(hitEn));
297 analysisManager->FillH1(21, hitNum);
298 }
299 }
300 }
301 for (size_t iHit = 0; iHit < physicalFastHitsCollection->entries(); iHit++) {
302 hit = static_cast<Par04Hit*>(physicalFastHitsCollection->GetHit(iHit));
303 hitLayer = hit->GetRhoId();
304 hitRow = hit->GetZid();
305 hitSlice = hit->GetPhiId();
306 hitEn = hit->GetEdep();
307 hitNum = hit->GetNdep();
308 if (hitEn > 0) {
309 totalPhysicalEnergy += hitEn;
310 totalNum += hitNum;
311 if (hitEn > 0.0005) { // e > 0.5 keV
312 fCalPhysicalEdep[numNonZeroThresholdCells] = hitEn;
313 fCalPhysicalLayer[numNonZeroThresholdCells] = hitLayer;
314 fCalPhysicalRow[numNonZeroThresholdCells] = hitRow;
315 fCalPhysicalSlice[numNonZeroThresholdCells] = hitSlice;
316 numNonZeroThresholdCells++;
317 analysisManager->FillH1(19, std::log10(hitEn));
318 analysisManager->FillH1(21, hitNum);
319 }
320 }
321 }
322 analysisManager->FillH1(16, totalPhysicalEnergy / GeV);
323 analysisManager->FillH1(17, totalPhysicalEnergy / primaryEnergy);
324 analysisManager->FillH1(18, numNonZeroThresholdCells);
325 analysisManager->FillH1(20, totalNum);
326 fCalPhysicalEdep.resize(numNonZeroThresholdCells);
327 fCalPhysicalLayer.resize(numNonZeroThresholdCells);
328 fCalPhysicalSlice.resize(numNonZeroThresholdCells);
329 fCalPhysicalRow.resize(numNonZeroThresholdCells);
330 analysisManager->AddNtupleRow(0);
331 analysisManager->AddNtupleRow(1);
332 analysisManager->AddNtupleRow(2);
333 }
334