source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeCrossSection.cc

Last change on this file was 1350, checked in by garnier, 15 years ago

update to last version 4.9.4
File size: 13.2 KB
Line 
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//
26// $Id: G4PenelopeCrossSection.cc,v 1.2 2010/12/15 07:39:14 gunter Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 18 Mar 2010 L Pandola First implementation
34//
35#include "G4PenelopeCrossSection.hh"
36#include "G4PhysicsTable.hh"
37#include "G4PhysicsFreeVector.hh"
38#include "G4DataVector.hh"
39
40//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
41G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
42 numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
43 hardCrossSections(0),shellCrossSections(0)
44{
45 //check the number of points is not zero
46 if (!numberOfEnergyPoints)
47 {
48 G4cout << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
49 G4Exception();
50 }
51
52 isNormalized = false;
53
54 // 1) soft XS table
55 softCrossSections = new G4PhysicsTable();
56 //the table contains 3 G4PhysicsFreeVectors,
57 //(softCrossSections)[0] --> log XS0 vs. log E
58 //(softCrossSections)[1] --> log XS1 vs. log E
59 //(softCrossSections)[2] --> log XS2 vs. log E
60
61 //everything is log-log
62 for (size_t i=0;i<3;i++)
63 softCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
64
65 //2) hard XS table
66 hardCrossSections = new G4PhysicsTable();
67 //the table contains 3 G4PhysicsFreeVectors,
68 //(hardCrossSections)[0] --> log XH0 vs. log E
69 //(hardCrossSections)[1] --> log XH1 vs. log E
70 //(hardCrossSections)[2] --> log XH2 vs. log E
71
72 //everything is log-log
73 for (size_t i=0;i<3;i++)
74 hardCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
75
76 //3) shell XS table, if it is the case
77 if (numberOfShells)
78 {
79 shellCrossSections = new G4PhysicsTable();
80 //the table has to contain numberofShells G4PhysicsFreeVectors,
81 //(theTable)[ishell] --> cross section for shell #ishell
82 for (size_t i=0;i<numberOfShells;i++)
83 shellCrossSections->push_back(new G4PhysicsFreeVector(numberOfEnergyPoints));
84 }
85}
86
87//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
88G4PenelopeCrossSection::~G4PenelopeCrossSection()
89{
90 //clean up tables
91 if (shellCrossSections)
92 {
93 shellCrossSections->clearAndDestroy();
94 delete shellCrossSections;
95 }
96 if (softCrossSections)
97 {
98 softCrossSections->clearAndDestroy();
99 delete softCrossSections;
100 }
101 if (hardCrossSections)
102 {
103 hardCrossSections->clearAndDestroy();
104 delete hardCrossSections;
105 }
106}
107
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
109void G4PenelopeCrossSection::AddCrossSectionPoint(size_t binNumber,G4double energy,
110 G4double XH0,
111 G4double XH1, G4double XH2,
112 G4double XS0, G4double XS1,
113 G4double XS2)
114{
115 if (!softCrossSections || !hardCrossSections)
116 {
117 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
118 G4endl;
119 G4cout << "Trying to fill un-initialized tables" << G4endl;
120 return;
121 }
122
123 //fill vectors
124 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
125
126 if (binNumber >= numberOfEnergyPoints)
127 {
128 G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
129 G4endl;
130 G4cout << "Trying to register more points than originally declared" << G4endl;
131 return;
132 }
133 G4double logEne = std::log(energy);
134
135 //XS0
136 G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
137 theVector->PutValue(binNumber,logEne,val);
138
139 //XS1
140 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
141 val = std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
142 theVector->PutValue(binNumber,logEne,val);
143
144 //XS2
145 theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
146 val = std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
147 theVector->PutValue(binNumber,logEne,val);
148
149 //XH0
150 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
151 val = std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
152 theVector->PutValue(binNumber,logEne,val);
153
154 //XH1
155 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
156 val = std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
157 theVector->PutValue(binNumber,logEne,val);
158
159 //XH2
160 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
161 val = std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
162 theVector->PutValue(binNumber,logEne,val);
163
164 return;
165}
166
167//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
168
169void G4PenelopeCrossSection::AddShellCrossSectionPoint(size_t binNumber,
170 size_t shellID,
171 G4double energy,
172 G4double xs)
173{
174 if (!shellCrossSections)
175 {
176 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
177 G4endl;
178 G4cout << "Trying to fill un-initialized table" << G4endl;
179 return;
180 }
181
182 if (shellID >= numberOfShells)
183 {
184 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
185 G4endl;
186 G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
187 << numberOfShells-1 << G4endl;
188 return;
189 }
190
191 //fill vector
192 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
193
194 if (binNumber >= numberOfEnergyPoints)
195 {
196 G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
197 G4endl;
198 G4cout << "Trying to register more points than originally declared" << G4endl;
199 return;
200 }
201 G4double logEne = std::log(energy);
202 G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
203 theVector->PutValue(binNumber,logEne,val);
204
205 return;
206}
207
208//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
209
210G4double G4PenelopeCrossSection::GetTotalCrossSection(G4double energy)
211{
212 G4double result = 0;
213 //take here XS0 + XH0
214 if (!softCrossSections || !hardCrossSections)
215 {
216 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
217 G4endl;
218 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
219 return result;
220 }
221
222 // 1) soft part
223 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[0];
224 if (theVector->GetVectorLength() < numberOfEnergyPoints)
225 {
226 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
227 G4endl;
228 G4cout << "Soft cross section table looks not filled" << G4endl;
229 return result;
230 }
231 G4double logene = std::log(energy);
232 G4double logXS = theVector->Value(logene);
233 G4double softXS = std::exp(logXS);
234
235 // 2) hard part
236 theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
237 if (theVector->GetVectorLength() < numberOfEnergyPoints)
238 {
239 G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
240 G4endl;
241 G4cout << "Hard cross section table looks not filled" << G4endl;
242 return result;
243 }
244 logXS = theVector->Value(logene);
245 G4double hardXS = std::exp(logXS);
246
247 result = hardXS + softXS;
248 return result;
249
250}
251
252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
253
254G4double G4PenelopeCrossSection::GetHardCrossSection(G4double energy)
255{
256 G4double result = 0;
257 //take here XH0
258 if (!hardCrossSections)
259 {
260 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
261 G4endl;
262 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
263 return result;
264 }
265
266 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
267 if (theVector->GetVectorLength() < numberOfEnergyPoints)
268 {
269 G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
270 G4endl;
271 G4cout << "Hard cross section table looks not filled" << G4endl;
272 return result;
273 }
274 G4double logene = std::log(energy);
275 G4double logXS = theVector->Value(logene);
276 result = std::exp(logXS);
277
278 return result;
279}
280
281
282//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
283
284G4double G4PenelopeCrossSection::GetSoftStoppingPower(G4double energy)
285{
286 G4double result = 0;
287 //take here XH0
288 if (!softCrossSections)
289 {
290 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
291 G4endl;
292 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
293 return result;
294 }
295
296 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
297 if (theVector->GetVectorLength() < numberOfEnergyPoints)
298 {
299 G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
300 G4endl;
301 G4cout << "Soft cross section table looks not filled" << G4endl;
302 return result;
303 }
304 G4double logene = std::log(energy);
305 G4double logXS = theVector->Value(logene);
306 result = std::exp(logXS);
307
308 return result;
309}
310
311//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
312
313G4double G4PenelopeCrossSection::GetShellCrossSection(size_t shellID,G4double energy)
314{
315 G4double result = 0;
316 if (!shellCrossSections)
317 {
318 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
319 G4endl;
320 G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
321 return result;
322 }
323 if (shellID >= numberOfShells)
324 {
325 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
326 G4endl;
327 G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
328 << numberOfShells-1 << G4endl;
329 return result;
330 }
331
332 G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
333
334 if (theVector->GetVectorLength() < numberOfEnergyPoints)
335 {
336 G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
337 G4endl;
338 G4cout << "Soft cross section table looks not filled" << G4endl;
339 return result;
340 }
341 G4double logene = std::log(energy);
342 G4double logXS = theVector->Value(logene);
343 result = std::exp(logXS);
344
345 return result;
346}
347
348//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
349
350void G4PenelopeCrossSection::NormalizeShellCrossSections()
351{
352 if (isNormalized) //already done!
353 {
354 G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
355 G4cout << "already invoked. Ignore it" << G4endl;
356 return;
357 }
358
359 for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
360 {
361 //energy grid is the same for all shells
362
363 //Recalculate manually the XS factor, to avoid problems with
364 //underflows
365 G4double normFactor = 0.;
366 for (size_t shellID=0;shellID<numberOfShells;shellID++)
367 {
368 G4PhysicsFreeVector* theVec =
369 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
370
371 normFactor += std::exp((*theVec)[i]);
372 }
373 G4double logNormFactor = std::log(normFactor);
374 //Normalize
375 for (size_t shellID=0;shellID<numberOfShells;shellID++)
376 {
377 G4PhysicsFreeVector* theVec =
378 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
379 G4double previousValue = (*theVec)[i]; //log(XS)
380 G4double logEnergy = theVec->GetLowEdgeEnergy(i);
381 //log(XS/normFactor) = log(XS) - log(normFactor)
382 theVec->PutValue(i,logEnergy,previousValue-logNormFactor);
383 }
384 }
385
386 isNormalized = true;
387
388
389 /*
390 //TESTING
391 for (size_t shellID=0;shellID<numberOfShells;shellID++)
392 {
393 G4cout << "SHELL " << shellID << G4endl;
394 G4PhysicsFreeVector* theVec =
395 (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
396 for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
397 {
398 G4double logene = theVec->GetLowEdgeEnergy(i);
399 G4cout << std::exp(logene)/MeV << " " << std::exp((*theVec)[i]) << G4endl;
400 }
401 }
402 */
403
404 return;
405}
Note: See TracBrowser for help on using the repository browser.