Geant4  10.03.p03
 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 96284 2016-04-04 07:19:26Z 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  fTargetMaterial = 0;
79 
80  fMinKinEnergy = 0.1*MeV;
81  fMaxKinEnergy = 10*TeV;
82  fMinMomentum = 1*MeV;
83  fMaxMomentum = 10*TeV;
84 
85  fBinsE = 800;
86  fBinsP = 700;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
92 {
93  delete fMessenger;
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
97 
99 {
100  G4double p1 = std::log10(fMinMomentum/GeV);
101  G4double p2 = std::log10(fMaxMomentum/GeV);
102  G4double e1 = std::log10(fMinKinEnergy/MeV);
103  G4double e2 = std::log10(fMaxKinEnergy/MeV);
104 
105  //G4cout<<"e1= "<<e1<<" e2= "<<e2<<" p1= "<<p1<<" p2= "<<p2<<G4endl;
106 
107  fAnalysisManager = G4AnalysisManager::Instance();
108  fAnalysisManager->OpenFile(fHistoName+".root");
109  fAnalysisManager->SetFirstHistoId(1);
110 
111  fAnalysisManager->CreateH1("h1",
112  "Elastic cross section (barn) as a functions of log10(p/GeV)",
113  fBinsP,p1,p2);
114  fAnalysisManager->CreateH1("h2",
115  "Elastic cross section (barn) as a functions of log10(E/MeV)",
116  fBinsE,e1,e2);
117  fAnalysisManager->CreateH1("h3",
118  "Inelastic cross section (barn) as a functions of log10(p/GeV)",
119  fBinsP,p1,p2);
120  fAnalysisManager->CreateH1("h4",
121  "Inelastic cross section (barn) as a functions of log10(E/MeV)",
122  fBinsE,e1,e2);
123  fAnalysisManager->CreateH1("h5",
124  "Capture cross section (barn) as a functions of log10(E/MeV)",
125  fBinsE,e1,e2);
126  fAnalysisManager->CreateH1("h6",
127  "Fission cross section (barn) as a functions of log10(E/MeV)",
128  fBinsE,e1,e2);
129  fAnalysisManager->CreateH1("h7",
130  "Charge exchange cross section (barn) as a functions of log10(E/MeV)",
131  fBinsE,e1,e2);
132  fAnalysisManager->CreateH1("h8",
133  "Total cross section (barn) as a functions of log10(E/MeV)",
134  fBinsE,e1,e2);
135  fAnalysisManager->CreateH1("h9",
136  "Inelastic cross section per volume as a functions of log10(E/MeV)",
137  fBinsE,e1,e2);
138  fAnalysisManager->CreateH1("h10",
139  "Elastic cross section per volume as a functions of log10(E/MeV)",
140  fBinsE,e1,e2);
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
146 {
147  if(fVerbose > 0) {
148  G4cout << "HistoManager: End of run actions are started" << G4endl;
149  }
150 
151  const G4Element* elm =
153  const G4Material* mat =
154  G4NistManager::Instance()->FindOrBuildMaterial("G4_"+fElementName);
155  const G4ParticleDefinition* particle =
157 
158  G4cout << "### Fill Cross Sections for " << fParticleName
159  << " off " << fElementName
160  << G4endl;
161  if(fVerbose > 0) {
162  G4cout << "-------------------------------------------------------------"
163  << G4endl;
164  G4cout << " N E(MeV) Elastic(b) Inelastic(b)";
165  if(particle == fNeutron) { G4cout << " Capture(b) Fission(b)"; }
166  G4cout << " Total(b)" << G4endl;
167  G4cout << "-------------------------------------------------------------"
168  << G4endl;
169  }
170  if(!particle || !elm) {
171  G4cout << "HistoManager WARNING Particle or element undefined" << G4endl;
172  return;
173  }
174 
175  G4int prec = G4cout.precision();
176  G4cout.precision(4);
177 
179  G4double mass = particle->GetPDGMass();
180 
181  // Build histograms
182 
183  G4double p1 = std::log10(fMinMomentum/GeV);
184  G4double p2 = std::log10(fMaxMomentum/GeV);
185  G4double e1 = std::log10(fMinKinEnergy/MeV);
186  G4double e2 = std::log10(fMaxKinEnergy/MeV);
187  G4double de = (e2 - e1)/G4double(fBinsE);
188  G4double dp = (p2 - p1)/G4double(fBinsP);
189 
190  G4double x = e1 - de*0.5;
191  G4double e, p, xs, xtot;
192  G4int i;
193 
194  G4double coeff = 1.0;
195  if(fTargetMaterial) { coeff = fTargetMaterial->GetDensity()*cm2/g; }
196 
197  for(i=0; i<fBinsE; i++) {
198  x += de;
199  e = std::pow(10.,x)*MeV;
200  if(fVerbose>0) G4cout << std::setw(5) << i << std::setw(12) << e;
201  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
202  xtot = xs;
203  if(fVerbose>0) G4cout << std::setw(12) << xs/barn;
204  fAnalysisManager->FillH1(2, x, xs/barn);
205  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
206  xtot += xs;
207  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
208  fAnalysisManager->FillH1(4, x, xs/barn);
209  if(fTargetMaterial) {
210  xs =
211  store->GetInelasticCrossSectionPerVolume(particle,e,fTargetMaterial);
212  fAnalysisManager->FillH1(9, x, xs/coeff);
213  xs =
214  store->GetElasticCrossSectionPerVolume(particle,e,fTargetMaterial);
215  fAnalysisManager->FillH1(10, x, xs/coeff);
216  }
217  if(particle == fNeutron) {
218  xs = store->GetCaptureCrossSectionPerAtom(particle,e,elm,mat);
219  xtot += xs;
220  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
221  fAnalysisManager->FillH1(5, x, xs/barn);
222  xs = store->GetFissionCrossSectionPerAtom(particle,e,elm,mat);
223  xtot += xs;
224  if(fVerbose>0) G4cout << " " << std::setw(12) << xs/barn;
225  fAnalysisManager->FillH1(6, x, xs/barn);
226  }
227  xs = store->GetChargeExchangeCrossSectionPerAtom(particle,e,elm,mat);
228  if(fVerbose>0) G4cout << " " << std::setw(12) << xtot/barn << G4endl;
229  fAnalysisManager->FillH1(7, x, xs/barn);
230  fAnalysisManager->FillH1(8, x, xtot/barn);
231  }
232 
233  x = p1 - dp*0.5;
234  for(i=0; i<fBinsP; i++) {
235  x += dp;
236  p = std::pow(10.,x)*GeV;
237  e = std::sqrt(p*p + mass*mass) - mass;
238  xs = store->GetElasticCrossSectionPerAtom(particle,e,elm,mat);
239  fAnalysisManager->FillH1(1, x, xs/barn);
240  xs = store->GetInelasticCrossSectionPerAtom(particle,e,elm,mat);
241  fAnalysisManager->FillH1(3, x, xs/barn);
242  }
243  if(fVerbose > 0) {
244  G4cout << "-------------------------------------------------------------"
245  << G4endl;
246  }
247  G4cout.precision(prec);
248  fAnalysisManager->Write();
249  fAnalysisManager->CloseFile();
250 
251  G4bool extra = true;
252  if(fTargetMaterial && extra) {
253  G4double E= 5*GeV;
254  G4double cross =
255  store->GetInelasticCrossSectionPerVolume(particle,E,fTargetMaterial);
256  if(cross <= 0.0) { cross = 1.e-100; }
257  G4cout << "### " << particle->GetParticleName() << " " << E/GeV
258  << " GeV on " << fTargetMaterial->GetName()
259  << " xs/X0= " << 1.0/(cross*fTargetMaterial->GetRadlen()) << G4endl;
260  }
261  delete fAnalysisManager;
262 }
263 
264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
265 
267 {
268  fVerbose = val;
269 }
270 
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
272 
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
static constexpr double cm2
Definition: G4SIunits.hh:120
static G4HadronicProcessStore * Instance()
G4double GetElasticCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetDensity() const
Definition: G4Material.hh:180
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()
const G4String & GetParticleName() const
static constexpr double TeV
Definition: G4SIunits.hh:218
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
static const double prec
Definition: RanecuEngine.cc:58
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
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 GetRadlen() const
Definition: G4Material.hh:220
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetVerbose(G4int val)
Definition: HistoManager.hh:95
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double GetElasticCrossSectionPerVolume(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Material *material)
double G4double
Definition: G4Types.hh:76
G4double GetInelasticCrossSectionPerVolume(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Material *material)
static constexpr double barn
Definition: G4SIunits.hh:105
G4Element * FindOrBuildElement(G4int Z, G4bool isotopes=true)
G4double GetCaptureCrossSectionPerAtom(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4Element *anElement, const G4Material *mat=0)