| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include "Par04EventAction.hh"
27
28 #include "Par04DetectorConstruction.hh" // fo
29 #include "Par04Hit.hh" // for Par04Hit, Par04
30 #include "Par04ParallelFullWorld.hh"
31
32 #include "G4AnalysisManager.hh" // for G4Anal
33 #include "G4Event.hh" // for G4Event
34 #include "G4EventManager.hh" // for G4EventMa
35 #include "G4Exception.hh" // for G4Exception,
36 #include "G4ExceptionSeverity.hh" // for Fata
37 #include "G4GenericAnalysisManager.hh" // for
38 #include "G4HCofThisEvent.hh" // for G4HCofTh
39 #include "G4PrimaryParticle.hh" // for G4Prim
40 #include "G4PrimaryVertex.hh" // for G4Primar
41 #include "G4SDManager.hh" // for G4SDManager
42 #include "G4SystemOfUnits.hh" // for GeV
43 #include "G4THitsCollection.hh" // for G4THit
44 #include "G4ThreeVector.hh" // for G4ThreeVec
45 #include "G4Timer.hh" // for G4Timer
46 #include "G4UserEventAction.hh" // for G4User
47
48 #include <CLHEP/Units/SystemOfUnits.h> // for
49 #include <CLHEP/Vector/ThreeVector.h> // for
50 #include <algorithm> // for max
51 #include <cmath> // for log10
52 #include <cstddef> // for size_t
53 #include <ostream> // for basic_ostream::oper
54
55 //....oooOO0OOooo........oooOO0OOooo........oo
56
57 Par04EventAction::Par04EventAction(Par04Detect
58 Par04Parall
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 * f
71 fCalRho.reserve(fCellNbRho * fCellNbPhi * fC
72 fCalPhi.reserve(fCellNbRho * fCellNbPhi * fC
73 fCalZ.reserve(fCellNbRho * fCellNbPhi * fCel
74 fCalPhysicalEdep.reserve(fPhysicalNbLayers *
75 fCalPhysicalLayer.reserve(fPhysicalNbLayers
76 fCalPhysicalSlice.reserve(fPhysicalNbLayers
77 fCalPhysicalRow.reserve(fPhysicalNbLayers *
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81
82 Par04EventAction::~Par04EventAction() = defaul
83
84 //....oooOO0OOooo........oooOO0OOooo........oo
85
86 void Par04EventAction::BeginOfEventAction(cons
87 {
88 StartTimer();
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93 void Par04EventAction::StartTimer()
94 {
95 fTimer.Start();
96 }
97
98 //....oooOO0OOooo........oooOO0OOooo........oo
99
100 void Par04EventAction::StopTimer()
101 {
102 fTimer.Stop();
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oo
106
107 void Par04EventAction::EndOfEventAction(const
108 {
109 G4SDManager::GetSDMpointer()->GetHCtable();
110 StopTimer();
111
112 // Get hits collection ID (only once)
113 if (fHitCollectionID == -1) {
114 fHitCollectionID = G4SDManager::GetSDMpoin
115 }
116 if (fPhysicalFullHitCollectionID == -1) {
117 fPhysicalFullHitCollectionID =
118 G4SDManager::GetSDMpointer()->GetCollect
119 }
120 if (fPhysicalFastHitCollectionID == -1) {
121 fPhysicalFastHitCollectionID =
122 G4SDManager::GetSDMpointer()->GetCollect
123 }
124 // Get hits collection
125 auto hitsCollection =
126 static_cast<Par04HitsCollection*>(aEvent->
127 auto physicalFullHitsCollection = static_cas
128 aEvent->GetHCofThisEvent()->GetHC(fPhysica
129 auto physicalFastHitsCollection = static_cas
130 aEvent->GetHCofThisEvent()->GetHC(fPhysica
131
132 if (hitsCollection == nullptr) {
133 G4ExceptionDescription msg;
134 msg << "Cannot access hitsCollection ID "
135 G4Exception("Par04EventAction::GetHitsColl
136 }
137 if (physicalFullHitsCollection == nullptr) {
138 G4ExceptionDescription msg;
139 msg << "Cannot access physical full sim hi
140 G4Exception("Par04EventAction::GetHitsColl
141 }
142 if (physicalFastHitsCollection == nullptr) {
143 G4ExceptionDescription msg;
144 msg << "Cannot access physical fast sim hi
145 G4Exception("Par04EventAction::GetHitsColl
146 }
147 // Get analysis manager
148 auto analysisManager = G4AnalysisManager::In
149 // Retrieve only once detector dimensions
150 if (fCellSizeZ == 0) {
151 fCellSizeZ = fDetector->GetMeshSizeOfCells
152 fCellSizePhi = fDetector->GetMeshSizeOfCel
153 fCellSizeRho = fDetector->GetMeshSizeOfCel
154 fCellNbRho = fDetector->GetMeshNbOfCells()
155 fCellNbPhi = fDetector->GetMeshNbOfCells()
156 fCellNbZ = fDetector->GetMeshNbOfCells().z
157 }
158 if (fPhysicalNbLayers == 0) {
159 fPhysicalNbLayers = fParallel->GetNbOfLaye
160 fPhysicalNbSlices = fParallel->GetNbOfSlic
161 fPhysicalNbRows = fParallel->GetNbOfRows()
162 }
163
164 // Retrieve information from primary vertex
165 // To calculate shower axis and entry point
166 auto primaryVertex =
167 G4EventManager::GetEventManager()->GetCons
168 auto primaryParticle = primaryVertex->GetPri
169 G4double primaryEnergy = primaryParticle->Ge
170 // Estimate from vertex and particle directi
171 // Calculate entrance point to the detector
172 auto primaryDirection = primaryParticle->Get
173 auto primaryEntrance =
174 primaryVertex->GetPosition() - primaryVert
175
176 // Resize back to initial mesh size
177 fCalEdep.resize(fCellNbRho * fCellNbPhi * fC
178 fCalRho.resize(fCellNbRho * fCellNbPhi * fCe
179 fCalPhi.resize(fCellNbRho * fCellNbPhi * fCe
180 fCalZ.resize(fCellNbRho * fCellNbPhi * fCell
181 fCalPhysicalEdep.resize(fPhysicalNbLayers *
182 fCalPhysicalLayer.resize(fPhysicalNbLayers *
183 fCalPhysicalSlice.resize(fPhysicalNbLayers *
184 fCalPhysicalRow.resize(fPhysicalNbLayers * f
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., phi
198 G4double tFirstMoment = 0., tSecondMoment =
199 G4double rFirstMoment = 0., rSecondMoment =
200 G4double phiMean = 0.;
201 for (size_t iHit = 0; iHit < hitsCollection-
202 hit = static_cast<Par04Hit*>(hitsCollectio
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, hi
219 analysisManager->FillH1(5, rDistance, hi
220 analysisManager->FillH1(10, hitType);
221 if (hitEn > 0.0005) { // e > 0.5 keV
222 fCalEdep[numNonZeroThresholdCells] = h
223 fCalRho[numNonZeroThresholdCells] = hi
224 fCalPhi[numNonZeroThresholdCells] = hi
225 fCalZ[numNonZeroThresholdCells] = hitZ
226 numNonZeroThresholdCells++;
227 analysisManager->FillH1(13, std::log10
228 analysisManager->FillH1(15, hitNum);
229 }
230 }
231 }
232 tFirstMoment /= totalEnergy;
233 rFirstMoment /= totalEnergy;
234 phiMean /= totalEnergy;
235 analysisManager->FillH1(0, primaryEnergy / G
236 analysisManager->FillH1(1, totalEnergy / GeV
237 analysisManager->FillH1(2, totalEnergy / pri
238 analysisManager->FillH1(3, fTimer.GetRealEla
239 analysisManager->FillH1(6, tFirstMoment);
240 analysisManager->FillH1(7, rFirstMoment);
241 analysisManager->FillH1(12, numNonZeroThresh
242 analysisManager->FillH1(14, totalNum);
243 // Resize to store only energy hits above th
244 fCalEdep.resize(numNonZeroThresholdCells);
245 fCalRho.resize(numNonZeroThresholdCells);
246 fCalPhi.resize(numNonZeroThresholdCells);
247 fCalZ.resize(numNonZeroThresholdCells);
248 analysisManager->FillNtupleDColumn(0, 0, pri
249 analysisManager->FillNtupleDColumn(0, 1, fTi
250 // Second loop over hits to calculate second
251 for (size_t iHit = 0; iHit < hitsCollection-
252 hit = static_cast<Par04Hit*>(hitsCollectio
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(tDista
262 rSecondMoment += hitEn * std::pow(rDista
263 analysisManager->FillH1(11, phiDistance
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 < physicalFullHit
281 hit = static_cast<Par04Hit*>(physicalFullH
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[numNonZeroThresholdCe
292 fCalPhysicalLayer[numNonZeroThresholdC
293 fCalPhysicalRow[numNonZeroThresholdCel
294 fCalPhysicalSlice[numNonZeroThresholdC
295 numNonZeroThresholdCells++;
296 analysisManager->FillH1(19, std::log10
297 analysisManager->FillH1(21, hitNum);
298 }
299 }
300 }
301 for (size_t iHit = 0; iHit < physicalFastHit
302 hit = static_cast<Par04Hit*>(physicalFastH
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[numNonZeroThresholdCe
313 fCalPhysicalLayer[numNonZeroThresholdC
314 fCalPhysicalRow[numNonZeroThresholdCel
315 fCalPhysicalSlice[numNonZeroThresholdC
316 numNonZeroThresholdCells++;
317 analysisManager->FillH1(19, std::log10
318 analysisManager->FillH1(21, hitNum);
319 }
320 }
321 }
322 analysisManager->FillH1(16, totalPhysicalEne
323 analysisManager->FillH1(17, totalPhysicalEne
324 analysisManager->FillH1(18, numNonZeroThresh
325 analysisManager->FillH1(20, totalNum);
326 fCalPhysicalEdep.resize(numNonZeroThresholdC
327 fCalPhysicalLayer.resize(numNonZeroThreshold
328 fCalPhysicalSlice.resize(numNonZeroThreshold
329 fCalPhysicalRow.resize(numNonZeroThresholdCe
330 analysisManager->AddNtupleRow(0);
331 analysisManager->AddNtupleRow(1);
332 analysisManager->AddNtupleRow(2);
333 }
334