Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 88154 2015-02-02 12:25:20Z 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 "G4RunManager.hh"
75 #include "DetectorConstruction.hh"
76 #include "G4NistManager.hh"
77 #include "G4Material.hh"
78 #include "G4BuilderType.hh"
79 #include "G4SystemOfUnits.hh"
80 
81 #include "IonUrQMDPhysics.hh"
82 #include "IonHIJINGPhysics.hh"
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
86 HistoManager* HistoManager::fManager = 0;
87 
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89 
91 {
92  if(!fManager) {
93  static HistoManager manager;
94  fManager = &manager;
95  }
96  return fManager;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
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 
152 {
153  delete fHisto;
154 }
155 
156 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
157 
159 {
160  fHisto->Add1D("0","Energy deposition (MeV/mm/event) in the target",
161  fNSlices,0.0,fLength/mm,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",
172  fNBinsE,-5.,5.,1.0);
173  fHisto->Add1D("11","Log10 Energy (GeV) of He3 and alpha",fNBinsE,-5.,5.,1.0);
174  fHisto->Add1D("12","Log10 Energy (GeV) of Generic Ions",fNBinsE,-5.,5.,1.0);
175  fHisto->Add1D("13","Log10 Energy (GeV) of muons",fNBinsE,-4.,6.,1.0);
176  fHisto->Add1D("14","Log10 Energy (GeV) of pi+",fNBinsE,-4.,6.,1.0);
177  fHisto->Add1D("15","Log10 Energy (GeV) of pi-",fNBinsE,-4.,6.,1.0);
178  fHisto->Add1D("16","X Section (mb) of Secondary Fragments Z with E>1 GeV (mb)"
179  ,25,0.5,25.5,1.0);
180  fHisto->Add1D("17","Secondary Fragment A E>1 GeV",50,0.5,50.5,1.0);
181  fHisto->Add1D("18","Secondary Fragment Z E<1 GeV",25,0.5,25.5,1.0);
182  fHisto->Add1D("19","Secondary Fragment A E<1 GeV",50,0.5,50.5,1.0);
183  fHisto->Add1D("20","X Section (mb) of Secondary Fragments Z (mb) ",
184  25,0.5,25.5,1.0);
185  fHisto->Add1D("21","Secondary Fragment A ",50,0.5,50.5,1.0);
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
189 
191 {
192  Initialise();
193  BookHisto();
194  fHisto->Book();
195 
196  if(fVerbose > 0) {
197  G4cout << "HistoManager: Histograms are booked and run has been started"
198  <<G4endl<<" BeginOfRun (After fHisto->book)"<< G4endl;
199  }
200 }
201 
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203 
205 {
206 
207  G4cout << "HistoManager: End of run actions are started" << G4endl;
208 
209  // Average values
210  G4cout<<"========================================================"<<G4endl;
211 
212  G4double x = (G4double)fNevt;
213  if(fNevt > 0) { x = 1.0/x; }
214 
215  G4double xe = x*fNelec;
216  G4double xg = x*fNgam;
217  G4double xp = x*fNposit;
218  G4double xs = x*fNstep;
219  G4double xn = x*fNneutron;
220  G4double xpn = x*fNproton;
221  G4double xap = x*fNaproton;
222  G4double xpc = x*fNcpions;
223  G4double xp0 = x*fNpi0;
224  G4double xpk = x*fNkaons;
225  G4double xpm = x*fNmuons;
226  G4double xid = x*fNdeut;
227  G4double xia = x*fNalpha;
228  G4double xio = x*fNions;
229 
230  fEdepSum *= x;
231  fEdepSum2 *= x;
232  fEdepSum2 -= fEdepSum*fEdepSum;
233  if(fEdepSum2 > 0.0) fEdepSum2 = std::sqrt(fEdepSum2);
234  else fEdepSum2 = 0.0;
235 
236  G4cout << "Beam particle "
237  << fPrimaryDef->GetParticleName() <<G4endl;
238  G4cout << "Beam Energy(GeV) "
239  << fPrimaryKineticEnergy/GeV <<G4endl;
240  G4cout << "Number of events "
241  << fNevt <<G4endl;
242  G4cout << std::setprecision(4) << "Average energy deposit (GeV) "
243  << fEdepSum/GeV
244  << " RMS(GeV) " << fEdepSum2/GeV << G4endl;
245  G4cout << std::setprecision(4) << "Average number of steps "
246  << xs << G4endl;
247  G4cout << std::setprecision(4) << "Average number of gamma "
248  << xg << G4endl;
249  G4cout << std::setprecision(4) << "Average number of e- "
250  << xe << G4endl;
251  G4cout << std::setprecision(4) << "Average number of e+ "
252  << xp << G4endl;
253  G4cout << std::setprecision(4) << "Average number of neutrons "
254  << xn << G4endl;
255  G4cout << std::setprecision(4) << "Average number of protons "
256  << xpn << G4endl;
257  G4cout << std::setprecision(4) << "Average number of antiprotons "
258  << xap << G4endl;
259  G4cout << std::setprecision(4) << "Average number of pi+ & pi- "
260  << xpc << G4endl;
261  G4cout << std::setprecision(4) << "Average number of pi0 "
262  << xp0 << G4endl;
263  G4cout << std::setprecision(4) << "Average number of kaons "
264  << xpk << G4endl;
265  G4cout << std::setprecision(4) << "Average number of muons "
266  << xpm << G4endl;
267  G4cout << std::setprecision(4) << "Average number of deuterons+tritons "
268  << xid << G4endl;
269  G4cout << std::setprecision(4) << "Average number of He3+alpha "
270  << xia << G4endl;
271  G4cout << std::setprecision(4) << "Average number of ions "
272  << xio << G4endl;
273  G4cout<<"========================================================"<<G4endl;
274  G4cout<<G4endl;
275 
276  // normalise histograms
277  for(G4int i=0; i<fNHisto; i++) {
278  fHisto->ScaleH1(i,x);
279  }
280  // will work only for pure material - 1 element
281  // G4cout << fMaterial << G4endl;
282  G4double F= 1000/(fLength*barn*fMaterial->GetTotNbOfAtomsPerVolume());
283  if(F > 0.0) {
284  fHisto->ScaleH1(16, F);
285  fHisto->ScaleH1(20, F);
286  }
287 
288  fHisto->Save();
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292 
294 {
295  ++fNevt;
296  fEdepEvt = 0.0;
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
300 
302 {
303  fEdepSum += fEdepEvt;
304  fEdepSum2 += fEdepEvt*fEdepEvt;
305 }
306 
307 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
308 
309 void HistoManager::ScoreNewTrack(const G4Track* track)
310 {
311  const G4ParticleDefinition* pd = track->GetDefinition();
312  G4String name = pd->GetParticleName();
313  G4double e = track->GetKineticEnergy();
314 
315  // Primary track
316  if(0 == track->GetParentID()) {
317 
318  fPrimaryKineticEnergy = e;
319  fPrimaryDef = pd;
320  G4ThreeVector dir = track->GetMomentumDirection();
321  if(1 < fVerbose)
322  G4cout << "### Primary " << name
323  << " kinE(GeV)= " << e/GeV
324  << "; m(GeV)= " << pd->GetPDGMass()/GeV
325  << "; pos(mm)= " << track->GetPosition()/mm
326  << "; dir= " << track->GetMomentumDirection()
327  << G4endl;
328 
329  // Secondary track
330  } else {
331  if(1 < fVerbose) {
332  G4cout << "=== Secondary " << name
333  << " kinE(GeV)= " << e/GeV
334  << "; m(GeV)= " << pd->GetPDGMass()/GeV
335  << "; pos(mm)= " << track->GetPosition()/mm
336  << "; dir= " << track->GetMomentumDirection()
337  << G4endl;
338  }
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;
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;
411  G4double edep = step->GetTotalEnergyDeposit();
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 == "HIJING") {
459 #ifdef G4_USE_HIJING
460  fIonPhysics = new IonHIJINGPhysics();
461  fPhysList->ReplacePhysics(fIonPhysics);
463  G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added"
464  << G4endl;
465 #else
466  G4cout << "Error: Ion Physics HIJING is requested but is not available"
467  <<G4endl;
468 #endif
469  } else if(nam == "UrQMD") {
470 #ifdef G4_USE_URQMD
471  fIonPhysics = new IonUrQMDPhysics(1);
472  fPhysList->ReplacePhysics(fIonPhysics);
474  G4cout << "### SetIonPhysics: Ion Physics UrQMD is added"
475  << G4endl;
476 #else
477  G4cout << "Error: Ion Physics UrQMD is requested but is not available"
478  <<G4endl;
479 #endif
480  } else {
481  G4cout << "### HistoManager WARNING: Ion Physics <"
482  << nam << "> is unknown!" << G4endl;
483  }
484 }
485 
486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487 
G4ParticleDefinition * GetDefinition() const
const XML_Char * name
Definition: expat.h:151
G4int GetParentID() const
void AddTargetStep(const G4Step *)
static G4MuonPlus * MuonPlus()
Definition: G4MuonPlus.cc:99
static constexpr double mm
Definition: G4SIunits.hh:115
Definition of the IonHIJINGPhysics class.
void BookHisto()
G4double GetStepLength() const
void EndOfEvent()
Definition: Histo.hh:56
Definition of the IonUrQMDPhysics class.
const G4ThreeVector & GetPosition() const
static G4KaonZeroLong * KaonZeroLong()
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
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
double A(double temperature)
const G4ThreeVector & GetPosition() const
void PhysicsHasBeenModified()
static G4KaonZeroShort * KaonZeroShort()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
void BeginOfEvent()
static constexpr double eplus
Definition: G4SIunits.hh:199
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:108
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4double GetTotalEnergyDeposit() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
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
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
void ScoreNewTrack(const G4Track *aTrack)
void SetVerbose(G4int val)
Definition: HistoManager.hh:95
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 constexpr double MeV
Definition: G4SIunits.hh:214
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static HistoManager * GetPointer()
Definition: HistoManager.cc:64
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4Track * GetTrack() const
static constexpr double barn
Definition: G4SIunits.hh:105
G4double GetPDGCharge() const
static G4He3 * He3()
Definition: G4He3.cc:94
static constexpr double keV
Definition: G4SIunits.hh:216
static const G4double pos