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 75840 2013-11-06 17:28:38Z 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 //
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"
56 
57 #include "G4NucleiProperties.hh"
58 #include "G4NistManager.hh"
59 #include "G4StableIsotopes.hh"
60 #include "G4SystemOfUnits.hh"
61 
62 #include "HistoManagerMessenger.hh"
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65 
67 {
68  fAnalysisManager = 0;
69  fHistoName = "hadr00";
70 
71  fNeutron = G4Neutron::Neutron();
72  fMessenger = new HistoManagerMessenger(this);
73  fVerbose = 1;
74 
75  fParticleName = "proton";
76  fElementName = "Al";
77 
78  fMinKinEnergy = 0.1*MeV;
79  fMaxKinEnergy = 10*TeV;
80  fMinMomentum = 1*MeV;
81  fMaxMomentum = 10*TeV;
82 
83  fBinsE = 800;
84  fBinsP = 700;
85 }
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 
90 {
91  delete fMessenger;
92 }
93 
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95 
97 {
98  G4double p1 = std::log10(fMinMomentum/GeV);
99  G4double p2 = std::log10(fMaxMomentum/GeV);
100  G4double e1 = std::log10(fMinKinEnergy/MeV);
101  G4double e2 = std::log10(fMaxKinEnergy/MeV);
102 
103  //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
104 
105  fAnalysisManager = G4AnalysisManager::Instance();
106  fAnalysisManager->OpenFile(fHistoName+".root");
107  fAnalysisManager->SetFirstHistoId(1);
108 
109  fAnalysisManager->CreateH1("h1",
110  "Elastic cross section (barn) as a functions of log10(p/GeV)",
111  fBinsP,p1,p2);
112  fAnalysisManager->CreateH1("h2",
113  "Elastic cross section (barn) as a functions of log10(E/MeV)",
114  fBinsE,e1,e2);
115  fAnalysisManager->CreateH1("h3",
116  "Inelastic cross section (barn) as a functions of log10(p/GeV)",
117  fBinsP,p1,p2);
118  fAnalysisManager->CreateH1("h4",
119  "Inelastic cross section (barn) as a functions of log10(E/MeV)",
120  fBinsE,e1,e2);
121  fAnalysisManager->CreateH1("h5",
122  "Capture cross section (barn) as a functions of log10(E/MeV)",
123  fBinsE,e1,e2);
124  fAnalysisManager->CreateH1("h6",
125  "Fission cross section (barn) as a functions of log10(E/MeV)",
126  fBinsE,e1,e2);
127  fAnalysisManager->CreateH1("h7",
128  "Charge exchange cross section (barn) as a functions of log10(E/MeV)",
129  fBinsE,e1,e2);
130  fAnalysisManager->CreateH1("h8",
131  "Total cross section (barn) as a functions of log10(E/MeV)",
132  fBinsE,e1,e2);
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 
138 {
139  if(fVerbose > 0) {
140  G4cout << "HistoManager: End of run actions are started" << G4endl;
141  }
142 
143  const G4Element* elm =
145  const G4Material* mat =
146  G4NistManager::Instance()->FindOrBuildMaterial("G4_"+fElementName);
147  const G4ParticleDefinition* particle =
149 
150  G4cout << "### Fill Cross Sections for " << fParticleName
151  << " off " << fElementName
152  << G4endl;
153  if(fVerbose > 0) {
154  G4cout << "-------------------------------------------------------------"
155  << G4endl;
156  G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
157  if(particle == fNeutron) { G4cout << " Capture(b) Fission(b)"; }
158  G4cout << " Total(b)" << G4endl;
159  G4cout << "-------------------------------------------------------------"
160  << G4endl;
161  }
162  if(!particle || !elm) {
163  G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
164  return;
165  }
166 
167  G4int prec = G4cout.precision();
168  G4cout.precision(4);
169 
171  G4double mass = particle->GetPDGMass();
172 
173  // Build histograms
174 
175  G4double p1 = std::log10(fMinMomentum/GeV);
176  G4double p2 = std::log10(fMaxMomentum/GeV);
177  G4double e1 = std::log10(fMinKinEnergy/MeV);
178  G4double e2 = std::log10(fMaxKinEnergy/MeV);
179  G4double de = (e2 - e1)/G4double(fBinsE);
180  G4double dp = (p2 - p1)/G4double(fBinsP);
181 
182  G4double x = e1 - de*0.5;
183  G4double e, p, xs, xtot;
184  G4int i;
185  for(i=0; i<fBinsE; i++) {
186  x += de;
187  e = std::pow(10.,x)*MeV;
188  if(fVerbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
189  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
190  xtot = xs;
191  if(fVerbose>0) G4cout << std::setw(12) << xs/barn;
192  fAnalysisManager->FillH1(2, x, xs/barn);
193  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
194  xtot += xs;
195  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
196  fAnalysisManager->FillH1(4, x, xs/barn);
197  if(particle == fNeutron) {
198  xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm,mat);
199  xtot += xs;
200  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
201  fAnalysisManager->FillH1(5, x, xs/barn);
202  xs = store->GetFissionCrossSectionPerAtom(particle,e,elm,mat);
203  xtot += xs;
204  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
205  fAnalysisManager->FillH1(6, x, xs/barn);
206  }
207  xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm,mat);
208  if(fVerbose>0) G4cout << " " << std::setw(12) << xtot/barn << G4endl;
209  fAnalysisManager->FillH1(7, x, xs/barn);
210  fAnalysisManager->FillH1(8, x, xtot/barn);
211  }
212 
213  x = p1 - dp*0.5;
214  for(i=0; i<fBinsP; i++) {
215  x += dp;
216  p = std::pow(10.,x)*GeV;
217  e = std::sqrt(p*p + mass*mass) - mass;
218  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
219  fAnalysisManager->FillH1(1, x, xs/barn);
220  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
221  fAnalysisManager->FillH1(3, x, xs/barn);
222  }
223  if(fVerbose > 0) {
224  G4cout << "-------------------------------------------------------------"
225  << G4endl;
226  }
227  G4cout.precision(prec);
228  fAnalysisManager->Write();
229  fAnalysisManager->CloseFile();
230 
231  delete fAnalysisManager;
232 }
233 
234 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
235 
237 {
238  fVerbose = val;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static G4HadronicProcessStore * Instance()
G4double GetElasticCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)
const char * p
Definition: xmltok.h:285
tuple x
Definition: test.py:50
G4double GetFissionCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
Float_t mat
Definition: plot.C:40
G4GLOB_DLL std::ostream G4cout
G4double GetInelasticCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)
void BeginOfRun()
G4double GetChargeExchangeCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void SetVerbose(G4int value)
G4Element * FindOrBuildElement(G4int Z, G4bool isotopes=true)
G4double GetCaptureCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)