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 //
41 //----------------------------------------------------------------------------
42 //
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
47 #include "HistoManager.hh"
48 #include "G4UnitsTable.hh"
49 #include "G4Neutron.hh"
50 #include "globals.hh"
51 #include "G4ios.hh"
52 #include "G4ParticleDefinition.hh"
53 #include "G4ParticleTable.hh"
54 #include "G4NistManager.hh"
55 #include "G4HadronicProcessStore.hh"
56 
57 #include "G4NucleiProperties.hh"
58 #include "G4NistManager.hh"
59 #include "G4StableIsotopes.hh"
60 #include "G4SystemOfUnits.hh"
61 
62 #include "Histo.hh"
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
66 HistoManager* HistoManager::fManager = 0;
67 
68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
69 
70 HistoManager* HistoManager::GetPointer()
71 {
72  if(!fManager) {
73  static HistoManager manager;
74  fManager = &manager;
75  }
76  return fManager;
77 }
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 HistoManager::HistoManager()
82 {
83  fHisto = new Histo();
84  fNeutron = G4Neutron::Neutron();
85  fVerbose = 1;
86 
87  fParticleName = "proton";
88  fElementName = "Al";
89 
90  fMinKinEnergy = 0.1*MeV;
91  fMaxKinEnergy = 10*TeV;
92  fMinMomentum = 1*MeV;
93  fMaxMomentum = 10*TeV;
94 
95  fBinsE = 800;
96  fBinsP = 700;
97 
98  fHisto->SetVerbose(fVerbose);
99 
100  fIsInitialised = false;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 
105 HistoManager::~HistoManager()
106 {
107  delete fHisto;
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
112 void HistoManager::BeginOfRun()
113 {
114  G4double p1 = std::log10(fMinMomentum/GeV);
115  G4double p2 = std::log10(fMaxMomentum/GeV);
116  G4double e1 = std::log10(fMinKinEnergy/MeV);
117  G4double e2 = std::log10(fMaxKinEnergy/MeV);
118 
119  //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
120 
121  if(fIsInitialised) {
122  fHisto->SetHisto1D(0,fBinsP,p1,p2,1.0);
123  fHisto->SetHisto1D(1,fBinsE,e1,e2,1.0);
124  fHisto->SetHisto1D(2,fBinsP,p1,p2,1.0);
125  fHisto->SetHisto1D(3,fBinsE,e1,e2,1.0);
126  fHisto->SetHisto1D(4,fBinsE,e1,e2,1.0);
127  fHisto->SetHisto1D(5,fBinsE,e1,e2,1.0);
128  fHisto->SetHisto1D(6,fBinsE,e1,e2,1.0);
129  fHisto->SetHisto1D(7,fBinsE,e1,e2,1.0);
130 
131  } else {
132  fHisto->Add1D("h1","Elastic cross section (barn,1.0) as a functions of log10(p/GeV)",
133  fBinsP,p1,p2,1.0);
134  fHisto->Add1D("h2","Elastic cross section (barn) as a functions of log10(E/MeV)",
135  fBinsE,e1,e2,1.0);
136  fHisto->Add1D("h3","Inelastic cross section (barn) as a functions of log10(p/GeV)",
137  fBinsP,p1,p2,1.0);
138  fHisto->Add1D("h4","Inelastic cross section (barn) as a functions of log10(E/MeV)",
139  fBinsE,e1,e2,1.0);
140  fHisto->Add1D("h5","Capture cross section (barn) as a functions of log10(E/MeV)",
141  fBinsE,e1,e2,1.0);
142  fHisto->Add1D("h6","Fission cross section (barn) as a functions of log10(E/MeV)",
143  fBinsE,e1,e2,1.0);
144  fHisto->Add1D("h7","Charge exchange cross section (barn) as a functions of log10(E/MeV)",
145  fBinsE,e1,e2,1.0);
146  fHisto->Add1D("h8","Total cross section (barn) as a functions of log10(E/MeV)",
147  fBinsE,e1,e2,1.0);
148  }
149 
150  fIsInitialised = true;
151  fHisto->Book();
152 }
153 
154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155 
156 void HistoManager::EndOfRun()
157 {
158  if(fVerbose > 0) {
159  G4cout << "HistoManager: End of run actions are started" << G4endl;
160  }
161 
162  const G4Element* elm =
163  G4NistManager::Instance()->FindOrBuildElement(fElementName);
164  const G4Material* mat =
165  G4NistManager::Instance()->FindOrBuildMaterial("G4_"+fElementName);
166  const G4ParticleDefinition* particle =
167  G4ParticleTable::GetParticleTable()->FindParticle(fParticleName);
168 
169  G4cout << "### Fill Cross Sections for " << fParticleName
170  << " off " << fElementName
171  << G4endl;
172  if(fVerbose > 0) {
173  G4cout << "------------------------------------------------------------------------"
174  << G4endl;
175  G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
176  if(particle == fNeutron) { G4cout << " Capture(b) Fission(b)"; }
177  G4cout << " Total(b)" << G4endl;
178  G4cout << "------------------------------------------------------------------------"
179  << G4endl;
180  }
181  if(!particle || !elm) {
182  G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
183  return;
184  }
185 
186  G4int prec = G4cout.precision();
187  G4cout.precision(4);
188 
189  G4HadronicProcessStore* store = G4HadronicProcessStore::Instance();
190  G4double mass = particle->GetPDGMass();
191 
192  // Build histograms
193 
194  G4double p1 = std::log10(fMinMomentum/GeV);
195  G4double p2 = std::log10(fMaxMomentum/GeV);
196  G4double e1 = std::log10(fMinKinEnergy/MeV);
197  G4double e2 = std::log10(fMaxKinEnergy/MeV);
198  G4double de = (e2 - e1)/G4double(fBinsE);
199  G4double dp = (p2 - p1)/G4double(fBinsP);
200 
201  G4double x = e1 - de*0.5;
202  G4double e, p, xs, xtot;
203  G4int i;
204  for(i=0; i<fBinsE; i++) {
205  x += de;
206  e = std::pow(10.,x)*MeV;
207  if(fVerbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
208  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
209  xtot = xs;
210  if(fVerbose>0) G4cout << std::setw(12) << xs/barn;
211  fHisto->Fill(1, x, xs/barn);
212  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
213  xtot += xs;
214  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
215  fHisto->Fill(3, x, xs/barn);
216  if(particle == fNeutron) {
217  xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm,mat);
218  xtot += xs;
219  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
220  fHisto->Fill(4, x, xs/barn);
221  xs = store->GetFissionCrossSectionPerAtom(particle,e,elm,mat);
222  xtot += xs;
223  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
224  fHisto->Fill(5, x, xs/barn);
225  }
226  xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm,mat);
227  if(fVerbose>0) G4cout << " " << std::setw(12) << xtot/barn << G4endl;
228  fHisto->Fill(6, x, xs/barn);
229  fHisto->Fill(7, x, xtot/barn);
230  }
231 
232  x = p1 - dp*0.5;
233  for(i=0; i<fBinsP; i++) {
234  x += dp;
235  p = std::pow(10.,x)*GeV;
236  e = std::sqrt(p*p + mass*mass) - mass;
237  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
238  fHisto->Fill(0, x, xs/barn);
239  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
240  fHisto->Fill(2, x, xs/barn);
241  }
242  if(fVerbose > 0) {
243  G4cout << "-----------------------------------------------------------------"
244  << G4endl;
245  }
246  G4cout.precision(prec);
247  fHisto->Save();
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251 
252 void HistoManager::SetVerbose(G4int val)
253 {
254  fVerbose = val;
255  fHisto->SetVerbose(val);
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
259 
Generated on Sat May 25 2013 14:32:04 for Geant4 by   doxygen 1.8.4