Geant4  10.01.p03
G4PenelopeCrossSection.cc
Go to the documentation of this file.
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 76220 2013-11-08 10:15:00Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 18 Mar 2010 L Pandola First implementation
33 // 09 Mar 2012 L. Pandola Add public method (and machinery) to return
34 // the absolute and the normalized shell cross
35 // sections independently.
36 //
38 #include "G4SystemOfUnits.hh"
39 #include "G4PhysicsTable.hh"
40 #include "G4PhysicsFreeVector.hh"
41 #include "G4DataVector.hh"
42 
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
44 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
45  numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
46  hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
47 {
48  //check the number of points is not zero
50  {
52  ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
53  G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
54  "em2017",FatalException,ed);
55  }
56 
57  isNormalized = false;
58 
59  // 1) soft XS table
61  //the table contains 3 G4PhysicsFreeVectors,
62  //(softCrossSections)[0] --> log XS0 vs. log E
63  //(softCrossSections)[1] --> log XS1 vs. log E
64  //(softCrossSections)[2] --> log XS2 vs. log E
65 
66  //everything is log-log
67  for (size_t i=0;i<3;i++)
69 
70  //2) hard XS table
72  //the table contains 3 G4PhysicsFreeVectors,
73  //(hardCrossSections)[0] --> log XH0 vs. log E
74  //(hardCrossSections)[1] --> log XH1 vs. log E
75  //(hardCrossSections)[2] --> log XH2 vs. log E
76 
77  //everything is log-log
78  for (size_t i=0;i<3;i++)
80 
81  //3) shell XS table, if it is the case
82  if (numberOfShells)
83  {
86  //the table has to contain numberofShells G4PhysicsFreeVectors,
87  //(theTable)[ishell] --> cross section for shell #ishell
88  for (size_t i=0;i<numberOfShells;i++)
89  {
92  }
93  }
94 }
95 
96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
98 {
99  //clean up tables
100  if (shellCrossSections)
101  {
102  //shellCrossSections->clearAndDestroy();
103  delete shellCrossSections;
104  }
106  {
107  //shellNormalizedCrossSections->clearAndDestroy();
109  }
110  if (softCrossSections)
111  {
112  //softCrossSections->clearAndDestroy();
113  delete softCrossSections;
114  }
115  if (hardCrossSections)
116  {
117  //hardCrossSections->clearAndDestroy();
118  delete hardCrossSections;
119  }
120 }
121 
122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
124  G4double XH0,
125  G4double XH1, G4double XH2,
126  G4double XS0, G4double XS1,
127  G4double XS2)
128 {
130  {
131  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
132  G4endl;
133  G4cout << "Trying to fill un-initialized tables" << G4endl;
134  return;
135  }
136 
137  //fill vectors
139 
140  if (binNumber >= numberOfEnergyPoints)
141  {
142  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
143  G4endl;
144  G4cout << "Trying to register more points than originally declared" << G4endl;
145  return;
146  }
147  G4double logEne = std::log(energy);
148 
149  //XS0
150  G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
151  theVector->PutValue(binNumber,logEne,val);
152 
153  //XS1
154  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
155  val = std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
156  theVector->PutValue(binNumber,logEne,val);
157 
158  //XS2
159  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
160  val = std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
161  theVector->PutValue(binNumber,logEne,val);
162 
163  //XH0
164  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
165  val = std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
166  theVector->PutValue(binNumber,logEne,val);
167 
168  //XH1
169  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
170  val = std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
171  theVector->PutValue(binNumber,logEne,val);
172 
173  //XH2
174  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
175  val = std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
176  theVector->PutValue(binNumber,logEne,val);
177 
178  return;
179 }
180 
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
182 
184  size_t shellID,
186  G4double xs)
187 {
188  if (!shellCrossSections)
189  {
190  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
191  G4endl;
192  G4cout << "Trying to fill un-initialized table" << G4endl;
193  return;
194  }
195 
196  if (shellID >= numberOfShells)
197  {
198  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
199  G4endl;
200  G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
201  << numberOfShells-1 << G4endl;
202  return;
203  }
204 
205  //fill vector
206  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
207 
208  if (binNumber >= numberOfEnergyPoints)
209  {
210  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
211  G4endl;
212  G4cout << "Trying to register more points than originally declared" << G4endl;
213  return;
214  }
215  G4double logEne = std::log(energy);
216  G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
217  theVector->PutValue(binNumber,logEne,val);
218 
219  return;
220 }
221 
222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
223 
225 {
226  G4double result = 0;
227  //take here XS0 + XH0
229  {
230  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
231  G4endl;
232  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
233  return result;
234  }
235 
236  // 1) soft part
238  if (theVector->GetVectorLength() < numberOfEnergyPoints)
239  {
240  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
241  G4endl;
242  G4cout << "Soft cross section table looks not filled" << G4endl;
243  return result;
244  }
245  G4double logene = std::log(energy);
246  G4double logXS = theVector->Value(logene);
247  G4double softXS = std::exp(logXS);
248 
249  // 2) hard part
250  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
251  if (theVector->GetVectorLength() < numberOfEnergyPoints)
252  {
253  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
254  G4endl;
255  G4cout << "Hard cross section table looks not filled" << G4endl;
256  return result;
257  }
258  logXS = theVector->Value(logene);
259  G4double hardXS = std::exp(logXS);
260 
261  result = hardXS + softXS;
262  return result;
263 
264 }
265 
266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
267 
269 {
270  G4double result = 0;
271  //take here XH0
272  if (!hardCrossSections)
273  {
274  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
275  G4endl;
276  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
277  return result;
278  }
279 
281  if (theVector->GetVectorLength() < numberOfEnergyPoints)
282  {
283  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
284  G4endl;
285  G4cout << "Hard cross section table looks not filled" << G4endl;
286  return result;
287  }
288  G4double logene = std::log(energy);
289  G4double logXS = theVector->Value(logene);
290  result = std::exp(logXS);
291 
292  return result;
293 }
294 
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
297 
299 {
300  G4double result = 0;
301  //take here XH0
302  if (!softCrossSections)
303  {
304  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
305  G4endl;
306  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
307  return result;
308  }
309 
311  if (theVector->GetVectorLength() < numberOfEnergyPoints)
312  {
313  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
314  G4endl;
315  G4cout << "Soft cross section table looks not filled" << G4endl;
316  return result;
317  }
318  G4double logene = std::log(energy);
319  G4double logXS = theVector->Value(logene);
320  result = std::exp(logXS);
321 
322  return result;
323 }
324 
325 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
326 
328 {
329  G4double result = 0;
330  if (!shellCrossSections)
331  {
332  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
333  G4endl;
334  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
335  return result;
336  }
337  if (shellID >= numberOfShells)
338  {
339  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
340  G4endl;
341  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
342  << numberOfShells-1 << G4endl;
343  return result;
344  }
345 
346  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
347 
348  if (theVector->GetVectorLength() < numberOfEnergyPoints)
349  {
350  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
351  G4endl;
352  G4cout << "Shell cross section table looks not filled" << G4endl;
353  return result;
354  }
355  G4double logene = std::log(energy);
356  G4double logXS = theVector->Value(logene);
357  result = std::exp(logXS);
358 
359  return result;
360 }
361 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
362 
364 {
365  G4double result = 0;
367  {
368  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
369  G4endl;
370  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
371  return result;
372  }
373 
374  if (!isNormalized)
375  {
376  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
377  G4cout << "The table of normalized cross section is not initialized" << G4endl;
378  }
379 
380 
381  if (shellID >= numberOfShells)
382  {
383  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
384  G4endl;
385  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
386  << numberOfShells-1 << G4endl;
387  return result;
388  }
389 
390  const G4PhysicsFreeVector* theVector =
392 
393  if (theVector->GetVectorLength() < numberOfEnergyPoints)
394  {
395  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
396  G4endl;
397  G4cout << "Shell cross section table looks not filled" << G4endl;
398  return result;
399  }
400  G4double logene = std::log(energy);
401  G4double logXS = theVector->Value(logene);
402  result = std::exp(logXS);
403 
404  return result;
405 }
406 
407 
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
410 
412 {
413  if (isNormalized) //already done!
414  {
415  G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
416  G4cout << "already invoked. Ignore it" << G4endl;
417  return;
418  }
419 
421  {
422  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
423  G4endl;
424  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
425  return;
426  }
427 
428  for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
429  {
430  //energy grid is the same for all shells
431 
432  //Recalculate manually the XS factor, to avoid problems with
433  //underflows
434  G4double normFactor = 0.;
435  for (size_t shellID=0;shellID<numberOfShells;shellID++)
436  {
437  G4PhysicsFreeVector* theVec =
439 
440  normFactor += std::exp((*theVec)[i]);
441  }
442  G4double logNormFactor = std::log(normFactor);
443  //Normalize
444  for (size_t shellID=0;shellID<numberOfShells;shellID++)
445  {
446  G4PhysicsFreeVector* theVec =
448  G4PhysicsFreeVector* theFullVec =
449  (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
450  G4double previousValue = (*theFullVec)[i]; //log(XS)
451  G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
452  //log(XS/normFactor) = log(XS) - log(normFactor)
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 << std::exp(logene)/MeV << " " << std::exp((*theVec)[i]) << G4endl;
471  }
472  }
473  */
474 
475  return;
476 }
static const double cm2
Definition: G4SIunits.hh:107
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void push_back(G4PhysicsVector *)
void AddShellCrossSectionPoint(size_t binNumber, size_t shellID, G4double energy, G4double xs)
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
size_t GetVectorLength() const
G4double GetLowEdgeEnergy(size_t binNumber) const
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
Public interface for the master thread.
G4GLOB_DLL std::ostream G4cout
G4double GetTotalCrossSection(G4double energy) const
Returns total cross section at the given energy.
G4PhysicsTable * shellNormalizedCrossSections
G4double Value(G4double theEnergy, size_t &lastidx) const
G4PenelopeCrossSection(size_t nOfEnergyPoints, size_t nOfShells=0)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const double eV
Definition: G4SIunits.hh:194
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetNormalizedShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (normalized to 1)
G4double GetShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (per molecule)