Geant4  9.6.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$
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 
91 HistoManager* HistoManager::GetPointer()
92 {
93  if(!fManager) {
94  static HistoManager manager;
95  fManager = &manager;
96  }
97  return fManager;
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
102 HistoManager::HistoManager()
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 
152 HistoManager::~HistoManager()
153 {
154  delete fHisto;
155 }
156 
157 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
158 
159 void HistoManager::BookHisto()
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",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)",25,0.5,25.5,1.0);
179  fHisto->Add1D("17","Secondary Fragment A E>1 GeV",50,0.5,50.5,1.0);
180  fHisto->Add1D("18","Secondary Fragment Z E<1 GeV",25,0.5,25.5,1.0);
181  fHisto->Add1D("19","Secondary Fragment A E<1 GeV",50,0.5,50.5,1.0);
182  fHisto->Add1D("20","X Section (mb) of Secondary Fragments Z (mb) ",25,0.5,25.5,1.0);
183  fHisto->Add1D("21","Secondary Fragment A ",50,0.5,50.5,1.0);
184 }
185 
186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
187 
188 void HistoManager::BeginOfRun()
189 {
190  Initialise();
191  BookHisto();
192  fHisto->Book();
193 
194  if(fVerbose > 0) {
195  G4cout << "HistoManager: Histograms are booked and run has been started"
196  <<G4endl<<" BeginOfRun (After fHisto->book)"<< G4endl;
197  }
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201 
202 void HistoManager::EndOfRun()
203 {
204 
205  G4cout << "HistoManager: End of run actions are started" << G4endl;
206 
207  // Average values
208  G4cout<<"========================================================"<<G4endl;
209 
210  G4double x = (G4double)fNevt;
211  if(fNevt > 0) { x = 1.0/x; }
212 
213  G4double xe = x*fNelec;
214  G4double xg = x*fNgam;
215  G4double xp = x*fNposit;
216  G4double xs = x*fNstep;
217  G4double xn = x*fNneutron;
218  G4double xpn = x*fNproton;
219  G4double xap = x*fNaproton;
220  G4double xpc = x*fNcpions;
221  G4double xp0 = x*fNpi0;
222  G4double xpk = x*fNkaons;
223  G4double xpm = x*fNmuons;
224  G4double xid = x*fNdeut;
225  G4double xia = x*fNalpha;
226  G4double xio = x*fNions;
227 
228  fEdepSum *= x;
229  fEdepSum2 *= x;
230  fEdepSum2 -= fEdepSum*fEdepSum;
231  if(fEdepSum2 > 0.0) fEdepSum2 = std::sqrt(fEdepSum2);
232  else fEdepSum2 = 0.0;
233 
234  G4cout << "Beam particle "
235  << fPrimaryDef->GetParticleName() <<G4endl;
236  G4cout << "Beam Energy(GeV) "
237  << fPrimaryKineticEnergy/GeV <<G4endl;
238  G4cout << "Number of events " << fNevt <<G4endl;
239  G4cout << std::setprecision(4) << "Average energy deposit (GeV) " << fEdepSum/GeV
240  << " RMS(GeV) " << fEdepSum2/GeV << G4endl;
241  G4cout << std::setprecision(4) << "Average number of steps " << xs << G4endl;
242  G4cout << std::setprecision(4) << "Average number of gamma " << xg << G4endl;
243  G4cout << std::setprecision(4) << "Average number of e- " << xe << G4endl;
244  G4cout << std::setprecision(4) << "Average number of e+ " << xp << G4endl;
245  G4cout << std::setprecision(4) << "Average number of neutrons " << xn << G4endl;
246  G4cout << std::setprecision(4) << "Average number of protons " << xpn << G4endl;
247  G4cout << std::setprecision(4) << "Average number of antiprotons " << xap << G4endl;
248  G4cout << std::setprecision(4) << "Average number of pi+ & pi- " << xpc << G4endl;
249  G4cout << std::setprecision(4) << "Average number of pi0 " << xp0 << G4endl;
250  G4cout << std::setprecision(4) << "Average number of kaons " << xpk << G4endl;
251  G4cout << std::setprecision(4) << "Average number of muons " << xpm << G4endl;
252  G4cout << std::setprecision(4) << "Average number of deuterons+tritons " << xid << G4endl;
253  G4cout << std::setprecision(4) << "Average number of He3+alpha " << xia << G4endl;
254  G4cout << std::setprecision(4) << "Average number of ions " << xio << G4endl;
255  G4cout<<"========================================================"<<G4endl;
256  G4cout<<G4endl;
257 
258  // normalise histograms
259  for(G4int i=0; i<fNHisto; i++) {
260  fHisto->ScaleH1(i,x);
261  }
262  // will work only for pure material - 1 element
263  // G4cout << fMaterial << G4endl;
264  G4double F= 1000/(fLength*barn*fMaterial->GetTotNbOfAtomsPerVolume());
265  if(F > 0.0) {
266  fHisto->ScaleH1(16, F);
267  fHisto->ScaleH1(20, F);
268  }
269 
270  fHisto->Save();
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
274 
275 void HistoManager::BeginOfEvent()
276 {
277  ++fNevt;
278  fEdepEvt = 0.0;
279 }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 
283 void HistoManager::EndOfEvent()
284 {
285  fEdepSum += fEdepEvt;
286  fEdepSum2 += fEdepEvt*fEdepEvt;
287 }
288 
289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
290 
291 void HistoManager::ScoreNewTrack(const G4Track* track)
292 {
293  const G4ParticleDefinition* pd = track->GetDefinition();
294  G4String name = pd->GetParticleName();
295  G4double e = track->GetKineticEnergy();
296 
297  // Primary track
298  if(0 == track->GetParentID()) {
299 
300  fPrimaryKineticEnergy = e;
301  fPrimaryDef = pd;
302  G4ThreeVector dir = track->GetMomentumDirection();
303  if(1 < fVerbose)
304  G4cout << "### Primary " << name
305  << " kinE(GeV)= " << e/GeV
306  << "; m(GeV)= " << pd->GetPDGMass()/GeV
307  << "; pos(mm)= " << track->GetPosition()/mm
308  << "; dir= " << track->GetMomentumDirection()
309  << G4endl;
310 
311  // Secondary track
312  } else {
313  if(1 < fVerbose)
314  G4cout << "=== Secondary " << name
315  << " kinE(GeV)= " << e/GeV
316  << "; m(GeV)= " << pd->GetPDGMass()/GeV
317  << "; pos(mm)= " << track->GetPosition()/mm
318  << "; dir= " << track->GetMomentumDirection()
319  << G4endl;
320  e = std::log10(e/GeV);
321  if(pd == G4Gamma::Gamma()) {
322  fNgam++;
323  fHisto->Fill(1,e,1.0);
324  } else if ( pd == G4Electron::Electron()) {
325  fNelec++;
326  fHisto->Fill(2,e,1.0);
327  } else if ( pd == G4Positron::Positron()) {
328  fNposit++;
329  fHisto->Fill(3,e,1.0);
330  } else if ( pd == G4Proton::Proton()) {
331  fNproton++;
332  fHisto->Fill(4,e,1.0);
333  } else if ( pd == fNeutron) {
334  fNneutron++;
335  fHisto->Fill(5,e,1.0);
336  } else if ( pd == G4AntiProton::AntiProton()) {
337  fNaproton++;
338  } else if ( pd == G4PionPlus::PionPlus() ) {
339  fNcpions++;
340  fHisto->Fill(6,e,1.0);
341  fHisto->Fill(14,e,1.0);
342 
343  } else if ( pd == G4PionMinus::PionMinus()) {
344  fNcpions++;
345  fHisto->Fill(6,e,1.0);
346  fHisto->Fill(15,e,1.0);
347 
348  } else if ( pd == G4PionZero::PionZero()) {
349  fNpi0++;
350  fHisto->Fill(7,e,1.0);
351  } else if ( pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
352  fNkaons++;
353  fHisto->Fill(8,e,1.0);
354  } else if ( pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
355  fNkaons++;
356  fHisto->Fill(9,e,1.0);
357  } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
358  fNdeut++;
359  fHisto->Fill(10,e,1.0);
360  } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
361  fNalpha++;
362  fHisto->Fill(11,e,1.0);
363  } else if ( pd->GetParticleType() == "nucleus") {
364  fNions++;
365  fHisto->Fill(12,e,1.0);
366  G4double Z = pd->GetPDGCharge()/eplus;
367  G4double A = (G4double)pd->GetBaryonNumber();
368  if(e > 0.0) {
369  fHisto->Fill(16,Z,1.0);
370  fHisto->Fill(17,A,1.0);
371  } else {
372  fHisto->Fill(18,Z,1.0);
373  fHisto->Fill(19,A,1.0);
374  }
375  fHisto->Fill(20,Z,1.0);
376  fHisto->Fill(21,A,1.0);
377  } else if ( pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
378  fNmuons++;
379  fHisto->Fill(13,e,1.0);
380  }
381  }
382 }
383 
384 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
385 
386 void HistoManager::AddTargetStep(const G4Step* step)
387 {
388  ++fNstep;
389  G4double edep = step->GetTotalEnergyDeposit();
390 
391  if(edep > 0.0) {
392  G4ThreeVector pos =
393  (step->GetPreStepPoint()->GetPosition() +
394  step->GetPostStepPoint()->GetPosition())*0.5;
395 
396  G4double z = pos.z() - fAbsZ0;
397 
398  // scoring
399  fEdepEvt += edep;
400  fHisto->Fill(0,z,edep);
401 
402  if(1 < fVerbose) {
403  const G4Track* track = step->GetTrack();
404  G4cout << "HistoManager::AddEnergy: e(keV)= " << edep/keV
405  << "; z(mm)= " << z/mm
406  << "; step(mm)= " << step->GetStepLength()/mm
407  << " by " << track->GetDefinition()->GetParticleName()
408  << " E(GeV)= " << track->GetKineticEnergy()/GeV
409  << G4endl;
410  }
411  }
412 }
413 
414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
415 
416 void HistoManager::SetVerbose(G4int val)
417 {
418  fVerbose = val;
419  fHisto->SetVerbose(val);
420 }
421 
422 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
423 
424 void HistoManager::Fill(G4int id, G4double x, G4double w)
425 {
426  fHisto->Fill(id, x, w);
427 }
428 
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
430 
431 void HistoManager::SetIonPhysics(const G4String& nam)
432 {
433  if(fIonPhysics) {
434  G4cout << "### HistoManager WARNING: Ion Physics is already defined: <"
435  << nam << "> is ignored!" << G4endl;
436  } else if(nam == "DPMJET") {
437  fIonPhysics = new IonDPMJETPhysics(false);
438  fPhysList->ReplacePhysics(fIonPhysics);
439  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
440  G4cout << "### SetIonPhysics: Ion Physics DPMJET/Binary is added"
441  << G4endl;
442  } else if(nam == "FTF") {
443  fIonPhysics = new G4IonPhysics();
444  fPhysList->ReplacePhysics(fIonPhysics);
445  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
446  G4cout << "### SetIonPhysics: Ion Physics FTFP/Binary is added"
447  << G4endl;
448  } else if(nam == "UrQMD") {
449 #ifdef G4_USE_URQMD
450  fIonPhysics = new IonUrQMDPhysics(1);
451  fPhysList->ReplacePhysics(fIonPhysics);
452  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
453  G4cout << "### SetIonPhysics: Ion Physics UrQMD is added"
454  << G4endl;
455 #else
456  G4cout << "Error: Ion Physics UrQMD is requested but is not available"
457  <<G4endl;
458 #endif
459  } else {
460  G4cout << "### HistoManager WARNING: Ion Physics <"
461  << nam << "> is unknown!" << G4endl;
462  }
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
466 
Generated on Sat May 25 2013 14:32:04 for Geant4 by   doxygen 1.8.4