| 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 // $Id: G4PenelopeCrossSection.cc 97613 2016-06-06 12:24:51Z gcosmo $
26 // 27 //
27 // Author: Luciano Pandola 28 // Author: Luciano Pandola
28 // 29 //
29 // History: 30 // History:
30 // -------- 31 // --------
31 // 18 Mar 2010 L Pandola First implementa 32 // 18 Mar 2010 L Pandola First implementation
32 // 09 Mar 2012 L. Pandola Add public metho 33 // 09 Mar 2012 L. Pandola Add public method (and machinery) to return
33 // the absolute and 34 // the absolute and the normalized shell cross
34 // sections indepen 35 // sections independently.
35 // 36 //
36 #include "G4PenelopeCrossSection.hh" 37 #include "G4PenelopeCrossSection.hh"
37 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh"
38 #include "G4PhysicsTable.hh" 39 #include "G4PhysicsTable.hh"
39 #include "G4PhysicsFreeVector.hh" 40 #include "G4PhysicsFreeVector.hh"
40 #include "G4DataVector.hh" 41 #include "G4DataVector.hh"
41 #include "G4Exp.hh" 42 #include "G4Exp.hh"
42 #include "G4Log.hh" <<
43 43
44 //....oooOO0OOooo........oooOO0OOooo........oo 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
45 G4PenelopeCrossSection::G4PenelopeCrossSection 45 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
46 fSoftCrossSections(nullptr), << 46 numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
47 fHardCrossSections(nullptr),fShellCrossSecti << 47 hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
48 fShellNormalizedCrossSections(nullptr), <<
49 fNumberOfEnergyPoints(nPointsE),fNumberOfShe <<
50 { 48 {
51 //check the number of points is not zero 49 //check the number of points is not zero
52 if (!fNumberOfEnergyPoints) << 50 if (!numberOfEnergyPoints)
53 { 51 {
54 G4ExceptionDescription ed; 52 G4ExceptionDescription ed;
55 ed << "G4PenelopeCrossSection: invalid n 53 ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
56 G4Exception("G4PenelopeCrossSection::G4P 54 G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
57 "em2017",FatalException,ed); 55 "em2017",FatalException,ed);
58 } 56 }
59 57
60 fIsNormalized = false; << 58 isNormalized = false;
61 59
62 // 1) soft XS table 60 // 1) soft XS table
63 fSoftCrossSections = new G4PhysicsTable(); << 61 softCrossSections = new G4PhysicsTable();
64 //the table contains 3 G4PhysicsFreeVectors, 62 //the table contains 3 G4PhysicsFreeVectors,
65 //(fSoftCrossSections)[0] --> log XS0 vs. l << 63 //(softCrossSections)[0] --> log XS0 vs. log E
66 //(fSoftCrossSections)[1] --> log XS1 vs. l << 64 //(softCrossSections)[1] --> log XS1 vs. log E
67 //(fSoftCrossSections)[2] --> log XS2 vs. l << 65 //(softCrossSections)[2] --> log XS2 vs. log E
68 66
69 //everything is log-log 67 //everything is log-log
70 for (size_t i=0;i<3;i++) 68 for (size_t i=0;i<3;i++)
71 fSoftCrossSections->push_back(new G4Physic << 69 softCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
72 70
73 //2) hard XS table 71 //2) hard XS table
74 fHardCrossSections = new G4PhysicsTable(); << 72 hardCrossSections = new G4PhysicsTable();
75 //the table contains 3 G4PhysicsFreeVectors, 73 //the table contains 3 G4PhysicsFreeVectors,
76 //(fHardCrossSections)[0] --> log XH0 vs. l << 74 //(hardCrossSections)[0] --> log XH0 vs. log E
77 //(fHardCrossSections)[1] --> log XH1 vs. l << 75 //(hardCrossSections)[1] --> log XH1 vs. log E
78 //(fHardCrossSections)[2] --> log XH2 vs. l << 76 //(hardCrossSections)[2] --> log XH2 vs. log E
79 77
80 //everything is log-log 78 //everything is log-log
81 for (size_t i=0;i<3;i++) 79 for (size_t i=0;i<3;i++)
82 fHardCrossSections->push_back(new G4Physic << 80 hardCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
83 81
84 //3) shell XS table, if it is the case 82 //3) shell XS table, if it is the case
85 if (fNumberOfShells) << 83 if (numberOfShells)
86 { 84 {
87 fShellCrossSections = new G4PhysicsTable << 85 shellCrossSections = new G4PhysicsTable();
88 fShellNormalizedCrossSections = new G4Ph << 86 shellNormalizedCrossSections = new G4PhysicsTable();
89 //the table has to contain numberofShell 87 //the table has to contain numberofShells G4PhysicsFreeVectors,
90 //(theTable)[ishell] --> cross section f 88 //(theTable)[ishell] --> cross section for shell #ishell
91 for (size_t i=0;i<fNumberOfShells;i++) << 89 for (size_t i=0;i<numberOfShells;i++)
92 { 90 {
93 fShellCrossSections->push_back(new G4Physi << 91 shellCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
94 fShellNormalizedCrossSections->push_back(n << 92 shellNormalizedCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
95 } 93 }
96 } 94 }
97 } 95 }
98 96
99 //....oooOO0OOooo........oooOO0OOooo........oo 97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
100 G4PenelopeCrossSection::~G4PenelopeCrossSectio 98 G4PenelopeCrossSection::~G4PenelopeCrossSection()
101 { 99 {
102 //clean up tables 100 //clean up tables
103 if (fShellCrossSections) << 101 if (shellCrossSections)
104 { 102 {
105 fShellCrossSections->clearAndDestroy(); << 103 //shellCrossSections->clearAndDestroy();
106 delete fShellCrossSections; << 104 delete shellCrossSections;
107 } 105 }
108 if (fShellNormalizedCrossSections) << 106 if (shellNormalizedCrossSections)
109 { 107 {
110 fShellNormalizedCrossSections->clearAndD << 108 //shellNormalizedCrossSections->clearAndDestroy();
111 delete fShellNormalizedCrossSections; << 109 delete shellNormalizedCrossSections;
112 } 110 }
113 if (fSoftCrossSections) << 111 if (softCrossSections)
114 { 112 {
115 fSoftCrossSections->clearAndDestroy(); << 113 //softCrossSections->clearAndDestroy();
116 delete fSoftCrossSections; << 114 delete softCrossSections;
117 } 115 }
118 if (fHardCrossSections) << 116 if (hardCrossSections)
119 { 117 {
120 fHardCrossSections->clearAndDestroy(); << 118 //hardCrossSections->clearAndDestroy();
121 delete fHardCrossSections; << 119 delete hardCrossSections;
122 } 120 }
123 } 121 }
124 122
125 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
126 void G4PenelopeCrossSection::AddCrossSectionPo 124 void G4PenelopeCrossSection::AddCrossSectionPoint(size_t binNumber,G4double energy,
127 G4double XH0, 125 G4double XH0,
128 G4double XH1, G4double XH2, 126 G4double XH1, G4double XH2,
129 G4double XS0, G4double XS1, 127 G4double XS0, G4double XS1,
130 G4double XS2) 128 G4double XS2)
131 { 129 {
132 if (!fSoftCrossSections || !fHardCrossSectio << 130 if (!softCrossSections || !hardCrossSections)
133 { 131 {
134 G4cout << "Something wrong in G4Penelope 132 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
135 G4endl; 133 G4endl;
136 G4cout << "Trying to fill un-initialized 134 G4cout << "Trying to fill un-initialized tables" << G4endl;
137 return; 135 return;
138 } 136 }
139 137
140 //fill vectors 138 //fill vectors
141 G4PhysicsFreeVector* theVector = (G4PhysicsF << 139 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
142 140
143 if (binNumber >= fNumberOfEnergyPoints) << 141 if (binNumber >= numberOfEnergyPoints)
144 { 142 {
145 G4cout << "Something wrong in G4Penelope 143 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
146 G4endl; 144 G4endl;
147 G4cout << "Trying to register more point 145 G4cout << "Trying to register more points than originally declared" << G4endl;
148 return; 146 return;
149 } 147 }
150 G4double logEne = G4Log(energy); << 148 G4double logEne = std::log(energy);
151 149
152 //XS0 150 //XS0
153 G4double val = G4Log(std::max(XS0,1e-42*cm2 << 151 G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
154 theVector->PutValues(binNumber,logEne,val); << 152 theVector->PutValue(binNumber,logEne,val);
155 153
156 //XS1 154 //XS1
157 theVector = (G4PhysicsFreeVector*) (*fSoftC << 155 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
158 val = G4Log(std::max(XS1,1e-42*eV*cm2)); / << 156 val = std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
159 theVector->PutValues(binNumber,logEne,val); << 157 theVector->PutValue(binNumber,logEne,val);
160 158
161 //XS2 159 //XS2
162 theVector = (G4PhysicsFreeVector*) (*fSoftC << 160 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
163 val = G4Log(std::max(XS2,1e-42*eV*eV*cm2)) << 161 val = std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
164 theVector->PutValues(binNumber,logEne,val); << 162 theVector->PutValue(binNumber,logEne,val);
165 163
166 //XH0 164 //XH0
167 theVector = (G4PhysicsFreeVector*) (*fHardC << 165 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
168 val = G4Log(std::max(XH0,1e-42*cm2)); //av << 166 val = std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
169 theVector->PutValues(binNumber,logEne,val); << 167 theVector->PutValue(binNumber,logEne,val);
170 168
171 //XH1 169 //XH1
172 theVector = (G4PhysicsFreeVector*) (*fHardC << 170 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
173 val = G4Log(std::max(XH1,1e-42*eV*cm2)); / << 171 val = std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
174 theVector->PutValues(binNumber,logEne,val); << 172 theVector->PutValue(binNumber,logEne,val);
175 173
176 //XH2 174 //XH2
177 theVector = (G4PhysicsFreeVector*) (*fHardC << 175 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
178 val = G4Log(std::max(XH2,1e-42*eV*eV*cm2)) << 176 val = std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
179 theVector->PutValues(binNumber,logEne,val); << 177 theVector->PutValue(binNumber,logEne,val);
180 178
181 return; 179 return;
182 } 180 }
183 181
184 //....oooOO0OOooo........oooOO0OOooo........oo 182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
185 183
186 void G4PenelopeCrossSection::AddShellCrossSect 184 void G4PenelopeCrossSection::AddShellCrossSectionPoint(size_t binNumber,
187 size_t shellID, 185 size_t shellID,
188 G4double energy, 186 G4double energy,
189 G4double xs) 187 G4double xs)
190 { 188 {
191 if (!fShellCrossSections) << 189 if (!shellCrossSections)
192 { 190 {
193 G4cout << "Something wrong in G4Penelope 191 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
194 G4endl; 192 G4endl;
195 G4cout << "Trying to fill un-initialized 193 G4cout << "Trying to fill un-initialized table" << G4endl;
196 return; 194 return;
197 } 195 }
198 196
199 if (shellID >= fNumberOfShells) << 197 if (shellID >= numberOfShells)
200 { 198 {
201 G4cout << "Something wrong in G4Penelope 199 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
202 G4endl; 200 G4endl;
203 G4cout << "Trying to fill shell #" << sh 201 G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
204 << fNumberOfShells-1 << G4endl; << 202 << numberOfShells-1 << G4endl;
205 return; 203 return;
206 } 204 }
207 205
208 //fill vector 206 //fill vector
209 G4PhysicsFreeVector* theVector = (G4PhysicsF << 207 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
210 208
211 if (binNumber >= fNumberOfEnergyPoints) << 209 if (binNumber >= numberOfEnergyPoints)
212 { 210 {
213 G4cout << "Something wrong in G4Penelope 211 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
214 G4endl; 212 G4endl;
215 G4cout << "Trying to register more point 213 G4cout << "Trying to register more points than originally declared" << G4endl;
216 return; 214 return;
217 } 215 }
218 G4double logEne = G4Log(energy); << 216 G4double logEne = std::log(energy);
219 G4double val = G4Log(std::max(xs,1e-42*cm2) << 217 G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
220 theVector->PutValues(binNumber,logEne,val); << 218 theVector->PutValue(binNumber,logEne,val);
221 219
222 return; 220 return;
223 } 221 }
224 222
225 //....oooOO0OOooo........oooOO0OOooo........oo 223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
226 224
227 G4double G4PenelopeCrossSection::GetTotalCross 225 G4double G4PenelopeCrossSection::GetTotalCrossSection(G4double energy) const
228 { 226 {
229 G4double result = 0; 227 G4double result = 0;
230 //take here XS0 + XH0 228 //take here XS0 + XH0
231 if (!fSoftCrossSections || !fHardCrossSectio << 229 if (!softCrossSections || !hardCrossSections)
232 { 230 {
233 G4cout << "Something wrong in G4Penelope 231 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
234 G4endl; 232 G4endl;
235 G4cout << "Trying to retrieve from un-in 233 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
236 return result; 234 return result;
237 } 235 }
238 236
239 // 1) soft part 237 // 1) soft part
240 G4PhysicsFreeVector* theVector = (G4PhysicsF << 238 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
241 if (theVector->GetVectorLength() < fNumberOf << 239 if (theVector->GetVectorLength() < numberOfEnergyPoints)
242 { 240 {
243 G4cout << "Something wrong in G4Penelope 241 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
244 G4endl; 242 G4endl;
245 G4cout << "Soft cross section table look 243 G4cout << "Soft cross section table looks not filled" << G4endl;
246 return result; 244 return result;
247 } 245 }
248 G4double logene = G4Log(energy); << 246 G4double logene = std::log(energy);
249 G4double logXS = theVector->Value(logene); 247 G4double logXS = theVector->Value(logene);
250 G4double softXS = G4Exp(logXS); 248 G4double softXS = G4Exp(logXS);
251 249
252 // 2) hard part 250 // 2) hard part
253 theVector = (G4PhysicsFreeVector*) (*fHardCr << 251 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
254 if (theVector->GetVectorLength() < fNumberOf << 252 if (theVector->GetVectorLength() < numberOfEnergyPoints)
255 { 253 {
256 G4cout << "Something wrong in G4Penelope 254 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
257 G4endl; 255 G4endl;
258 G4cout << "Hard cross section table look 256 G4cout << "Hard cross section table looks not filled" << G4endl;
259 return result; 257 return result;
260 } 258 }
261 logXS = theVector->Value(logene); 259 logXS = theVector->Value(logene);
262 G4double hardXS = G4Exp(logXS); 260 G4double hardXS = G4Exp(logXS);
263 261
264 result = hardXS + softXS; 262 result = hardXS + softXS;
265 return result; 263 return result;
>> 264
266 } 265 }
267 266
268 //....oooOO0OOooo........oooOO0OOooo........oo 267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
269 268
270 G4double G4PenelopeCrossSection::GetHardCrossS 269 G4double G4PenelopeCrossSection::GetHardCrossSection(G4double energy) const
271 { 270 {
272 G4double result = 0; 271 G4double result = 0;
273 //take here XH0 272 //take here XH0
274 if (!fHardCrossSections) << 273 if (!hardCrossSections)
275 { 274 {
276 G4cout << "Something wrong in G4Penelope 275 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
277 G4endl; 276 G4endl;
278 G4cout << "Trying to retrieve from un-in 277 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
279 return result; 278 return result;
280 } 279 }
281 280
282 G4PhysicsFreeVector* theVector = (G4PhysicsF << 281 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
283 if (theVector->GetVectorLength() < fNumberOf << 282 if (theVector->GetVectorLength() < numberOfEnergyPoints)
284 { 283 {
285 G4cout << "Something wrong in G4Penelope 284 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
286 G4endl; 285 G4endl;
287 G4cout << "Hard cross section table look 286 G4cout << "Hard cross section table looks not filled" << G4endl;
288 return result; 287 return result;
289 } 288 }
290 G4double logene = G4Log(energy); << 289 G4double logene = std::log(energy);
291 G4double logXS = theVector->Value(logene); 290 G4double logXS = theVector->Value(logene);
292 result = G4Exp(logXS); 291 result = G4Exp(logXS);
293 292
294 return result; 293 return result;
295 } 294 }
296 295
>> 296
297 //....oooOO0OOooo........oooOO0OOooo........oo 297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
298 298
299 G4double G4PenelopeCrossSection::GetSoftStoppi 299 G4double G4PenelopeCrossSection::GetSoftStoppingPower(G4double energy) const
300 { 300 {
301 G4double result = 0; 301 G4double result = 0;
302 //take here XH0 302 //take here XH0
303 if (!fSoftCrossSections) << 303 if (!softCrossSections)
304 { 304 {
305 G4cout << "Something wrong in G4Penelope 305 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
306 G4endl; 306 G4endl;
307 G4cout << "Trying to retrieve from un-in 307 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
308 return result; 308 return result;
309 } 309 }
310 310
311 G4PhysicsFreeVector* theVector = (G4PhysicsF << 311 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
312 if (theVector->GetVectorLength() < fNumberOf << 312 if (theVector->GetVectorLength() < numberOfEnergyPoints)
313 { 313 {
314 G4cout << "Something wrong in G4Penelope 314 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
315 G4endl; 315 G4endl;
316 G4cout << "Soft cross section table look 316 G4cout << "Soft cross section table looks not filled" << G4endl;
317 return result; 317 return result;
318 } 318 }
319 G4double logene = G4Log(energy); << 319 G4double logene = std::log(energy);
320 G4double logXS = theVector->Value(logene); 320 G4double logXS = theVector->Value(logene);
321 result = G4Exp(logXS); 321 result = G4Exp(logXS);
322 322
323 return result; 323 return result;
324 } 324 }
325 325
326 //....oooOO0OOooo........oooOO0OOooo........oo 326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
327 327
328 G4double G4PenelopeCrossSection::GetShellCross 328 G4double G4PenelopeCrossSection::GetShellCrossSection(size_t shellID,G4double energy) const
329 { 329 {
330 G4double result = 0; 330 G4double result = 0;
331 if (!fShellCrossSections) << 331 if (!shellCrossSections)
332 { 332 {
333 G4cout << "Something wrong in G4Penelope 333 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
334 G4endl; 334 G4endl;
335 G4cout << "Trying to retrieve from un-in 335 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
336 return result; 336 return result;
337 } 337 }
338 if (shellID >= fNumberOfShells) << 338 if (shellID >= numberOfShells)
339 { 339 {
340 G4cout << "Something wrong in G4Penelope 340 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
341 G4endl; 341 G4endl;
342 G4cout << "Trying to retrieve shell #" < 342 G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
343 << fNumberOfShells-1 << G4endl; << 343 << numberOfShells-1 << G4endl;
344 return result; 344 return result;
345 } 345 }
346 346
347 G4PhysicsFreeVector* theVector = (G4PhysicsF << 347 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
348 348
349 if (theVector->GetVectorLength() < fNumberOf << 349 if (theVector->GetVectorLength() < numberOfEnergyPoints)
350 { 350 {
351 G4cout << "Something wrong in G4Penelope 351 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
352 G4endl; 352 G4endl;
353 G4cout << "Shell cross section table loo 353 G4cout << "Shell cross section table looks not filled" << G4endl;
354 return result; 354 return result;
355 } 355 }
356 G4double logene = G4Log(energy); << 356 G4double logene = std::log(energy);
357 G4double logXS = theVector->Value(logene); 357 G4double logXS = theVector->Value(logene);
358 result = G4Exp(logXS); 358 result = G4Exp(logXS);
359 359
360 return result; 360 return result;
361 } 361 }
362 //....oooOO0OOooo........oooOO0OOooo........oo 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
363 363
364 G4double G4PenelopeCrossSection::GetNormalized 364 G4double G4PenelopeCrossSection::GetNormalizedShellCrossSection(size_t shellID,G4double energy) const
365 { 365 {
366 G4double result = 0; 366 G4double result = 0;
367 if (!fShellNormalizedCrossSections) << 367 if (!shellNormalizedCrossSections)
368 { 368 {
369 G4cout << "Something wrong in G4Penelope 369 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
370 G4endl; 370 G4endl;
371 G4cout << "Trying to retrieve from un-in 371 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
372 return result; 372 return result;
373 } 373 }
374 374
375 if (!fIsNormalized) << 375 if (!isNormalized)
376 { 376 {
377 G4cout << "Something wrong in G4Penelope 377 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
378 G4cout << "The table of normalized cross 378 G4cout << "The table of normalized cross section is not initialized" << G4endl;
379 } 379 }
380 380
381 if (shellID >= fNumberOfShells) << 381
>> 382 if (shellID >= numberOfShells)
382 { 383 {
383 G4cout << "Something wrong in G4Penelope 384 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
384 G4endl; 385 G4endl;
385 G4cout << "Trying to retrieve shell #" < 386 G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
386 << fNumberOfShells-1 << G4endl; << 387 << numberOfShells-1 << G4endl;
387 return result; 388 return result;
388 } 389 }
389 390
390 const G4PhysicsFreeVector* theVector = 391 const G4PhysicsFreeVector* theVector =
391 (G4PhysicsFreeVector*) (*fShellNormalizedC << 392 (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
392 393
393 if (theVector->GetVectorLength() < fNumberOf << 394 if (theVector->GetVectorLength() < numberOfEnergyPoints)
394 { 395 {
395 G4cout << "Something wrong in G4Penelope 396 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
396 G4endl; 397 G4endl;
397 G4cout << "Shell cross section table loo 398 G4cout << "Shell cross section table looks not filled" << G4endl;
398 return result; 399 return result;
399 } 400 }
400 G4double logene = G4Log(energy); << 401 G4double logene = std::log(energy);
401 G4double logXS = theVector->Value(logene); 402 G4double logXS = theVector->Value(logene);
402 result = G4Exp(logXS); 403 result = G4Exp(logXS);
403 404
404 return result; 405 return result;
405 } 406 }
406 407
>> 408
>> 409
407 //....oooOO0OOooo........oooOO0OOooo........oo 410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
408 411
409 void G4PenelopeCrossSection::NormalizeShellCro 412 void G4PenelopeCrossSection::NormalizeShellCrossSections()
410 { 413 {
411 if (fIsNormalized) //already done! << 414 if (isNormalized) //already done!
412 { 415 {
413 G4cout << "G4PenelopeCrossSection::Norma 416 G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
414 G4cout << "already invoked. Ignore it" < 417 G4cout << "already invoked. Ignore it" << G4endl;
415 return; 418 return;
416 } 419 }
417 420
418 if (!fShellNormalizedCrossSections) << 421 if (!shellNormalizedCrossSections)
419 { 422 {
420 G4cout << "Something wrong in G4Penelope 423 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
421 G4endl; 424 G4endl;
422 G4cout << "Trying to retrieve from un-in 425 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
423 return; 426 return;
424 } 427 }
425 428
426 for (size_t i=0;i<fNumberOfEnergyPoints;i++) << 429 for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
427 { 430 {
428 //energy grid is the same for all shells 431 //energy grid is the same for all shells
429 432
430 //Recalculate manually the XS factor, to 433 //Recalculate manually the XS factor, to avoid problems with
431 //underflows 434 //underflows
432 G4double normFactor = 0.; 435 G4double normFactor = 0.;
433 for (size_t shellID=0;shellID<fNumberOfS << 436 for (size_t shellID=0;shellID<numberOfShells;shellID++)
434 { 437 {
435 G4PhysicsFreeVector* theVec = 438 G4PhysicsFreeVector* theVec =
436 (G4PhysicsFreeVector*) (*fShellCrossSect << 439 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
437 440
438 normFactor += G4Exp((*theVec)[i]); 441 normFactor += G4Exp((*theVec)[i]);
439 } 442 }
440 G4double logNormFactor = G4Log(normFacto << 443 G4double logNormFactor = std::log(normFactor);
441 //Normalize 444 //Normalize
442 for (size_t shellID=0;shellID<fNumberOfS << 445 for (size_t shellID=0;shellID<numberOfShells;shellID++)
443 { 446 {
444 G4PhysicsFreeVector* theVec = 447 G4PhysicsFreeVector* theVec =
445 (G4PhysicsFreeVector*) (*fShellNormalize << 448 (G4PhysicsFreeVector*) (*shellNormalizedCrossSections)[shellID];
446 G4PhysicsFreeVector* theFullVec = 449 G4PhysicsFreeVector* theFullVec =
447 (G4PhysicsFreeVector*) (*fShellCrossSecti << 450 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
448 G4double previousValue = (*theFullVec)[i]; 451 G4double previousValue = (*theFullVec)[i]; //log(XS)
449 G4double logEnergy = theFullVec->GetLowEdge 452 G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
450 //log(XS/normFactor) = log(XS) - log(normFa 453 //log(XS/normFactor) = log(XS) - log(normFactor)
451 theVec->PutValues(i,logEnergy,previousValue << 454 theVec->PutValue(i,logEnergy,previousValue-logNormFactor);
>> 455 }
>> 456 }
>> 457
>> 458 isNormalized = true;
>> 459
>> 460
>> 461 /*
>> 462 //TESTING
>> 463 for (size_t shellID=0;shellID<numberOfShells;shellID++)
>> 464 {
>> 465 G4cout << "SHELL " << shellID << G4endl;
>> 466 G4PhysicsFreeVector* theVec =
>> 467 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
>> 468 for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
>> 469 {
>> 470 G4double logene = theVec->GetLowEdgeEnergy(i);
>> 471 G4cout << G4Exp(logene)/MeV << " " << G4Exp((*theVec)[i]) << G4endl;
452 } 472 }
453 } 473 }
454 fIsNormalized = true; << 474 */
455 475
456 return; 476 return;
457 } 477 }
458 478