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

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

update
File size: 15.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.14 2007/05/24 13:52:31 vnivanch Exp $
27// GEANT4 tag $Name:  $
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    = 22;
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 (GeV) in the target",
136               nBinsE,0.0,edepMax,GeV);
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141void HistoManager::BeginOfRun()
142{
143  absZ0       = -0.5*length;
144  n_evt       = 0;
145  n_elec      = 0;
146  n_posit     = 0;
147  n_gam       = 0;
148  n_step      = 0;
149  n_prot_leak = 0;
150  n_pion_leak = 0;
151  n_ions      = 0;
152  n_deut      = 0;
153  n_alpha     = 0;
154  n_kaons     = 0;
155  n_muons     = 0;
156  n_cpions    = 0;
157  n_pi0       = 0;
158  n_neutron   = 0;
159  n_proton    = 0;
160  n_aproton   = 0;
161  n_neu_forw  = 0;
162  n_neu_leak  = 0;
163  n_neu_back  = 0;
164
165  edepSum     = 0.0;
166  edepSum2    = 0.0;
167
168  bookHisto();
169  histo->book();
170
171  if(verbose > 0) 
172    G4cout << "HistoManager: Histograms are booked and run has been started"
173           <<G4endl<<"  BeginOfRun (After histo->book)"<< G4endl;
174}
175
176//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
177
178void HistoManager::EndOfRun()
179{
180
181  G4cout << "HistoManager: End of run actions are started" << G4endl;
182
183  // Average values
184  G4cout<<"========================================================"<<G4endl;
185
186  G4double x = (G4double)n_evt;
187  if(n_evt > 0) x = 1.0/x;
188
189  G4double xe = x*(G4double)n_elec;
190  G4double xg = x*(G4double)n_gam;
191  G4double xp = x*(G4double)n_posit;
192  G4double xs = x*(G4double)n_step;
193  G4double xn = x*(G4double)n_neutron;
194  G4double xpn = x*(G4double)n_proton;
195  G4double xap = x*(G4double)n_aproton;
196  G4double xnf = x*(G4double)n_neu_forw;
197  G4double xnb = x*(G4double)n_neu_leak;
198  G4double xnbw= x*(G4double)n_neu_back;
199  G4double xpl = x*(G4double)n_prot_leak;
200  G4double xal = x*(G4double)n_pion_leak;
201  G4double xpc = x*(G4double)n_cpions;
202  G4double xp0 = x*(G4double)n_pi0;
203  G4double xpk = x*(G4double)n_kaons;
204  G4double xpm = x*(G4double)n_muons;
205  G4double xid = x*(G4double)n_deut;
206  G4double xia = x*(G4double)n_alpha;
207  G4double xio = x*(G4double)n_ions;
208
209  edepSum  *= x;
210  edepSum2 *= x;
211  edepSum2 -= edepSum*edepSum;
212  if(edepSum2 > 0.0) edepSum2 = std::sqrt(edepSum2);
213  else               edepSum2 = 0.0;
214
215  G4cout                         << "Beam particle                        "
216                                 << primaryDef->GetParticleName() <<G4endl;
217  G4cout                         << "Beam Energy(MeV)                     " 
218                                 << primaryKineticEnergy/MeV <<G4endl;
219  G4cout                         << "Number of events                     " << n_evt <<G4endl;
220  G4cout << std::setprecision(4) << "Average energy deposit (MeV)         " << edepSum/MeV
221         << "   RMS(MeV) " << edepSum2/MeV << G4endl;
222  G4cout << std::setprecision(4) << "Average number of steps              " << xs << G4endl;
223  G4cout << std::setprecision(4) << "Average number of gamma              " << xg << G4endl;
224  G4cout << std::setprecision(4) << "Average number of e-                 " << xe << G4endl;
225  G4cout << std::setprecision(4) << "Average number of e+                 " << xp << G4endl;
226  G4cout << std::setprecision(4) << "Average number of neutrons           " << xn << G4endl;
227  G4cout << std::setprecision(4) << "Average number of protons            " << xpn << G4endl;
228  G4cout << std::setprecision(4) << "Average number of antiprotons        " << xap << G4endl;
229  G4cout << std::setprecision(4) << "Average number of pi+ & pi-          " << xpc << G4endl;
230  G4cout << std::setprecision(4) << "Average number of pi0                " << xp0 << G4endl;
231  G4cout << std::setprecision(4) << "Average number of kaons              " << xpk << G4endl;
232  G4cout << std::setprecision(4) << "Average number of muons              " << xpm << G4endl;
233  G4cout << std::setprecision(4) << "Average number of deuterons+tritons  " << xid << G4endl;
234  G4cout << std::setprecision(4) << "Average number of He3+alpha          " << xia << G4endl;
235  G4cout << std::setprecision(4) << "Average number of ions               " << xio << G4endl;
236  G4cout << std::setprecision(4) << "Average number of forward neutrons   " << xnf << G4endl;
237  G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl;
238  G4cout << std::setprecision(4) << "Average number of leaked neutrons    " << xnbw << G4endl;
239  G4cout << std::setprecision(4) << "Average number of proton leak        " << xpl << G4endl;
240  G4cout << std::setprecision(4) << "Average number of pion leak          " << xal << G4endl;
241  G4cout<<"========================================================"<<G4endl;
242  G4cout<<G4endl;
243
244  // normalise histograms
245  for(G4int i=0; i<nHisto; i++) {
246    histo->scale(i,x);
247  }
248
249  if(verbose > 1) histo->print(0);
250  histo->save();
251}
252
253//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
254
255void HistoManager::BeginOfEvent()
256{
257  edepEvt = 0.0;
258}
259
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
261
262void HistoManager::EndOfEvent()
263{
264  edepSum  += edepEvt;
265  edepSum2 += edepEvt*edepEvt;
266  histo->fill(21,edepEvt,1.0);
267}
268
269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
270
271void HistoManager::ScoreNewTrack(const G4Track* track)
272{
273  const G4ParticleDefinition* pd = track->GetDefinition();
274  G4String name = pd->GetParticleName();
275  G4double e = track->GetKineticEnergy();
276
277  // Primary track
278  if(0 == track->GetParentID()) {
279
280    n_evt++;
281    primaryKineticEnergy = e;
282    primaryDef = pd;
283    G4ThreeVector dir = track->GetMomentumDirection();
284    if(1 < verbose) 
285      G4cout << "### Primary " << name
286             << " kinE(MeV)= " << e/MeV
287             << "; m(MeV)= " << pd->GetPDGMass()/MeV
288             << "; pos(mm)= " << track->GetPosition()/mm
289             << ";  dir= " << track->GetMomentumDirection() 
290             << G4endl;
291
292    // Secondary track
293  } else {
294    if(1 < verbose) 
295      G4cout << "=== Secondary " << name
296             << " kinE(MeV)= " << e/MeV
297             << "; m(MeV)= " << pd->GetPDGMass()/MeV
298             << "; pos(mm)= " << track->GetPosition()/mm
299             << ";  dir= " << track->GetMomentumDirection() 
300             << G4endl;
301    e = std::log10(e/MeV);
302    if(pd == G4Gamma::Gamma()) {
303      n_gam++;
304      histo->fill(1,e,1.0);
305    } else if ( pd == G4Electron::Electron()) {
306      n_elec++;
307      histo->fill(2,e,1.0);
308    } else if ( pd == G4Positron::Positron()) {
309      n_posit++;
310      histo->fill(3,e,1.0);
311    } else if ( pd == G4Proton::Proton()) {
312      n_proton++;
313      histo->fill(4,e,1.0);
314    } else if ( pd == neutron) {
315      n_neutron++;
316      histo->fill(5,e,1.0);
317    } else if ( pd == G4AntiProton::AntiProton()) {
318      n_aproton++;
319    } else if ( pd == G4PionPlus::PionPlus() ) {
320      n_cpions++;
321      histo->fill(6,e,1.0);
322      histo->fill(19,e,1.0);
323
324    } else if ( pd == G4PionMinus::PionMinus()) {
325      n_cpions++;
326      histo->fill(6,e,1.0);
327      histo->fill(20,e,1.0);
328
329    } else if ( pd == G4PionZero::PionZero()) {
330      n_pi0++;
331      histo->fill(7,e,1.0);
332    } else if ( pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
333      n_kaons++;
334      histo->fill(8,e,1.0);
335    } else if ( pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
336      n_kaons++;
337      histo->fill(9,e,1.0);
338    } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
339      n_deut++;
340      histo->fill(10,e,1.0);
341    } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
342      n_alpha++;
343      histo->fill(11,e,1.0);
344    } else if ( pd->GetParticleType() == "nucleus") {
345      n_ions++;
346      histo->fill(12,e,1.0);
347    } else if ( pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
348      n_muons++;
349      histo->fill(13,e,1.0);   
350    }
351  }
352}
353
354//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
355
356void HistoManager::AddTargetStep(const G4Step* step)
357{
358  n_step++;
359  G4double edep = step->GetTotalEnergyDeposit();
360  if(edep >= DBL_MIN) { 
361    const G4Track* track = step->GetTrack();
362    currentDef = track->GetDefinition(); 
363    currentKinEnergy = track->GetKineticEnergy();
364
365    G4ThreeVector pos = 
366      (step->GetPreStepPoint()->GetPosition() +
367       step->GetPostStepPoint()->GetPosition())*0.5;
368
369    G4double z = pos.z() - absZ0;
370
371    // scoring
372    edepEvt += edep;
373    histo->fill(0,z,edep);
374
375    if(1 < verbose) 
376      G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
377             << "; z(mm)= " << z/mm
378             << "; step(mm)= " << step->GetStepLength()/mm
379             << " by " << currentDef->GetParticleName()
380             << " E(MeV)= " << currentKinEnergy/MeV
381             << G4endl;
382  }
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387void HistoManager::AddLeakingParticle(const G4Track* track)
388{
389  const G4ParticleDefinition* pd = track->GetDefinition(); 
390  G4double e = std::log10(track->GetKineticEnergy()/MeV);
391
392  G4ThreeVector pos = track->GetPosition();
393  G4ThreeVector dir = track->GetMomentumDirection();
394  G4double x = pos.x();
395  G4double y = pos.y();
396  G4double z = pos.z();
397 
398  G4bool isLeaking = false;
399
400  // Forward
401  if(z > -absZ0 && dir.z() > 0.0) {
402    if(pd == neutron) {
403      n_neu_forw++;
404      histo->fill(15,e,1.0);
405    } else isLeaking = true;
406
407    // Backward
408  } else if (z < absZ0 && dir.z() < 0.0) {
409    if(pd == neutron) {
410      n_neu_leak++;
411      histo->fill(16,e,1.0);
412    } else isLeaking = true;
413
414    // Side
415  } else if (std::abs(z) <= -absZ0 && x*dir.x() + y*dir.y() > 0.0) {
416    isLeaking = true;
417    if(pd == neutron) {
418      n_neu_back++;
419      histo->fill(14,e,1.0);
420    } else isLeaking = true;
421  }
422
423  // protons and pions
424  if(isLeaking) {
425    if(pd == G4Proton::Proton()) {
426      histo->fill(17,e,1.0);
427      n_prot_leak++;
428    } else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
429      histo->fill(18,e,1.0);
430      n_pion_leak++;
431    }
432  }
433}
434
435//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
436
437void HistoManager::SetVerbose(G4int val)       
438{
439  verbose = val; 
440  histo->setVerbose(val);
441}
442
443//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444
445void HistoManager::SetTargetMaterial(const G4Material* mat)         
446{
447  if(mat) {
448    material = mat;
449    elm = (*(material->GetElementVector()))[0];
450  }
451}
452
453//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
454
455void HistoManager::Fill(G4int id, G4double x, G4double w)
456{
457  histo->fill(id, x, w);
458}
459
460//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
461
Note: See TracBrowser for help on using the repository browser.