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

Last change on this file since 1353 was 1350, checked in by garnier, 14 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.