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