| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file hadronic/Hadr00/src/HistoManager.cc
27 /// \brief Implementation of the HistoManager
28 //
29 //
30 //--------------------------------------------
31 //
32 // ClassName: HistoManager
33 //
34 //
35 // Author: V.Ivanchenko 30/01/01
36 //
37 // Modified:
38 // 04.06.2006 Adoptation of hadr01 (V.Ivanchen
39 //
40 //--------------------------------------------
41 //
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 #include "HistoManager.hh"
47
48 #include "HistoManagerMessenger.hh"
49
50 #include "G4HadronicProcessStore.hh"
51 #include "G4Neutron.hh"
52 #include "G4NistManager.hh"
53 #include "G4NucleiProperties.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4ParticleTable.hh"
56 #include "G4StableIsotopes.hh"
57 #include "G4SystemOfUnits.hh"
58 #include "G4UnitsTable.hh"
59 #include "G4ios.hh"
60 #include "globals.hh"
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63
64 HistoManager::HistoManager()
65 {
66 fAnalysisManager = 0;
67 fHistoName = "hadr00";
68
69 fNeutron = G4Neutron::Neutron();
70 fMessenger = new HistoManagerMessenger(this)
71 fVerbose = 1;
72
73 fParticleName = "proton";
74 fElementName = "Al";
75
76 fTargetMaterial = 0;
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........oo
88
89 HistoManager::~HistoManager()
90 {
91 delete fMessenger;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oo
95
96 void HistoManager::BeginOfRun()
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= "<
104
105 fAnalysisManager = G4AnalysisManager::Instan
106 fAnalysisManager->OpenFile(fHistoName + ".ro
107 fAnalysisManager->SetFirstHistoId(1);
108
109 fAnalysisManager->CreateH1("h1", "Elastic cr
110 fBinsP, p1, p2);
111 fAnalysisManager->CreateH1("h2", "Elastic cr
112 fBinsE, e1, e2);
113 fAnalysisManager->CreateH1("h3", "Inelastic
114 fBinsP, p1, p2);
115 fAnalysisManager->CreateH1("h4", "Inelastic
116 fBinsE, e1, e2);
117 fAnalysisManager->CreateH1("h5", "Capture cr
118 fBinsE, e1, e2);
119 fAnalysisManager->CreateH1("h6", "Fission cr
120 fBinsE, e1, e2);
121 fAnalysisManager->CreateH1(
122 "h7", "Charge exchange cross section (barn
123 fAnalysisManager->CreateH1("h8", "Total cros
124 fBinsE, e1, e2);
125 fAnalysisManager->CreateH1(
126 "h9", "Inelastic cross section per volume
127 fAnalysisManager->CreateH1(
128 "h10", "Elastic cross section per volume a
129 }
130
131 //....oooOO0OOooo........oooOO0OOooo........oo
132
133 void HistoManager::EndOfRun()
134 {
135 if (fVerbose > 0) {
136 G4cout << "HistoManager: End of run action
137 }
138
139 const G4Element* elm = G4NistManager::Instan
140 const G4Material* mat = G4NistManager::Insta
141 const G4ParticleDefinition* particle =
142 G4ParticleTable::GetParticleTable()->FindP
143
144 G4cout << "### Fill Cross Sections for " <<
145 if (fVerbose > 0) {
146 G4cout << "-------------------------------
147 G4cout << " N E(MeV) Elastic(b)
148 if (particle == fNeutron) {
149 G4cout << " Capture(b) Fission(b)";
150 }
151 G4cout << " Total(b)" << G4endl;
152 G4cout << "-------------------------------
153 }
154 if (!particle || !elm) {
155 G4cout << "HistoManager WARNING Particle o
156 return;
157 }
158
159 G4int prec = G4cout.precision();
160 G4cout.precision(4);
161
162 G4HadronicProcessStore* store = G4HadronicPr
163 G4double mass = particle->GetPDGMass();
164
165 // Build histograms
166
167 G4double p1 = std::log10(fMinMomentum / GeV)
168 G4double p2 = std::log10(fMaxMomentum / GeV)
169 G4double e1 = std::log10(fMinKinEnergy / MeV
170 G4double e2 = std::log10(fMaxKinEnergy / MeV
171 G4double de = (e2 - e1) / G4double(fBinsE);
172 G4double dp = (p2 - p1) / G4double(fBinsP);
173
174 G4double x = e1 - de * 0.5;
175 G4double e, p, xs, xtot;
176 G4int i;
177
178 G4double coeff = 1.0;
179 if (fTargetMaterial) {
180 coeff = fTargetMaterial->GetDensity() * cm
181 }
182
183 for (i = 0; i < fBinsE; i++) {
184 x += de;
185 e = std::pow(10., x) * MeV;
186 if (fVerbose > 0) G4cout << std::setw(5) <
187 xs = store->GetElasticCrossSectionPerAtom(
188 xtot = xs;
189 if (fVerbose > 0) G4cout << std::setw(12)
190 fAnalysisManager->FillH1(2, x, xs / barn);
191 xs = store->GetInelasticCrossSectionPerAto
192 xtot += xs;
193 if (fVerbose > 0) G4cout << " " << std::se
194 fAnalysisManager->FillH1(4, x, xs / barn);
195 if (fTargetMaterial) {
196 xs = store->GetInelasticCrossSectionPerV
197 fAnalysisManager->FillH1(9, x, xs / coef
198 xs = store->GetElasticCrossSectionPerVol
199 fAnalysisManager->FillH1(10, x, xs / coe
200 }
201 if (particle == fNeutron) {
202 xs = store->GetCaptureCrossSectionPerAto
203 xtot += xs;
204 if (fVerbose > 0) G4cout << " " << std::
205 fAnalysisManager->FillH1(5, x, xs / barn
206 xs = store->GetFissionCrossSectionPerAto
207 xtot += xs;
208 if (fVerbose > 0) G4cout << " " << std::
209 fAnalysisManager->FillH1(6, x, xs / barn
210 }
211 xs = store->GetChargeExchangeCrossSectionP
212 if (fVerbose > 0) G4cout << " " << std::se
213 fAnalysisManager->FillH1(7, x, xs / barn);
214 fAnalysisManager->FillH1(8, x, xtot / barn
215 }
216
217 x = p1 - dp * 0.5;
218 for (i = 0; i < fBinsP; i++) {
219 x += dp;
220 p = std::pow(10., x) * GeV;
221 e = std::sqrt(p * p + mass * mass) - mass;
222 xs = store->GetElasticCrossSectionPerAtom(
223 fAnalysisManager->FillH1(1, x, xs / barn);
224 xs = store->GetInelasticCrossSectionPerAto
225 fAnalysisManager->FillH1(3, x, xs / barn);
226 }
227 if (fVerbose > 0) {
228 G4cout << "-------------------------------
229 }
230 G4cout.precision(prec);
231 fAnalysisManager->Write();
232 fAnalysisManager->CloseFile();
233 fAnalysisManager->Clear();
234
235 G4bool extra = true;
236 if (fTargetMaterial && extra) {
237 G4double E = 5 * GeV;
238 G4double cross = store->GetInelasticCrossS
239 if (cross <= 0.0) {
240 cross = 1.e-100;
241 }
242 G4cout << "### " << particle->GetParticleN
243 << fTargetMaterial->GetName()
244 << " xs/X0= " << 1.0 / (cross * fTa
245 }
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 void HistoManager::SetVerbose(G4int val)
251 {
252 fVerbose = val;
253 }
254
255 //....oooOO0OOooo........oooOO0OOooo........oo
256