Hits collection is retrieved, analysed, and histograms are filled.
105{
106 G4SDManager::GetSDMpointer()->GetHCtable();
108
109
111 {
113 }
115 {
117 G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFullSim");
118 }
120 {
122 G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFastSim");
123 }
124
125 auto hitsCollection =
127 auto physicalFullHitsCollection =
130 auto physicalFastHitsCollection =
133
134 if(hitsCollection == nullptr)
135 {
136 G4ExceptionDescription msg;
138 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
139 }
140 if(physicalFullHitsCollection == nullptr)
141 {
142 G4ExceptionDescription msg;
144 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
145 }
146 if(physicalFastHitsCollection == nullptr)
147 {
148 G4ExceptionDescription msg;
150 G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
151 }
152
153 auto analysisManager = G4AnalysisManager::Instance();
154
156 {
163 }
165 {
169 }
170
171
172
173 auto primaryVertex =
174 G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetPrimaryVertex();
175 auto primaryParticle = primaryVertex->GetPrimary(0);
176 G4double primaryEnergy = primaryParticle->GetTotalEnergy();
177
178
179 auto primaryDirection = primaryParticle->GetMomentumDirection();
180 auto primaryEntrance =
181 primaryVertex->GetPosition() - primaryVertex->GetPosition().z() * primaryDirection;
182
183
192
193
195 G4double hitEn = 0;
196 G4double totalEnergy = 0;
197 G4int hitNum = 0;
198 G4int totalNum = 0;
199 G4int hitZ = -1;
200 G4int hitRho = -1;
201 G4int hitPhi = -1;
202 G4int hitType = -1;
203 G4int numNonZeroThresholdCells = 0;
204 G4double tDistance = 0., rDistance = 0., phiDistance = 0.;
205 G4double tFirstMoment = 0., tSecondMoment = 0.;
206 G4double rFirstMoment = 0., rSecondMoment = 0.;
207 G4double phiMean = 0.;
208 for(size_t iHit = 0; iHit < hitsCollection->entries(); iHit++)
209 {
210 hit =
static_cast<Par04Hit*
>(hitsCollection->GetHit(iHit));
217 if(hitEn > 0)
218 {
219 totalEnergy += hitEn;
220 totalNum += hitNum;
224 tFirstMoment += hitEn * tDistance;
225 rFirstMoment += hitEn * rDistance;
226 phiMean += hitEn * phiDistance;
227 analysisManager->FillH1(4, tDistance, hitEn);
228 analysisManager->FillH1(5, rDistance, hitEn);
229 analysisManager->FillH1(10, hitType);
230 if(hitEn > 0.0005)
231 {
232 fCalEdep[numNonZeroThresholdCells] = hitEn;
233 fCalRho[numNonZeroThresholdCells] = hitRho;
234 fCalPhi[numNonZeroThresholdCells] = hitPhi;
235 fCalZ[numNonZeroThresholdCells] = hitZ;
236 numNonZeroThresholdCells++;
237 analysisManager->FillH1(13, std::log10(hitEn));
238 analysisManager->FillH1(15, hitNum);
239 }
240 }
241 }
242 tFirstMoment /= totalEnergy;
243 rFirstMoment /= totalEnergy;
244 phiMean /= totalEnergy;
245 analysisManager->FillH1(0, primaryEnergy / GeV);
246 analysisManager->FillH1(1, totalEnergy / GeV);
247 analysisManager->FillH1(2, totalEnergy / primaryEnergy);
248 analysisManager->FillH1(3,
fTimer.GetRealElapsed());
249 analysisManager->FillH1(6, tFirstMoment);
250 analysisManager->FillH1(7, rFirstMoment);
251 analysisManager->FillH1(12, numNonZeroThresholdCells);
252 analysisManager->FillH1(14, totalNum);
253
254 fCalEdep.resize(numNonZeroThresholdCells);
255 fCalRho.resize(numNonZeroThresholdCells);
256 fCalPhi.resize(numNonZeroThresholdCells);
257 fCalZ.resize(numNonZeroThresholdCells);
258 analysisManager->FillNtupleDColumn(0, 0, primaryEnergy);
259 analysisManager->FillNtupleDColumn(0, 1,
fTimer.GetRealElapsed());
260
261 for(size_t iHit = 0; iHit < hitsCollection->entries(); iHit++)
262 {
263 hit =
static_cast<Par04Hit*
>(hitsCollection->GetHit(iHit));
268 if(hitEn > 0)
269 {
273 tSecondMoment += hitEn * std::pow(tDistance - tFirstMoment, 2);
274 rSecondMoment += hitEn * std::pow(rDistance - rFirstMoment, 2);
275 analysisManager->FillH1(11, phiDistance - phiMean, hitEn);
276 }
277 }
278 tSecondMoment /= totalEnergy;
279 rSecondMoment /= totalEnergy;
280 analysisManager->FillH1(8, tSecondMoment);
281 analysisManager->FillH1(9, rSecondMoment);
282
283
284 G4double totalPhysicalEnergy = 0;
285 totalNum = 0;
286 hitEn = 0;
287 hitNum = 0;
288 G4int hitLayer = -1;
289 G4int hitRow = -1;
290 G4int hitSlice = -1;
291 numNonZeroThresholdCells = 0;
292 for(size_t iHit = 0; iHit < physicalFullHitsCollection->entries(); iHit++)
293 {
294 hit =
static_cast<Par04Hit*
>(physicalFullHitsCollection->GetHit(iHit));
300 if(hitEn > 0)
301 {
302 totalPhysicalEnergy += hitEn;
303 totalNum += hitNum;
304 if(hitEn > 0.0005)
305 {
310 numNonZeroThresholdCells++;
311 analysisManager->FillH1(19, std::log10(hitEn));
312 analysisManager->FillH1(21, hitNum);
313 }
314 }
315 }for(size_t iHit = 0; iHit < physicalFastHitsCollection->entries(); iHit++)
316 {
317 hit =
static_cast<Par04Hit*
>(physicalFastHitsCollection->GetHit(iHit));
323 if(hitEn > 0)
324 {
325 totalPhysicalEnergy += hitEn;
326 totalNum += hitNum;
327 if(hitEn > 0.0005)
328 {
333 numNonZeroThresholdCells++;
334 analysisManager->FillH1(19, std::log10(hitEn));
335 analysisManager->FillH1(21, hitNum);
336 }
337 }
338 }
339 analysisManager->FillH1(16, totalPhysicalEnergy / GeV);
340 analysisManager->FillH1(17, totalPhysicalEnergy / primaryEnergy);
341 analysisManager->FillH1(18, numNonZeroThresholdCells);
342 analysisManager->FillH1(20, totalNum);
347 analysisManager->AddNtupleRow(0);
348 analysisManager->AddNtupleRow(1);
349 analysisManager->AddNtupleRow(2);
350}
G4THitsCollection< Par04Hit > Par04HitsCollection
G4ThreeVector GetMeshSizeOfCells() const
Get size of Mesh cells in cylindrical coordinates (r, phi, z)
G4double fCellSizePhi
Size of cell in azimuthal angle.
G4double fCellSizeRho
Size of cell along radius of cylinder.
G4double fCellSizeZ
Size of cell along Z axis.
Hit class to store energy deposited in the sensitive detector.
G4int GetZid() const
Get Z id of the cell in the readout segmentation.
G4double GetEdep() const
Get energy.
G4int GetType() const
Get type (0 = full sim, 1 = fast sim)
G4int GetRhoId() const
Get rho id of the cell in the readout segmentation.
G4int GetNdep() const
Get number of deposits per hit/cell.
G4int GetPhiId() const
Get phi id of the cell in the readout segmentation.
G4int GetNbOfRows() const
Get number of rows.
G4int GetNbOfLayers() const
Get number of layers.
G4int GetNbOfSlices() const
Get number of slices.