source: trunk/examples/extended/hadronic/Hadr01/src/HistoManager.cc@ 1330

Last change on this file since 1330 was 1230, checked in by garnier, 16 years ago

update to geant4.9.3
File size: 16.4 KB
Line 
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// $Id: HistoManager.cc,v 1.16 2009/09/02 10:21:32 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-03-cand-01 $
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#include "G4UnitsTable.hh"
48#include "G4Neutron.hh"
49#include "G4Proton.hh"
50#include "G4AntiProton.hh"
51#include "G4Gamma.hh"
52#include "G4Electron.hh"
53#include "G4Positron.hh"
54#include "G4MuonPlus.hh"
55#include "G4MuonMinus.hh"
56#include "G4PionPlus.hh"
57#include "G4PionMinus.hh"
58#include "G4PionZero.hh"
59#include "G4KaonPlus.hh"
60#include "G4KaonMinus.hh"
61#include "G4KaonZeroShort.hh"
62#include "G4KaonZeroLong.hh"
63#include "G4Deuteron.hh"
64#include "G4Triton.hh"
65#include "G4He3.hh"
66#include "G4Alpha.hh"
67#include "Histo.hh"
68#include "globals.hh"
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72HistoManager* HistoManager::fManager = 0;
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76HistoManager* HistoManager::GetPointer()
77{
78 if(!fManager) {
79 static HistoManager manager;
80 fManager = &manager;
81 }
82 return fManager;
83}
84
85//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87HistoManager::HistoManager()
88{
89 verbose= 0;
90 nSlices = 300;
91 nBinsE = 100;
92 nHisto = 25;
93 length = 300.*mm;
94 edepMax = 1.0*GeV;
95 beamFlag = true;
96 material = 0;
97 elm = 0;
98 histo = new Histo(verbose);
99 neutron = G4Neutron::Neutron();
100}
101
102//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103
104HistoManager::~HistoManager()
105{
106 delete histo;
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
111void HistoManager::bookHisto()
112{
113 histo->add1D("1","Energy deposition (MeV/mm/event) in the target",
114 nSlices,0.0,length/mm,MeV/mm);
115 histo->add1D("2","Log10 Energy (MeV) of gammas",nBinsE,-5.,5.,1.0);
116 histo->add1D("3","Log10 Energy (MeV) of electrons",nBinsE,-5.,5.,1.0);
117 histo->add1D("4","Log10 Energy (MeV) of positrons",nBinsE,-5.,5.,1.0);
118 histo->add1D("5","Log10 Energy (MeV) of protons",nBinsE,-5.,5.,1.0);
119 histo->add1D("6","Log10 Energy (MeV) of neutrons",nBinsE,-5.,5.,1.0);
120 histo->add1D("7","Log10 Energy (MeV) of charged pions",nBinsE,-4.,6.,1.0);
121 histo->add1D("8","Log10 Energy (MeV) of pi0",nBinsE,-4.,6.,1.0);
122 histo->add1D("9","Log10 Energy (MeV) of charged kaons",nBinsE,-4.,6.,1.0);
123 histo->add1D("10","Log10 Energy (MeV) of neutral kaons",nBinsE,-4.,6.,1.0);
124 histo->add1D("11","Log10 Energy (MeV) of deuterons and tritons",nBinsE,-5.,5.,1.0);
125 histo->add1D("12","Log10 Energy (MeV) of He3 and alpha",nBinsE,-5.,5.,1.0);
126 histo->add1D("13","Log10 Energy (MeV) of Generic Ions",nBinsE,-5.,5.,1.0);
127 histo->add1D("14","Log10 Energy (MeV) of muons",nBinsE,-4.,6.,1.0);
128 histo->add1D("15","log10 Energy (MeV) of side-leaked neutrons",nBinsE,-5.,5.,1.0);
129 histo->add1D("16","log10 Energy (MeV) of forward-leaked neutrons",nBinsE,-5.,5.,1.0);
130 histo->add1D("17","log10 Energy (MeV) of backward-leaked neutrons",nBinsE,-5.,5.,1.0);
131 histo->add1D("18","log10 Energy (MeV) of leaking protons",nBinsE,-4.,6.,1.0);
132 histo->add1D("19","log10 Energy (MeV) of leaking charged pions",nBinsE,-4.,6.,1.0);
133 histo->add1D("20","Log10 Energy (MeV) of pi+",nBinsE,-4.,6.,1.0);
134 histo->add1D("21","Log10 Energy (MeV) of pi-",nBinsE,-4.,6.,1.0);
135 histo->add1D("22","Energy deposition in the target normalized to beam energy",
136 110,0.0,1.1,1.0);
137 histo->add1D("23","EM energy deposition in the target normalized to beam energy",
138 110,0.0,1.1,1.0);
139 histo->add1D("24","Pion energy deposition in the target normalized to beam energy",
140 110,0.0,1.1,1.0);
141 histo->add1D("25","Proton energy deposition in the target normalized to beam energy",
142 110,0.0,1.1,1.0);
143}
144
145//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
146
147void HistoManager::BeginOfRun()
148{
149 absZ0 = -0.5*length;
150 n_evt = 0;
151 n_elec = 0;
152 n_posit = 0;
153 n_gam = 0;
154 n_step = 0;
155 n_prot_leak = 0;
156 n_pion_leak = 0;
157 n_ions = 0;
158 n_deut = 0;
159 n_alpha = 0;
160 n_kaons = 0;
161 n_muons = 0;
162 n_cpions = 0;
163 n_pi0 = 0;
164 n_neutron = 0;
165 n_proton = 0;
166 n_aproton = 0;
167 n_neu_forw = 0;
168 n_neu_leak = 0;
169 n_neu_back = 0;
170
171 edepSum = 0.0;
172 edepSum2 = 0.0;
173
174 bookHisto();
175 histo->book();
176
177 if(verbose > 0)
178 G4cout << "HistoManager: Histograms are booked and run has been started"
179 <<G4endl<<" BeginOfRun (After histo->book)"<< G4endl;
180}
181
182//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183
184void HistoManager::EndOfRun()
185{
186
187 G4cout << "HistoManager: End of run actions are started" << G4endl;
188
189 // Average values
190 G4cout<<"========================================================"<<G4endl;
191
192 G4double x = (G4double)n_evt;
193 if(n_evt > 0) x = 1.0/x;
194
195 G4double xe = x*(G4double)n_elec;
196 G4double xg = x*(G4double)n_gam;
197 G4double xp = x*(G4double)n_posit;
198 G4double xs = x*(G4double)n_step;
199 G4double xn = x*(G4double)n_neutron;
200 G4double xpn = x*(G4double)n_proton;
201 G4double xap = x*(G4double)n_aproton;
202 G4double xnf = x*(G4double)n_neu_forw;
203 G4double xnb = x*(G4double)n_neu_leak;
204 G4double xnbw= x*(G4double)n_neu_back;
205 G4double xpl = x*(G4double)n_prot_leak;
206 G4double xal = x*(G4double)n_pion_leak;
207 G4double xpc = x*(G4double)n_cpions;
208 G4double xp0 = x*(G4double)n_pi0;
209 G4double xpk = x*(G4double)n_kaons;
210 G4double xpm = x*(G4double)n_muons;
211 G4double xid = x*(G4double)n_deut;
212 G4double xia = x*(G4double)n_alpha;
213 G4double xio = x*(G4double)n_ions;
214
215 edepSum *= x;
216 edepSum2 *= x;
217 edepSum2 -= edepSum*edepSum;
218 if(edepSum2 > 0.0) edepSum2 = std::sqrt(edepSum2);
219 else edepSum2 = 0.0;
220
221 G4cout << "Beam particle "
222 << primaryDef->GetParticleName() <<G4endl;
223 G4cout << "Beam Energy(MeV) "
224 << primaryKineticEnergy/MeV <<G4endl;
225 G4cout << "Number of events " << n_evt <<G4endl;
226 G4cout << std::setprecision(4) << "Average energy deposit (MeV) " << edepSum/MeV
227 << " RMS(MeV) " << edepSum2/MeV << G4endl;
228 G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
229 G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
230 G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
231 G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
232 G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
233 G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
234 G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
235 G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
236 G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
237 G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
238 G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
239 G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
240 G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
241 G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
242 G4cout << std::setprecision(4) << "Average number of forward neutrons " << xnf << G4endl;
243 G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl;
244 G4cout << std::setprecision(4) << "Average number of leaked neutrons " << xnbw << G4endl;
245 G4cout << std::setprecision(4) << "Average number of proton leak " << xpl << G4endl;
246 G4cout << std::setprecision(4) << "Average number of pion leak " << xal << G4endl;
247 G4cout<<"========================================================"<<G4endl;
248 G4cout<<G4endl;
249
250 // normalise histograms
251 for(G4int i=0; i<nHisto; i++) {
252 histo->scale(i,x);
253 }
254
255 if(verbose > 1) histo->print(0);
256 histo->save();
257}
258
259//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
260
261void HistoManager::BeginOfEvent()
262{
263 edepEvt = 0.0;
264 edepEM = 0.0;
265 edepPI = 0.0;
266 edepP = 0.0;
267}
268
269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270
271void HistoManager::EndOfEvent()
272{
273 edepSum += edepEvt;
274 edepSum2 += edepEvt*edepEvt;
275 histo->fill(21,edepEvt/primaryKineticEnergy,1.0);
276 histo->fill(22,edepEM/primaryKineticEnergy,1.0);
277 histo->fill(23,edepPI/primaryKineticEnergy,1.0);
278 histo->fill(24,edepP/primaryKineticEnergy,1.0);
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282
283void HistoManager::ScoreNewTrack(const G4Track* track)
284{
285 const G4ParticleDefinition* pd = track->GetDefinition();
286 G4String name = pd->GetParticleName();
287 G4double e = track->GetKineticEnergy();
288
289 // Primary track
290 if(0 == track->GetParentID()) {
291
292 n_evt++;
293 primaryKineticEnergy = e;
294 primaryDef = pd;
295 G4ThreeVector dir = track->GetMomentumDirection();
296 if(1 < verbose)
297 G4cout << "### Primary " << name
298 << " kinE(MeV)= " << e/MeV
299 << "; m(MeV)= " << pd->GetPDGMass()/MeV
300 << "; pos(mm)= " << track->GetPosition()/mm
301 << "; dir= " << track->GetMomentumDirection()
302 << G4endl;
303
304 // Secondary track
305 } else {
306 if(1 < verbose)
307 G4cout << "=== Secondary " << name
308 << " kinE(MeV)= " << e/MeV
309 << "; m(MeV)= " << pd->GetPDGMass()/MeV
310 << "; pos(mm)= " << track->GetPosition()/mm
311 << "; dir= " << track->GetMomentumDirection()
312 << G4endl;
313 e = std::log10(e/MeV);
314 if(pd == G4Gamma::Gamma()) {
315 n_gam++;
316 histo->fill(1,e,1.0);
317 } else if ( pd == G4Electron::Electron()) {
318 n_elec++;
319 histo->fill(2,e,1.0);
320 } else if ( pd == G4Positron::Positron()) {
321 n_posit++;
322 histo->fill(3,e,1.0);
323 } else if ( pd == G4Proton::Proton()) {
324 n_proton++;
325 histo->fill(4,e,1.0);
326 } else if ( pd == neutron) {
327 n_neutron++;
328 histo->fill(5,e,1.0);
329 } else if ( pd == G4AntiProton::AntiProton()) {
330 n_aproton++;
331 } else if ( pd == G4PionPlus::PionPlus() ) {
332 n_cpions++;
333 histo->fill(6,e,1.0);
334 histo->fill(19,e,1.0);
335
336 } else if ( pd == G4PionMinus::PionMinus()) {
337 n_cpions++;
338 histo->fill(6,e,1.0);
339 histo->fill(20,e,1.0);
340
341 } else if ( pd == G4PionZero::PionZero()) {
342 n_pi0++;
343 histo->fill(7,e,1.0);
344 } else if ( pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
345 n_kaons++;
346 histo->fill(8,e,1.0);
347 } else if ( pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
348 n_kaons++;
349 histo->fill(9,e,1.0);
350 } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
351 n_deut++;
352 histo->fill(10,e,1.0);
353 } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
354 n_alpha++;
355 histo->fill(11,e,1.0);
356 } else if ( pd->GetParticleType() == "nucleus") {
357 n_ions++;
358 histo->fill(12,e,1.0);
359 } else if ( pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
360 n_muons++;
361 histo->fill(13,e,1.0);
362 }
363 }
364}
365
366//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
367
368void HistoManager::AddTargetStep(const G4Step* step)
369{
370 n_step++;
371 G4double edep = step->GetTotalEnergyDeposit();
372 if(edep >= DBL_MIN) {
373 const G4Track* track = step->GetTrack();
374 currentDef = track->GetDefinition();
375 currentKinEnergy = track->GetKineticEnergy();
376
377 G4ThreeVector pos =
378 (step->GetPreStepPoint()->GetPosition() +
379 step->GetPostStepPoint()->GetPosition())*0.5;
380
381 G4double z = pos.z() - absZ0;
382
383 // scoring
384 edepEvt += edep;
385 histo->fill(0,z,edep);
386 const G4ParticleDefinition* pd = currentDef;
387
388 if(pd == G4Gamma::Gamma() || pd == G4Electron::Electron()
389 || pd == G4Positron::Positron()) {
390 edepEM += edep;
391 } else if ( pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
392 edepPI += edep;
393 } else if ( pd == G4Proton::Proton() || pd == G4AntiProton::AntiProton()) {
394 edepP += edep;
395 }
396
397 if(1 < verbose)
398 G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
399 << "; z(mm)= " << z/mm
400 << "; step(mm)= " << step->GetStepLength()/mm
401 << " by " << currentDef->GetParticleName()
402 << " E(MeV)= " << currentKinEnergy/MeV
403 << G4endl;
404 }
405}
406
407//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408
409void HistoManager::AddLeakingParticle(const G4Track* track)
410{
411 const G4ParticleDefinition* pd = track->GetDefinition();
412 G4double e = std::log10(track->GetKineticEnergy()/MeV);
413
414 G4ThreeVector pos = track->GetPosition();
415 G4ThreeVector dir = track->GetMomentumDirection();
416 G4double x = pos.x();
417 G4double y = pos.y();
418 G4double z = pos.z();
419
420 G4bool isLeaking = false;
421
422 // Forward
423 if(z > -absZ0 && dir.z() > 0.0) {
424 if(pd == neutron) {
425 n_neu_forw++;
426 histo->fill(15,e,1.0);
427 } else isLeaking = true;
428
429 // Backward
430 } else if (z < absZ0 && dir.z() < 0.0) {
431 if(pd == neutron) {
432 n_neu_leak++;
433 histo->fill(16,e,1.0);
434 } else isLeaking = true;
435
436 // Side
437 } else if (std::abs(z) <= -absZ0 && x*dir.x() + y*dir.y() > 0.0) {
438 isLeaking = true;
439 if(pd == neutron) {
440 n_neu_back++;
441 histo->fill(14,e,1.0);
442 } else isLeaking = true;
443 }
444
445 // protons and pions
446 if(isLeaking) {
447 if(pd == G4Proton::Proton()) {
448 histo->fill(17,e,1.0);
449 n_prot_leak++;
450 } else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
451 histo->fill(18,e,1.0);
452 n_pion_leak++;
453 }
454 }
455}
456
457//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
458
459void HistoManager::SetVerbose(G4int val)
460{
461 verbose = val;
462 histo->setVerbose(val);
463}
464
465//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466
467void HistoManager::SetTargetMaterial(const G4Material* mat)
468{
469 if(mat) {
470 material = mat;
471 elm = (*(material->GetElementVector()))[0];
472 }
473}
474
475//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476
477void HistoManager::Fill(G4int id, G4double x, G4double w)
478{
479 histo->fill(id, x, w);
480}
481
482//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
483
Note: See TracBrowser for help on using the repository browser.