| 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/Hadr02/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager
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.Ivanchen
39 // 16.11.2006 Add beamFlag (V.Ivanchenko)
40 // 16.10.2012 Renamed G4IonFTFPBinaryCascadePh
41 //
42 //--------------------------------------------
43 //
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
85
86 HistoManager* HistoManager::fManager = 0;
87
88 //....oooOO0OOooo........oooOO0OOooo........oo
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........oo
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........oo
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........oo
150
151 HistoManager::~HistoManager()
152 {
153 delete fHisto;
154 }
155
156 //....oooOO0OOooo........oooOO0OOooo........oo
157
158 void HistoManager::BookHisto()
159 {
160 fHisto->Add1D("0", "Energy deposition (MeV/m
161 MeV / mm);
162 fHisto->Add1D("1", "Log10 Energy (GeV) of ga
163 fHisto->Add1D("2", "Log10 Energy (GeV) of el
164 fHisto->Add1D("3", "Log10 Energy (GeV) of po
165 fHisto->Add1D("4", "Log10 Energy (GeV) of pr
166 fHisto->Add1D("5", "Log10 Energy (GeV) of ne
167 fHisto->Add1D("6", "Log10 Energy (GeV) of ch
168 fHisto->Add1D("7", "Log10 Energy (GeV) of pi
169 fHisto->Add1D("8", "Log10 Energy (GeV) of ch
170 fHisto->Add1D("9", "Log10 Energy (GeV) of ne
171 fHisto->Add1D("10", "Log10 Energy (GeV) of d
172 fHisto->Add1D("11", "Log10 Energy (GeV) of H
173 fHisto->Add1D("12", "Log10 Energy (GeV) of G
174 fHisto->Add1D("13", "Log10 Energy (GeV) of m
175 fHisto->Add1D("14", "Log10 Energy (GeV) of p
176 fHisto->Add1D("15", "Log10 Energy (GeV) of p
177 fHisto->Add1D("16", "X Section (mb) of Secon
178 1.0);
179 fHisto->Add1D("17", "Secondary Fragment A E>
180 fHisto->Add1D("18", "Secondary Fragment Z E<
181 fHisto->Add1D("19", "Secondary Fragment A E<
182 fHisto->Add1D("20", "X Section (mb) of Secon
183 fHisto->Add1D("21", "Secondary Fragment A ",
184 }
185
186 //....oooOO0OOooo........oooOO0OOooo........oo
187
188 void HistoManager::BeginOfRun()
189 {
190 Initialise();
191 BookHisto();
192 fHisto->Book();
193
194 if (fVerbose > 0) {
195 G4cout << "HistoManager: Histograms are bo
196 << " BeginOfRun (After fHisto->boo
197 }
198 }
199
200 //....oooOO0OOooo........oooOO0OOooo........oo
201
202 void HistoManager::EndOfRun()
203 {
204 G4cout << "HistoManager: End of run actions
205
206 // Average values
207 G4cout << "=================================
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
238 G4cout << "Beam Energy(GeV)
239 G4cout << "Number of events
240 G4cout << std::setprecision(4) << "Average e
241 << " RMS(GeV) " << fEdepSum2 / GeV
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 << "=================================
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 elem
264 // G4cout << fMaterial << G4endl;
265 G4double F = 1000 / (fLength * barn * fMater
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........oo
275
276 void HistoManager::BeginOfEvent()
277 {
278 ++fNevt;
279 fEdepEvt = 0.0;
280 }
281
282 //....oooOO0OOooo........oooOO0OOooo........oo
283
284 void HistoManager::EndOfEvent()
285 {
286 fEdepSum += fEdepEvt;
287 fEdepSum2 += fEdepEvt * fEdepEvt;
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oo
291
292 void HistoManager::ScoreNewTrack(const G4Track
293 {
294 const G4ParticleDefinition* pd = track->GetD
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->GetMomentumDire
303 if (1 < fVerbose)
304 G4cout << "### Primary " << name << " ki
305 << "; m(GeV)= " << pd->GetPDGMass
306 << "; dir= " << track->GetMoment
307
308 // Secondary track
309 }
310 else {
311 if (1 < fVerbose) {
312 G4cout << "=== Secondary " << name << "
313 << "; m(GeV)= " << pd->GetPDGMass
314 << "; dir= " << track->GetMoment
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() || p
355 fNkaons++;
356 fHisto->Fill(8, e, 1.0);
357 }
358 else if (pd == G4KaonZeroShort::KaonZeroSh
359 fNkaons++;
360 fHisto->Fill(9, e, 1.0);
361 }
362 else if (pd == G4Deuteron::Deuteron() || p
363 fNdeut++;
364 fHisto->Fill(10, e, 1.0);
365 }
366 else if (pd == G4He3::He3() || pd == G4Alp
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->GetBaryonNumb
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() || p
387 fNmuons++;
388 fHisto->Fill(13, e, 1.0);
389 }
390 }
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oo
394
395 void HistoManager::AddTargetStep(const G4Step*
396 {
397 ++fNstep;
398 G4double edep = step->GetTotalEnergyDeposit(
399
400 if (edep > 0.0) {
401 G4ThreeVector pos =
402 (step->GetPreStepPoint()->GetPosition()
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(ke
413 << "; step(mm)= " << step->GetSte
414 << track->GetDefinition()->GetPar
415 << " E(GeV)= " << track->GetKinet
416 }
417 }
418 }
419
420 //....oooOO0OOooo........oooOO0OOooo........oo
421
422 void HistoManager::SetVerbose(G4int val)
423 {
424 fVerbose = val;
425 fHisto->SetVerbose(val);
426 }
427
428 //....oooOO0OOooo........oooOO0OOooo........oo
429
430 void HistoManager::Fill(G4int id, G4double x,
431 {
432 fHisto->Fill(id, x, w);
433 }
434
435 //....oooOO0OOooo........oooOO0OOooo........oo
436
437 void HistoManager::SetIonPhysics(const G4Strin
438 {
439 if (fIonPhysics) {
440 G4cout << "### HistoManager WARNING: Ion P
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()->PhysicsHasB
448 G4cout << "### SetIonPhysics: Ion Physics
449 #else
450 G4cout << "Error: Ion Physics HIJING is re
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()->PhysicsHasB
458 G4cout << "### SetIonPhysics: Ion Physics
459 #else
460 G4cout << "Error: Ion Physics UrQMD is req
461 #endif
462 }
463 else {
464 G4cout << "### HistoManager WARNING: Ion P
465 }
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oo
469