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