| 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 /// \file hadronic/Hadr01/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager
28 //
29 //--------------------------------------------
30 //
31 // ClassName: HistoManager
32 //
33 //
34 // Author: V.Ivanchenko 30/01/01
35 //
36 // Modified:
37 // 04.06.2006 Adoptation of Hadr01 (V.Ivanchen
38 // 16.11.2006 Add beamFlag (V.Ivanchenko)
39 //
40 //--------------------------------------------
41 //
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 #include "HistoManager.hh"
47
48 #include "Histo.hh"
49
50 #include "G4Alpha.hh"
51 #include "G4AntiProton.hh"
52 #include "G4Deuteron.hh"
53 #include "G4Electron.hh"
54 #include "G4Gamma.hh"
55 #include "G4He3.hh"
56 #include "G4KaonMinus.hh"
57 #include "G4KaonPlus.hh"
58 #include "G4KaonZeroLong.hh"
59 #include "G4KaonZeroShort.hh"
60 #include "G4MuonMinus.hh"
61 #include "G4MuonPlus.hh"
62 #include "G4Neutron.hh"
63 #include "G4PhysicalConstants.hh"
64 #include "G4PionMinus.hh"
65 #include "G4PionPlus.hh"
66 #include "G4PionZero.hh"
67 #include "G4Positron.hh"
68 #include "G4Proton.hh"
69 #include "G4Step.hh"
70 #include "G4SystemOfUnits.hh"
71 #include "G4Track.hh"
72 #include "G4Triton.hh"
73 #include "G4UnitsTable.hh"
74 #include "globals.hh"
75
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 HistoManager* HistoManager::fManager = nullptr
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81
82 HistoManager* HistoManager::GetPointer()
83 {
84 if (nullptr == fManager) {
85 static HistoManager manager;
86 fManager = &manager;
87 }
88 return fManager;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oo
92
93 HistoManager::HistoManager()
94 : fPrimaryDef(nullptr),
95 fRadius(10 * CLHEP::cm),
96 fLength(300. * CLHEP::mm),
97 fEdepMax(1.0 * CLHEP::GeV)
98 {
99 fHisto = new Histo();
100 fHisto->SetVerbose(fVerbose);
101 fNeutron = G4Neutron::Neutron();
102 }
103
104 //....oooOO0OOooo........oooOO0OOooo........oo
105
106 HistoManager::~HistoManager()
107 {
108 delete fHisto;
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oo
112
113 void HistoManager::BookHisto()
114 {
115 fHistoBooked = true;
116 fHisto->Add1D("1", "Energy deposition (MeV/m
117 MeV / mm);
118 fHisto->Add1D("2", "Log10 Energy (MeV) of ga
119 fHisto->Add1D("3", "Log10 Energy (MeV) of el
120 fHisto->Add1D("4", "Log10 Energy (MeV) of po
121 fHisto->Add1D("5", "Log10 Energy (MeV) of pr
122 fHisto->Add1D("6", "Log10 Energy (MeV) of ne
123 fHisto->Add1D("7", "Log10 Energy (MeV) of ch
124 fHisto->Add1D("8", "Log10 Energy (MeV) of pi
125 fHisto->Add1D("9", "Log10 Energy (MeV) of ch
126 fHisto->Add1D("10", "Log10 Energy (MeV) of n
127 fHisto->Add1D("11", "Log10 Energy (MeV) of d
128 fHisto->Add1D("12", "Log10 Energy (MeV) of H
129 fHisto->Add1D("13", "Log10 Energy (MeV) of G
130 fHisto->Add1D("14", "Log10 Energy (MeV) of m
131 fHisto->Add1D("15", "log10 Energy (MeV) of s
132 fHisto->Add1D("16", "log10 Energy (MeV) of f
133 fHisto->Add1D("17", "log10 Energy (MeV) of b
134 fHisto->Add1D("18", "log10 Energy (MeV) of l
135 fHisto->Add1D("19", "log10 Energy (MeV) of l
136 fHisto->Add1D("20", "Log10 Energy (MeV) of p
137 fHisto->Add1D("21", "Log10 Energy (MeV) of p
138 fHisto->Add1D("22", "Energy deposition in th
139 1.0);
140 fHisto->Add1D("23", "EM energy deposition in
141 1.0);
142 fHisto->Add1D("24", "Pion energy deposition
143 1.1, 1.0);
144 fHisto->Add1D("25", "Proton energy depositio
145 1.1, 1.0);
146 fHisto->Add1D("26", "Energy (MeV) of pi+", f
147 fHisto->Add1D("27", "Energy (MeV) of pi-", f
148 fHisto->Add1D("28", "Energy (MeV) of pi0", f
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oo
152
153 void HistoManager::BeginOfRun()
154 {
155 fR2 = fRadius * fRadius;
156 fAbsZ0 = 0.5 * fLength;
157 fNevt = 0;
158 fNelec = 0;
159 fNposit = 0;
160 fNgam = 0;
161 fNstep = 0;
162 fNprot_leak = 0;
163 fNpiofNleak = 0;
164 fNions = 0;
165 fNdeut = 0;
166 fNalpha = 0;
167 fNkaons = 0;
168 fNmuons = 0;
169 fNcpions = 0;
170 fNpi0 = 0;
171 fNneutron = 0;
172 fNproton = 0;
173 fNaproton = 0;
174 fNneu_forw = 0;
175 fNneu_leak = 0;
176 fNneu_back = 0;
177
178 fEdepSum = 0.0;
179 fEdepSum2 = 0.0;
180
181 if (!fHistoBooked) {
182 BookHisto();
183 }
184 fHisto->Book();
185
186 if (fVerbose > 0) {
187 G4cout << "HistoManager: Histograms are bo
188 }
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oo
192
193 void HistoManager::EndOfRun()
194 {
195 G4cout << "HistoManager: End of run actions
196
197 // Average values
198 G4cout << "=================================
199
200 G4double x = (G4double)fNevt;
201 if (fNevt > 0) {
202 x = 1.0 / x;
203 }
204
205 G4double xe = x * (G4double)fNelec;
206 G4double xg = x * (G4double)fNgam;
207 G4double xp = x * (G4double)fNposit;
208 G4double xs = x * (G4double)fNstep;
209 G4double xn = x * (G4double)fNneutron;
210 G4double xpn = x * (G4double)fNproton;
211 G4double xap = x * (G4double)fNaproton;
212 G4double xnf = x * (G4double)fNneu_forw;
213 G4double xnb = x * (G4double)fNneu_leak;
214 G4double xnbw = x * (G4double)fNneu_back;
215 G4double xpl = x * (G4double)fNprot_leak;
216 G4double xal = x * (G4double)fNpiofNleak;
217 G4double xpc = x * (G4double)fNcpions;
218 G4double xp0 = x * (G4double)fNpi0;
219 G4double xpk = x * (G4double)fNkaons;
220 G4double xpm = x * (G4double)fNmuons;
221 G4double xid = x * (G4double)fNdeut;
222 G4double xia = x * (G4double)fNalpha;
223 G4double xio = x * (G4double)fNions;
224
225 fEdepSum *= x;
226 fEdepSum2 *= x;
227 fEdepSum2 -= fEdepSum * fEdepSum;
228 if (fEdepSum2 > 0.0) {
229 fEdepSum2 = std::sqrt(fEdepSum2);
230 }
231 else {
232 fEdepSum2 = 0.0;
233 }
234
235 G4cout << "Beam particle
236 G4cout << "Beam Energy(MeV)
237 G4cout << "Number of events
238 G4cout << std::setprecision(4) << "Average e
239 << " RMS(MeV) " << fEdepSum2 / MeV
240 G4cout << std::setprecision(4) << "Average n
241 G4cout << std::setprecision(4) << "Average n
242 G4cout << std::setprecision(4) << "Average n
243 G4cout << std::setprecision(4) << "Average n
244 G4cout << std::setprecision(4) << "Average n
245 G4cout << std::setprecision(4) << "Average n
246 G4cout << std::setprecision(4) << "Average n
247 G4cout << std::setprecision(4) << "Average n
248 G4cout << std::setprecision(4) << "Average n
249 G4cout << std::setprecision(4) << "Average n
250 G4cout << std::setprecision(4) << "Average n
251 G4cout << std::setprecision(4) << "Average n
252 G4cout << std::setprecision(4) << "Average n
253 G4cout << std::setprecision(4) << "Average n
254 G4cout << std::setprecision(4) << "Average n
255 G4cout << std::setprecision(4) << "Average n
256 G4cout << std::setprecision(4) << "Average n
257 G4cout << std::setprecision(4) << "Average n
258 G4cout << std::setprecision(4) << "Average n
259 G4cout << "=================================
260 G4cout << G4endl;
261
262 // normalise histograms
263 for (G4int i = 0; i < fNHisto; i++) {
264 fHisto->ScaleH1(i, x);
265 }
266
267 fHisto->Save();
268 }
269
270 //....oooOO0OOooo........oooOO0OOooo........oo
271
272 void HistoManager::BeginOfEvent()
273 {
274 fEdepEvt = 0.0;
275 fEdepEM = 0.0;
276 fEdepPI = 0.0;
277 fEdepP = 0.0;
278 }
279
280 //....oooOO0OOooo........oooOO0OOooo........oo
281
282 void HistoManager::EndOfEvent()
283 {
284 fEdepSum += fEdepEvt;
285 fEdepSum2 += fEdepEvt * fEdepEvt;
286 fHisto->Fill(21, fEdepEvt / fPrimaryKineticE
287 fHisto->Fill(22, fEdepEM / fPrimaryKineticEn
288 fHisto->Fill(23, fEdepPI / fPrimaryKineticEn
289 fHisto->Fill(24, fEdepP / fPrimaryKineticEne
290 }
291
292 //....oooOO0OOooo........oooOO0OOooo........oo
293
294 void HistoManager::ScoreNewTrack(const G4Track
295 {
296 const G4ParticleDefinition* pd = track->GetD
297 const G4String name = pd->GetParticleName();
298 const G4double ee = track->GetKineticEnergy(
299 G4double e = ee;
300
301 // Primary track
302 if (0 == track->GetParentID()) {
303 fNevt++;
304 fPrimaryKineticEnergy = e;
305 fPrimaryDef = pd;
306 G4ThreeVector dir = track->GetMomentumDire
307 if (1 < fVerbose)
308 G4cout << "### Primary " << name << " ki
309 << "; m(MeV)= " << pd->GetPDGMass
310 << "; dir= " << track->GetMoment
311
312 // Secondary track
313 }
314 else {
315 if (1 < fVerbose)
316 G4cout << "=== Secondary " << name << "
317 << "; m(MeV)= " << pd->GetPDGMass
318 << "; dir= " << track->GetMoment
319 e = (e > 0.0) ? std::log10(e / MeV) : -100
320 if (pd == G4Gamma::Gamma()) {
321 fNgam++;
322 fHisto->Fill(1, e, 1.0);
323 }
324 else if (pd == G4Electron::Electron()) {
325 fNelec++;
326 fHisto->Fill(2, e, 1.0);
327 }
328 else if (pd == G4Positron::Positron()) {
329 fNposit++;
330 fHisto->Fill(3, e, 1.0);
331 }
332 else if (pd == G4Proton::Proton()) {
333 fNproton++;
334 fHisto->Fill(4, e, 1.0);
335 }
336 else if (pd == fNeutron) {
337 fNneutron++;
338 fHisto->Fill(5, e, 1.0);
339 }
340 else if (pd == G4AntiProton::AntiProton())
341 fNaproton++;
342 }
343 else if (pd == G4PionPlus::PionPlus()) {
344 fNcpions++;
345 fHisto->Fill(6, e, 1.0);
346 fHisto->Fill(19, e, 1.0);
347 fHisto->Fill(25, ee, 1.0);
348 }
349 else if (pd == G4PionMinus::PionMinus()) {
350 fNcpions++;
351 fHisto->Fill(6, e, 1.0);
352 fHisto->Fill(20, e, 1.0);
353 fHisto->Fill(26, ee, 1.0);
354 }
355 else if (pd == G4PionZero::PionZero()) {
356 fNpi0++;
357 fHisto->Fill(7, e, 1.0);
358 fHisto->Fill(27, ee, 1.0);
359 }
360 else if (pd == G4KaonPlus::KaonPlus() || p
361 fNkaons++;
362 fHisto->Fill(8, e, 1.0);
363 }
364 else if (pd == G4KaonZeroShort::KaonZeroSh
365 fNkaons++;
366 fHisto->Fill(9, e, 1.0);
367 }
368 else if (pd == G4Deuteron::Deuteron() || p
369 fNdeut++;
370 fHisto->Fill(10, e, 1.0);
371 }
372 else if (pd == G4He3::He3() || pd == G4Alp
373 fNalpha++;
374 fHisto->Fill(11, e, 1.0);
375 }
376 else if (pd->GetParticleType() == "nucleus
377 fNions++;
378 fHisto->Fill(12, e, 1.0);
379 }
380 else if (pd == G4MuonPlus::MuonPlus() || p
381 fNmuons++;
382 fHisto->Fill(13, e, 1.0);
383 }
384 }
385 }
386
387 //....oooOO0OOooo........oooOO0OOooo........oo
388
389 void HistoManager::AddTargetStep(const G4Step*
390 {
391 ++fNstep;
392 const G4double fEdep = step->GetTotalEnergyD
393 const G4Track* track = step->GetTrack();
394 const G4ParticleDefinition* pd = track->GetD
395 if (1 < fVerbose) {
396 G4cout << "TargetSD::ProcessHits: " << pd-
397 << " E(MeV)=" << track->GetKineticE
398 << " beta1= " << step->GetPreStepPo
399 << " beta2= " << step->GetPostStepP
400 << " weight= " << step->GetTrack()-
401 << " t(ns)=" << track->GetGlobalTim
402 }
403 if (fEdep > 0.0) {
404 G4ThreeVector pos =
405 (step->GetPreStepPoint()->GetPosition()
406
407 G4double z = pos.z() + fAbsZ0;
408
409 // scoring
410 fEdepEvt += fEdep;
411 fHisto->Fill(0, z, fEdep);
412
413 if (pd == G4Gamma::Gamma() || pd == G4Elec
414 fEdepEM += fEdep;
415 }
416 else if (pd == G4PionPlus::PionPlus() || p
417 fEdepPI += fEdep;
418 }
419 else if (pd == G4Proton::Proton() || pd ==
420 fEdepP += fEdep;
421 }
422
423 if (1 < fVerbose) {
424 G4cout << "HistoManager::AddEnergy: e(ke
425 << "; step(mm)= " << step->GetSte
426 << " E(MeV)= " << track->GetKinet
427 }
428 }
429 }
430
431 //....oooOO0OOooo........oooOO0OOooo........oo
432
433 void HistoManager::AddLeakingParticle(const G4
434 {
435 const G4ParticleDefinition* pd = track->GetD
436 const G4StepPoint* sp = track->GetStep()->Ge
437 const G4double ekin = sp->GetKineticEnergy()
438 G4double e = -100.;
439 if (ekin > 0.0) {
440 e = std::log10(ekin / CLHEP::MeV);
441 }
442 else {
443 G4cout << "### HistoManager::AddLeakingPar
444 << " Eprestep(MeV)=" << ekin << "
445 return;
446 }
447
448 const G4ThreeVector& pos = sp->GetPosition()
449 const G4ThreeVector& dir = sp->GetMomentumDi
450 G4double x = pos.x();
451 G4double y = pos.y();
452 G4double z = pos.z();
453 G4double vx = dir.x();
454 G4double vy = dir.y();
455 G4double vz = dir.z();
456
457 G4bool isLeaking = false;
458
459 // Forward
460 const G4double del = 0.001 * CLHEP::mm;
461 if (std::abs(z - fAbsZ0) < del && vz > 0.0)
462 isLeaking = true;
463 if (pd == fNeutron) {
464 ++fNneu_forw;
465 fHisto->Fill(15, e, 1.0);
466 }
467 else
468 isLeaking = true;
469
470 // Backward
471 }
472 else if (std::abs(z + fAbsZ0) < del && vz <
473 isLeaking = true;
474 if (pd == fNeutron) {
475 ++fNneu_back;
476 fHisto->Fill(16, e, 1.0);
477 }
478 else
479 isLeaking = true;
480
481 // Side
482 }
483 else if (std::abs(z) <= fAbsZ0 + del && x *
484 && std::abs(x * x + y * y - fR2) <
485 {
486 isLeaking = true;
487 if (pd == fNeutron) {
488 ++fNneu_leak;
489 fHisto->Fill(14, e, 1.0);
490 }
491 else
492 isLeaking = true;
493 }
494
495 // protons and pions
496 if (isLeaking) {
497 if (pd == G4Proton::Proton()) {
498 fHisto->Fill(17, e, 1.0);
499 ++fNprot_leak;
500 }
501 else if (pd == G4PionPlus::PionPlus() || p
502 fHisto->Fill(18, e, 1.0);
503 ++fNpiofNleak;
504 }
505 }
506 }
507
508 //....oooOO0OOooo........oooOO0OOooo........oo
509
510 void HistoManager::SetVerbose(G4int val)
511 {
512 fVerbose = val;
513 fHisto->SetVerbose(val);
514 }
515
516 //....oooOO0OOooo........oooOO0OOooo........oo
517
518 void HistoManager::Fill(G4int id, G4double x,
519 {
520 fHisto->Fill(id, x, w);
521 }
522
523 //....oooOO0OOooo........oooOO0OOooo........oo
524