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