| 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 /// \file hadronic/Hadr01/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager class
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.Ivanchenko)
38 // 16.11.2006 Add beamFlag (V.Ivanchenko)
39 //
40 //----------------------------------------------------------------------------
41 //
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
77
78 HistoManager* HistoManager::fManager = nullptr;
79
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
105
106 HistoManager::~HistoManager()
107 {
108 delete fHisto;
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112
113 void HistoManager::BookHisto()
114 {
115 fHistoBooked = true;
116 fHisto->Add1D("1", "Energy deposition (MeV/mm/event) in the target", fNSlices, 0.0, fLength / mm,
117 MeV / mm);
118 fHisto->Add1D("2", "Log10 Energy (MeV) of gammas", fNBinsE, -5., 5., 1.0);
119 fHisto->Add1D("3", "Log10 Energy (MeV) of electrons", fNBinsE, -5., 5., 1.0);
120 fHisto->Add1D("4", "Log10 Energy (MeV) of positrons", fNBinsE, -5., 5., 1.0);
121 fHisto->Add1D("5", "Log10 Energy (MeV) of protons", fNBinsE, -5., 5., 1.0);
122 fHisto->Add1D("6", "Log10 Energy (MeV) of neutrons", fNBinsE, -5., 5., 1.0);
123 fHisto->Add1D("7", "Log10 Energy (MeV) of charged pions", fNBinsE, -4., 6., 1.0);
124 fHisto->Add1D("8", "Log10 Energy (MeV) of pi0", fNBinsE, -4., 6., 1.0);
125 fHisto->Add1D("9", "Log10 Energy (MeV) of charged kaons", fNBinsE, -4., 6., 1.0);
126 fHisto->Add1D("10", "Log10 Energy (MeV) of neutral kaons", fNBinsE, -4., 6., 1.0);
127 fHisto->Add1D("11", "Log10 Energy (MeV) of deuterons and tritons", fNBinsE, -5., 5., 1.0);
128 fHisto->Add1D("12", "Log10 Energy (MeV) of He3 and alpha", fNBinsE, -5., 5., 1.0);
129 fHisto->Add1D("13", "Log10 Energy (MeV) of Generic Ions", fNBinsE, -5., 5., 1.0);
130 fHisto->Add1D("14", "Log10 Energy (MeV) of muons", fNBinsE, -4., 6., 1.0);
131 fHisto->Add1D("15", "log10 Energy (MeV) of side-leaked neutrons", fNBinsE, -5., 5., 1.0);
132 fHisto->Add1D("16", "log10 Energy (MeV) of forward-leaked neutrons", fNBinsE, -5., 5., 1.0);
133 fHisto->Add1D("17", "log10 Energy (MeV) of backward-leaked neutrons", fNBinsE, -5., 5., 1.0);
134 fHisto->Add1D("18", "log10 Energy (MeV) of leaking protons", fNBinsE, -4., 6., 1.0);
135 fHisto->Add1D("19", "log10 Energy (MeV) of leaking charged pions", fNBinsE, -4., 6., 1.0);
136 fHisto->Add1D("20", "Log10 Energy (MeV) of pi+", fNBinsE, -4., 6., 1.0);
137 fHisto->Add1D("21", "Log10 Energy (MeV) of pi-", fNBinsE, -4., 6., 1.0);
138 fHisto->Add1D("22", "Energy deposition in the target normalized to beam energy", 110, 0.0, 1.1,
139 1.0);
140 fHisto->Add1D("23", "EM energy deposition in the target normalized to beam energy", 110, 0.0, 1.1,
141 1.0);
142 fHisto->Add1D("24", "Pion energy deposition in the target normalized to beam energy", 110, 0.0,
143 1.1, 1.0);
144 fHisto->Add1D("25", "Proton energy deposition in the target normalized to beam energy", 110, 0.0,
145 1.1, 1.0);
146 fHisto->Add1D("26", "Energy (MeV) of pi+", fNBinsE, 0., 500., 1.0);
147 fHisto->Add1D("27", "Energy (MeV) of pi-", fNBinsE, 0., 500., 1.0);
148 fHisto->Add1D("28", "Energy (MeV) of pi0", fNBinsE, 0., 500., 1.0);
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
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 booked and run has been started" << G4endl;
188 }
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192
193 void HistoManager::EndOfRun()
194 {
195 G4cout << "HistoManager: End of run actions are started" << G4endl;
196
197 // Average values
198 G4cout << "========================================================" << G4endl;
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 " << fPrimaryDef->GetParticleName() << G4endl;
236 G4cout << "Beam Energy(MeV) " << fPrimaryKineticEnergy / MeV << G4endl;
237 G4cout << "Number of events " << fNevt << G4endl;
238 G4cout << std::setprecision(4) << "Average energy deposit (MeV) " << fEdepSum / MeV
239 << " RMS(MeV) " << fEdepSum2 / MeV << G4endl;
240 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
241 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
242 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
243 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
244 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
245 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
246 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
247 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
248 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
249 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
250 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
251 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
252 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
253 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
254 G4cout << std::setprecision(4) << "Average number of forward neutrons " << xnf << G4endl;
255 G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl;
256 G4cout << std::setprecision(4) << "Average number of leaked neutrons " << xnbw << G4endl;
257 G4cout << std::setprecision(4) << "Average number of proton leak " << xpl << G4endl;
258 G4cout << std::setprecision(4) << "Average number of pion leak " << xal << G4endl;
259 G4cout << "========================================================" << G4endl;
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........oooOO0OOooo........oooOO0OOooo....
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........oooOO0OOooo........oooOO0OOooo....
281
282 void HistoManager::EndOfEvent()
283 {
284 fEdepSum += fEdepEvt;
285 fEdepSum2 += fEdepEvt * fEdepEvt;
286 fHisto->Fill(21, fEdepEvt / fPrimaryKineticEnergy, 1.0);
287 fHisto->Fill(22, fEdepEM / fPrimaryKineticEnergy, 1.0);
288 fHisto->Fill(23, fEdepPI / fPrimaryKineticEnergy, 1.0);
289 fHisto->Fill(24, fEdepP / fPrimaryKineticEnergy, 1.0);
290 }
291
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
293
294 void HistoManager::ScoreNewTrack(const G4Track* track)
295 {
296 const G4ParticleDefinition* pd = track->GetDefinition();
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->GetMomentumDirection();
307 if (1 < fVerbose)
308 G4cout << "### Primary " << name << " kinE(MeV)= " << e / MeV
309 << "; m(MeV)= " << pd->GetPDGMass() / MeV << "; pos(mm)= " << track->GetPosition() / mm
310 << "; dir= " << track->GetMomentumDirection() << G4endl;
311
312 // Secondary track
313 }
314 else {
315 if (1 < fVerbose)
316 G4cout << "=== Secondary " << name << " kinE(MeV)= " << e / MeV
317 << "; m(MeV)= " << pd->GetPDGMass() / MeV << "; pos(mm)= " << track->GetPosition() / mm
318 << "; dir= " << track->GetMomentumDirection() << G4endl;
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() || pd == G4KaonMinus::KaonMinus()) {
361 fNkaons++;
362 fHisto->Fill(8, e, 1.0);
363 }
364 else if (pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
365 fNkaons++;
366 fHisto->Fill(9, e, 1.0);
367 }
368 else if (pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
369 fNdeut++;
370 fHisto->Fill(10, e, 1.0);
371 }
372 else if (pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
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() || pd == G4MuonMinus::MuonMinus()) {
381 fNmuons++;
382 fHisto->Fill(13, e, 1.0);
383 }
384 }
385 }
386
387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
388
389 void HistoManager::AddTargetStep(const G4Step* step)
390 {
391 ++fNstep;
392 const G4double fEdep = step->GetTotalEnergyDeposit();
393 const G4Track* track = step->GetTrack();
394 const G4ParticleDefinition* pd = track->GetDefinition();
395 if (1 < fVerbose) {
396 G4cout << "TargetSD::ProcessHits: " << pd->GetParticleName()
397 << " E(MeV)=" << track->GetKineticEnergy() / CLHEP::MeV
398 << " beta1= " << step->GetPreStepPoint()->GetVelocity() / CLHEP::c_light
399 << " beta2= " << step->GetPostStepPoint()->GetVelocity() / CLHEP::c_light
400 << " weight= " << step->GetTrack()->GetWeight()
401 << " t(ns)=" << track->GetGlobalTime() / CLHEP::ns << G4endl;
402 }
403 if (fEdep > 0.0) {
404 G4ThreeVector pos =
405 (step->GetPreStepPoint()->GetPosition() + step->GetPostStepPoint()->GetPosition()) * 0.5;
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 == G4Electron::Electron() || pd == G4Positron::Positron()) {
414 fEdepEM += fEdep;
415 }
416 else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
417 fEdepPI += fEdep;
418 }
419 else if (pd == G4Proton::Proton() || pd == G4AntiProton::AntiProton()) {
420 fEdepP += fEdep;
421 }
422
423 if (1 < fVerbose) {
424 G4cout << "HistoManager::AddEnergy: e(keV)= " << fEdep / keV << "; z(mm)= " << z / mm
425 << "; step(mm)= " << step->GetStepLength() / mm << " by " << pd->GetParticleName()
426 << " E(MeV)= " << track->GetKineticEnergy() / MeV << G4endl;
427 }
428 }
429 }
430
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432
433 void HistoManager::AddLeakingParticle(const G4Track* track)
434 {
435 const G4ParticleDefinition* pd = track->GetDefinition();
436 const G4StepPoint* sp = track->GetStep()->GetPreStepPoint();
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::AddLeakingParticle " << pd->GetParticleName()
444 << " Eprestep(MeV)=" << ekin << " trackID=" << track->GetTrackID() << G4endl;
445 return;
446 }
447
448 const G4ThreeVector& pos = sp->GetPosition();
449 const G4ThreeVector& dir = sp->GetMomentumDirection();
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 < 0.0) {
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 * vx + y * vy > 0.0
484 && std::abs(x * x + y * y - fR2) < del * fRadius)
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() || pd == G4PionMinus::PionMinus()) {
502 fHisto->Fill(18, e, 1.0);
503 ++fNpiofNleak;
504 }
505 }
506 }
507
508 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
509
510 void HistoManager::SetVerbose(G4int val)
511 {
512 fVerbose = val;
513 fHisto->SetVerbose(val);
514 }
515
516 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
517
518 void HistoManager::Fill(G4int id, G4double x, G4double w)
519 {
520 fHisto->Fill(id, x, w);
521 }
522
523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
524