| 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/Hadr02/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager class
28 //
29 //
30 //---------------------------------------------------------------------------
31 //
32 // ClassName: HistoManager
33 //
34 //
35 // Author: V.Ivanchenko 30/01/01
36 //
37 // Modified:
38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchenko)
39 // 16.11.2006 Add beamFlag (V.Ivanchenko)
40 // 16.10.2012 Renamed G4IonFTFPBinaryCascadePhysics as G4IonPhysics (A.Ribon)
41 //
42 //----------------------------------------------------------------------------
43 //
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48 #include "HistoManager.hh"
49
50 #include "DetectorConstruction.hh"
51 #include "Histo.hh"
52 #include "IonHIJINGPhysics.hh"
53 #include "IonUrQMDPhysics.hh"
54
55 #include "G4Alpha.hh"
56 #include "G4AntiProton.hh"
57 #include "G4BuilderType.hh"
58 #include "G4Deuteron.hh"
59 #include "G4Electron.hh"
60 #include "G4Gamma.hh"
61 #include "G4He3.hh"
62 #include "G4KaonMinus.hh"
63 #include "G4KaonPlus.hh"
64 #include "G4KaonZeroLong.hh"
65 #include "G4KaonZeroShort.hh"
66 #include "G4Material.hh"
67 #include "G4MuonMinus.hh"
68 #include "G4MuonPlus.hh"
69 #include "G4Neutron.hh"
70 #include "G4NistManager.hh"
71 #include "G4PionMinus.hh"
72 #include "G4PionPlus.hh"
73 #include "G4PionZero.hh"
74 #include "G4Positron.hh"
75 #include "G4Proton.hh"
76 #include "G4RunManager.hh"
77 #include "G4SystemOfUnits.hh"
78 #include "G4Triton.hh"
79 #include "G4UnitsTable.hh"
80 #include "G4VModularPhysicsList.hh"
81 #include "G4VPhysicsConstructor.hh"
82 #include "globals.hh"
83
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85
86 HistoManager* HistoManager::fManager = 0;
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90 HistoManager* HistoManager::GetPointer()
91 {
92 if (!fManager) {
93 static HistoManager manager;
94 fManager = &manager;
95 }
96 return fManager;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100
101 HistoManager::HistoManager()
102 {
103 fVerbose = 0;
104 fNSlices = 100;
105 fNBinsE = 100;
106 fNHisto = 20;
107 fLength = 300. * mm;
108 fEdepMax = 1.0 * GeV;
109
110 fPrimaryDef = 0;
111 fPrimaryKineticEnergy = 0.0;
112 fMaterial = 0;
113 fBeamFlag = true;
114
115 fHisto = new Histo();
116 fHisto->SetVerbose(fVerbose);
117 fNeutron = G4Neutron::Neutron();
118 fPhysList = 0;
119 fIonPhysics = 0;
120 Initialise();
121 }
122
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124
125 void HistoManager::Initialise()
126 {
127 fAbsZ0 = -0.5 * fLength;
128 fNevt = 0;
129 fNelec = 0;
130 fNposit = 0;
131 fNgam = 0;
132 fNstep = 0;
133 fNions = 0;
134 fNdeut = 0;
135 fNalpha = 0;
136 fNkaons = 0;
137 fNmuons = 0;
138 fNcpions = 0;
139 fNpi0 = 0;
140 fNneutron = 0;
141 fNproton = 0;
142 fNaproton = 0;
143
144 fEdepEvt = 0.0;
145 fEdepSum = 0.0;
146 fEdepSum2 = 0.0;
147 }
148
149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150
151 HistoManager::~HistoManager()
152 {
153 delete fHisto;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157
158 void HistoManager::BookHisto()
159 {
160 fHisto->Add1D("0", "Energy deposition (MeV/mm/event) in the target", fNSlices, 0.0, fLength / mm,
161 MeV / mm);
162 fHisto->Add1D("1", "Log10 Energy (GeV) of gammas", fNBinsE, -5., 5., 1.0);
163 fHisto->Add1D("2", "Log10 Energy (GeV) of electrons", fNBinsE, -5., 5., 1.0);
164 fHisto->Add1D("3", "Log10 Energy (GeV) of positrons", fNBinsE, -5., 5., 1.0);
165 fHisto->Add1D("4", "Log10 Energy (GeV) of protons", fNBinsE, -5., 5., 1.0);
166 fHisto->Add1D("5", "Log10 Energy (GeV) of neutrons", fNBinsE, -5., 5., 1.0);
167 fHisto->Add1D("6", "Log10 Energy (GeV) of charged pions", fNBinsE, -4., 6., 1.0);
168 fHisto->Add1D("7", "Log10 Energy (GeV) of pi0", fNBinsE, -4., 6., 1.0);
169 fHisto->Add1D("8", "Log10 Energy (GeV) of charged kaons", fNBinsE, -4., 6., 1.0);
170 fHisto->Add1D("9", "Log10 Energy (GeV) of neutral kaons", fNBinsE, -4., 6., 1.0);
171 fHisto->Add1D("10", "Log10 Energy (GeV) of deuterons and tritons", fNBinsE, -5., 5., 1.0);
172 fHisto->Add1D("11", "Log10 Energy (GeV) of He3 and alpha", fNBinsE, -5., 5., 1.0);
173 fHisto->Add1D("12", "Log10 Energy (GeV) of Generic Ions", fNBinsE, -5., 5., 1.0);
174 fHisto->Add1D("13", "Log10 Energy (GeV) of muons", fNBinsE, -4., 6., 1.0);
175 fHisto->Add1D("14", "Log10 Energy (GeV) of pi+", fNBinsE, -4., 6., 1.0);
176 fHisto->Add1D("15", "Log10 Energy (GeV) of pi-", fNBinsE, -4., 6., 1.0);
177 fHisto->Add1D("16", "X Section (mb) of Secondary Fragments Z with E>1 GeV (mb)", 25, 0.5, 25.5,
178 1.0);
179 fHisto->Add1D("17", "Secondary Fragment A E>1 GeV", 50, 0.5, 50.5, 1.0);
180 fHisto->Add1D("18", "Secondary Fragment Z E<1 GeV", 25, 0.5, 25.5, 1.0);
181 fHisto->Add1D("19", "Secondary Fragment A E<1 GeV", 50, 0.5, 50.5, 1.0);
182 fHisto->Add1D("20", "X Section (mb) of Secondary Fragments Z (mb) ", 25, 0.5, 25.5, 1.0);
183 fHisto->Add1D("21", "Secondary Fragment A ", 50, 0.5, 50.5, 1.0);
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187
188 void HistoManager::BeginOfRun()
189 {
190 Initialise();
191 BookHisto();
192 fHisto->Book();
193
194 if (fVerbose > 0) {
195 G4cout << "HistoManager: Histograms are booked and run has been started" << G4endl
196 << " BeginOfRun (After fHisto->book)" << G4endl;
197 }
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201
202 void HistoManager::EndOfRun()
203 {
204 G4cout << "HistoManager: End of run actions are started" << G4endl;
205
206 // Average values
207 G4cout << "========================================================" << G4endl;
208
209 G4double x = (G4double)fNevt;
210 if (fNevt > 0) {
211 x = 1.0 / x;
212 }
213
214 G4double xe = x * fNelec;
215 G4double xg = x * fNgam;
216 G4double xp = x * fNposit;
217 G4double xs = x * fNstep;
218 G4double xn = x * fNneutron;
219 G4double xpn = x * fNproton;
220 G4double xap = x * fNaproton;
221 G4double xpc = x * fNcpions;
222 G4double xp0 = x * fNpi0;
223 G4double xpk = x * fNkaons;
224 G4double xpm = x * fNmuons;
225 G4double xid = x * fNdeut;
226 G4double xia = x * fNalpha;
227 G4double xio = x * fNions;
228
229 fEdepSum *= x;
230 fEdepSum2 *= x;
231 fEdepSum2 -= fEdepSum * fEdepSum;
232 if (fEdepSum2 > 0.0)
233 fEdepSum2 = std::sqrt(fEdepSum2);
234 else
235 fEdepSum2 = 0.0;
236
237 G4cout << "Beam particle " << fPrimaryDef->GetParticleName() << G4endl;
238 G4cout << "Beam Energy(GeV) " << fPrimaryKineticEnergy / GeV << G4endl;
239 G4cout << "Number of events " << fNevt << G4endl;
240 G4cout << std::setprecision(4) << "Average energy deposit (GeV) " << fEdepSum / GeV
241 << " RMS(GeV) " << fEdepSum2 / GeV << G4endl;
242 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
243 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
244 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
245 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
246 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
247 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
248 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
249 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
250 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
251 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
252 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
253 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
254 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
255 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
256 G4cout << "========================================================" << G4endl;
257 G4cout << G4endl;
258
259 // normalise histograms
260 for (G4int i = 0; i < fNHisto; i++) {
261 fHisto->ScaleH1(i, x);
262 }
263 // will work only for pure material - 1 element
264 // G4cout << fMaterial << G4endl;
265 G4double F = 1000 / (fLength * barn * fMaterial->GetTotNbOfAtomsPerVolume());
266 if (F > 0.0) {
267 fHisto->ScaleH1(16, F);
268 fHisto->ScaleH1(20, F);
269 }
270
271 fHisto->Save();
272 }
273
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
275
276 void HistoManager::BeginOfEvent()
277 {
278 ++fNevt;
279 fEdepEvt = 0.0;
280 }
281
282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283
284 void HistoManager::EndOfEvent()
285 {
286 fEdepSum += fEdepEvt;
287 fEdepSum2 += fEdepEvt * fEdepEvt;
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
291
292 void HistoManager::ScoreNewTrack(const G4Track* track)
293 {
294 const G4ParticleDefinition* pd = track->GetDefinition();
295 G4String name = pd->GetParticleName();
296 G4double e = track->GetKineticEnergy();
297
298 // Primary track
299 if (0 == track->GetParentID()) {
300 fPrimaryKineticEnergy = e;
301 fPrimaryDef = pd;
302 G4ThreeVector dir = track->GetMomentumDirection();
303 if (1 < fVerbose)
304 G4cout << "### Primary " << name << " kinE(GeV)= " << e / GeV
305 << "; m(GeV)= " << pd->GetPDGMass() / GeV << "; pos(mm)= " << track->GetPosition() / mm
306 << "; dir= " << track->GetMomentumDirection() << G4endl;
307
308 // Secondary track
309 }
310 else {
311 if (1 < fVerbose) {
312 G4cout << "=== Secondary " << name << " kinE(GeV)= " << e / GeV
313 << "; m(GeV)= " << pd->GetPDGMass() / GeV << "; pos(mm)= " << track->GetPosition() / mm
314 << "; dir= " << track->GetMomentumDirection() << G4endl;
315 }
316 e = std::log10(e / GeV);
317 if (pd == G4Gamma::Gamma()) {
318 fNgam++;
319 fHisto->Fill(1, e, 1.0);
320 }
321 else if (pd == G4Electron::Electron()) {
322 fNelec++;
323 fHisto->Fill(2, e, 1.0);
324 }
325 else if (pd == G4Positron::Positron()) {
326 fNposit++;
327 fHisto->Fill(3, e, 1.0);
328 }
329 else if (pd == G4Proton::Proton()) {
330 fNproton++;
331 fHisto->Fill(4, e, 1.0);
332 }
333 else if (pd == fNeutron) {
334 fNneutron++;
335 fHisto->Fill(5, e, 1.0);
336 }
337 else if (pd == G4AntiProton::AntiProton()) {
338 fNaproton++;
339 }
340 else if (pd == G4PionPlus::PionPlus()) {
341 fNcpions++;
342 fHisto->Fill(6, e, 1.0);
343 fHisto->Fill(14, e, 1.0);
344 }
345 else if (pd == G4PionMinus::PionMinus()) {
346 fNcpions++;
347 fHisto->Fill(6, e, 1.0);
348 fHisto->Fill(15, e, 1.0);
349 }
350 else if (pd == G4PionZero::PionZero()) {
351 fNpi0++;
352 fHisto->Fill(7, e, 1.0);
353 }
354 else if (pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
355 fNkaons++;
356 fHisto->Fill(8, e, 1.0);
357 }
358 else if (pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
359 fNkaons++;
360 fHisto->Fill(9, e, 1.0);
361 }
362 else if (pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
363 fNdeut++;
364 fHisto->Fill(10, e, 1.0);
365 }
366 else if (pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
367 fNalpha++;
368 fHisto->Fill(11, e, 1.0);
369 }
370 else if (pd->GetParticleType() == "nucleus") {
371 fNions++;
372 fHisto->Fill(12, e, 1.0);
373 G4double Z = pd->GetPDGCharge() / eplus;
374 G4double A = (G4double)pd->GetBaryonNumber();
375 if (e > 0.0) {
376 fHisto->Fill(16, Z, 1.0);
377 fHisto->Fill(17, A, 1.0);
378 }
379 else {
380 fHisto->Fill(18, Z, 1.0);
381 fHisto->Fill(19, A, 1.0);
382 }
383 fHisto->Fill(20, Z, 1.0);
384 fHisto->Fill(21, A, 1.0);
385 }
386 else if (pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
387 fNmuons++;
388 fHisto->Fill(13, e, 1.0);
389 }
390 }
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394
395 void HistoManager::AddTargetStep(const G4Step* step)
396 {
397 ++fNstep;
398 G4double edep = step->GetTotalEnergyDeposit();
399
400 if (edep > 0.0) {
401 G4ThreeVector pos =
402 (step->GetPreStepPoint()->GetPosition() + step->GetPostStepPoint()->GetPosition()) * 0.5;
403
404 G4double z = pos.z() - fAbsZ0;
405
406 // scoring
407 fEdepEvt += edep;
408 fHisto->Fill(0, z, edep);
409
410 if (1 < fVerbose) {
411 const G4Track* track = step->GetTrack();
412 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep / keV << "; z(mm)= " << z / mm
413 << "; step(mm)= " << step->GetStepLength() / mm << " by "
414 << track->GetDefinition()->GetParticleName()
415 << " E(GeV)= " << track->GetKineticEnergy() / GeV << G4endl;
416 }
417 }
418 }
419
420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
421
422 void HistoManager::SetVerbose(G4int val)
423 {
424 fVerbose = val;
425 fHisto->SetVerbose(val);
426 }
427
428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
429
430 void HistoManager::Fill(G4int id, G4double x, G4double w)
431 {
432 fHisto->Fill(id, x, w);
433 }
434
435 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436
437 void HistoManager::SetIonPhysics(const G4String& nam)
438 {
439 if (fIonPhysics) {
440 G4cout << "### HistoManager WARNING: Ion Physics is already defined: <" << nam
441 << "> is ignored!" << G4endl;
442 }
443 else if (nam == "HIJING") {
444 #ifdef G4_USE_HIJING
445 fIonPhysics = new IonHIJINGPhysics();
446 fPhysList->ReplacePhysics(fIonPhysics);
447 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
448 G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added" << G4endl;
449 #else
450 G4cout << "Error: Ion Physics HIJING is requested but is not available" << G4endl;
451 #endif
452 }
453 else if (nam == "UrQMD") {
454 #ifdef G4_USE_URQMD
455 fIonPhysics = new IonUrQMDPhysics(1);
456 fPhysList->ReplacePhysics(fIonPhysics);
457 G4RunManager::GetRunManager()->PhysicsHasBeenModified();
458 G4cout << "### SetIonPhysics: Ion Physics UrQMD is added" << G4endl;
459 #else
460 G4cout << "Error: Ion Physics UrQMD is requested but is not available" << G4endl;
461 #endif
462 }
463 else {
464 G4cout << "### HistoManager WARNING: Ion Physics <" << nam << "> is unknown!" << G4endl;
465 }
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469