Geant4  10.02.p02
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 97613 2016-06-06 12:24:51Z 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 #include "G4Exp.hh"
43 
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
45 G4PenelopeCrossSection::G4PenelopeCrossSection(size_t nPointsE,size_t nShells) :
46  numberOfEnergyPoints(nPointsE),numberOfShells(nShells),softCrossSections(0),
47  hardCrossSections(0),shellCrossSections(0),shellNormalizedCrossSections(0)
48 {
49  //check the number of points is not zero
51  {
53  ed << "G4PenelopeCrossSection: invalid number of energy points " << G4endl;
54  G4Exception("G4PenelopeCrossSection::G4PenelopeCrossSection()",
55  "em2017",FatalException,ed);
56  }
57 
58  isNormalized = false;
59 
60  // 1) soft XS table
62  //the table contains 3 G4PhysicsFreeVectors,
63  //(softCrossSections)[0] --> log XS0 vs. log E
64  //(softCrossSections)[1] --> log XS1 vs. log E
65  //(softCrossSections)[2] --> log XS2 vs. log E
66 
67  //everything is log-log
68  for (size_t i=0;i<3;i++)
70 
71  //2) hard XS table
73  //the table contains 3 G4PhysicsFreeVectors,
74  //(hardCrossSections)[0] --> log XH0 vs. log E
75  //(hardCrossSections)[1] --> log XH1 vs. log E
76  //(hardCrossSections)[2] --> log XH2 vs. log E
77 
78  //everything is log-log
79  for (size_t i=0;i<3;i++)
81 
82  //3) shell XS table, if it is the case
83  if (numberOfShells)
84  {
87  //the table has to contain numberofShells G4PhysicsFreeVectors,
88  //(theTable)[ishell] --> cross section for shell #ishell
89  for (size_t i=0;i<numberOfShells;i++)
90  {
93  }
94  }
95 }
96 
97 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
99 {
100  //clean up tables
101  if (shellCrossSections)
102  {
103  //shellCrossSections->clearAndDestroy();
104  delete shellCrossSections;
105  }
107  {
108  //shellNormalizedCrossSections->clearAndDestroy();
110  }
111  if (softCrossSections)
112  {
113  //softCrossSections->clearAndDestroy();
114  delete softCrossSections;
115  }
116  if (hardCrossSections)
117  {
118  //hardCrossSections->clearAndDestroy();
119  delete hardCrossSections;
120  }
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
125  G4double XH0,
126  G4double XH1, G4double XH2,
127  G4double XS0, G4double XS1,
128  G4double XS2)
129 {
131  {
132  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
133  G4endl;
134  G4cout << "Trying to fill un-initialized tables" << G4endl;
135  return;
136  }
137 
138  //fill vectors
140 
141  if (binNumber >= numberOfEnergyPoints)
142  {
143  G4cout << "Something wrong in G4PenelopeCrossSection::AddCrossSectionPoint" <<
144  G4endl;
145  G4cout << "Trying to register more points than originally declared" << G4endl;
146  return;
147  }
148  G4double logEne = std::log(energy);
149 
150  //XS0
151  G4double val = std::log(std::max(XS0,1e-42*cm2)); //avoid log(0)
152  theVector->PutValue(binNumber,logEne,val);
153 
154  //XS1
155  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[1];
156  val = std::log(std::max(XS1,1e-42*eV*cm2)); //avoid log(0)
157  theVector->PutValue(binNumber,logEne,val);
158 
159  //XS2
160  theVector = (G4PhysicsFreeVector*) (*softCrossSections)[2];
161  val = std::log(std::max(XS2,1e-42*eV*eV*cm2)); //avoid log(0)
162  theVector->PutValue(binNumber,logEne,val);
163 
164  //XH0
165  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
166  val = std::log(std::max(XH0,1e-42*cm2)); //avoid log(0)
167  theVector->PutValue(binNumber,logEne,val);
168 
169  //XH1
170  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[1];
171  val = std::log(std::max(XH1,1e-42*eV*cm2)); //avoid log(0)
172  theVector->PutValue(binNumber,logEne,val);
173 
174  //XH2
175  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[2];
176  val = std::log(std::max(XH2,1e-42*eV*eV*cm2)); //avoid log(0)
177  theVector->PutValue(binNumber,logEne,val);
178 
179  return;
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
183 
185  size_t shellID,
187  G4double xs)
188 {
189  if (!shellCrossSections)
190  {
191  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
192  G4endl;
193  G4cout << "Trying to fill un-initialized table" << G4endl;
194  return;
195  }
196 
197  if (shellID >= numberOfShells)
198  {
199  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
200  G4endl;
201  G4cout << "Trying to fill shell #" << shellID << " while the maximum is "
202  << numberOfShells-1 << G4endl;
203  return;
204  }
205 
206  //fill vector
207  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
208 
209  if (binNumber >= numberOfEnergyPoints)
210  {
211  G4cout << "Something wrong in G4PenelopeCrossSection::AddShellCrossSectionPoint" <<
212  G4endl;
213  G4cout << "Trying to register more points than originally declared" << G4endl;
214  return;
215  }
216  G4double logEne = std::log(energy);
217  G4double val = std::log(std::max(xs,1e-42*cm2)); //avoid log(0)
218  theVector->PutValue(binNumber,logEne,val);
219 
220  return;
221 }
222 
223 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
224 
226 {
227  G4double result = 0;
228  //take here XS0 + XH0
230  {
231  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
232  G4endl;
233  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
234  return result;
235  }
236 
237  // 1) soft part
239  if (theVector->GetVectorLength() < numberOfEnergyPoints)
240  {
241  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
242  G4endl;
243  G4cout << "Soft cross section table looks not filled" << G4endl;
244  return result;
245  }
246  G4double logene = std::log(energy);
247  G4double logXS = theVector->Value(logene);
248  G4double softXS = G4Exp(logXS);
249 
250  // 2) hard part
251  theVector = (G4PhysicsFreeVector*) (*hardCrossSections)[0];
252  if (theVector->GetVectorLength() < numberOfEnergyPoints)
253  {
254  G4cout << "Something wrong in G4PenelopeCrossSection::GetTotalCrossSection" <<
255  G4endl;
256  G4cout << "Hard cross section table looks not filled" << G4endl;
257  return result;
258  }
259  logXS = theVector->Value(logene);
260  G4double hardXS = G4Exp(logXS);
261 
262  result = hardXS + softXS;
263  return result;
264 
265 }
266 
267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
268 
270 {
271  G4double result = 0;
272  //take here XH0
273  if (!hardCrossSections)
274  {
275  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
276  G4endl;
277  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
278  return result;
279  }
280 
282  if (theVector->GetVectorLength() < numberOfEnergyPoints)
283  {
284  G4cout << "Something wrong in G4PenelopeCrossSection::GetHardCrossSection" <<
285  G4endl;
286  G4cout << "Hard cross section table looks not filled" << G4endl;
287  return result;
288  }
289  G4double logene = std::log(energy);
290  G4double logXS = theVector->Value(logene);
291  result = G4Exp(logXS);
292 
293  return result;
294 }
295 
296 
297 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
298 
300 {
301  G4double result = 0;
302  //take here XH0
303  if (!softCrossSections)
304  {
305  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
306  G4endl;
307  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
308  return result;
309  }
310 
312  if (theVector->GetVectorLength() < numberOfEnergyPoints)
313  {
314  G4cout << "Something wrong in G4PenelopeCrossSection::GetSoftStoppingPower" <<
315  G4endl;
316  G4cout << "Soft cross section table looks not filled" << G4endl;
317  return result;
318  }
319  G4double logene = std::log(energy);
320  G4double logXS = theVector->Value(logene);
321  result = G4Exp(logXS);
322 
323  return result;
324 }
325 
326 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
327 
329 {
330  G4double result = 0;
331  if (!shellCrossSections)
332  {
333  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
334  G4endl;
335  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
336  return result;
337  }
338  if (shellID >= numberOfShells)
339  {
340  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
341  G4endl;
342  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
343  << numberOfShells-1 << G4endl;
344  return result;
345  }
346 
347  G4PhysicsFreeVector* theVector = (G4PhysicsFreeVector*) (*shellCrossSections)[shellID];
348 
349  if (theVector->GetVectorLength() < numberOfEnergyPoints)
350  {
351  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
352  G4endl;
353  G4cout << "Shell cross section table looks not filled" << G4endl;
354  return result;
355  }
356  G4double logene = std::log(energy);
357  G4double logXS = theVector->Value(logene);
358  result = G4Exp(logXS);
359 
360  return result;
361 }
362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
363 
365 {
366  G4double result = 0;
368  {
369  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
370  G4endl;
371  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
372  return result;
373  }
374 
375  if (!isNormalized)
376  {
377  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" << G4endl;
378  G4cout << "The table of normalized cross section is not initialized" << G4endl;
379  }
380 
381 
382  if (shellID >= numberOfShells)
383  {
384  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
385  G4endl;
386  G4cout << "Trying to retrieve shell #" << shellID << " while the maximum is "
387  << numberOfShells-1 << G4endl;
388  return result;
389  }
390 
391  const G4PhysicsFreeVector* theVector =
393 
394  if (theVector->GetVectorLength() < numberOfEnergyPoints)
395  {
396  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
397  G4endl;
398  G4cout << "Shell cross section table looks not filled" << G4endl;
399  return result;
400  }
401  G4double logene = std::log(energy);
402  G4double logXS = theVector->Value(logene);
403  result = G4Exp(logXS);
404 
405  return result;
406 }
407 
408 
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
411 
413 {
414  if (isNormalized) //already done!
415  {
416  G4cout << "G4PenelopeCrossSection::NormalizeShellCrossSections()" << G4endl;
417  G4cout << "already invoked. Ignore it" << G4endl;
418  return;
419  }
420 
422  {
423  G4cout << "Something wrong in G4PenelopeCrossSection::GetShellCrossSection" <<
424  G4endl;
425  G4cout << "Trying to retrieve from un-initialized tables" << G4endl;
426  return;
427  }
428 
429  for (size_t i=0;i<numberOfEnergyPoints;i++) //loop on energy
430  {
431  //energy grid is the same for all shells
432 
433  //Recalculate manually the XS factor, to avoid problems with
434  //underflows
435  G4double normFactor = 0.;
436  for (size_t shellID=0;shellID<numberOfShells;shellID++)
437  {
438  G4PhysicsFreeVector* theVec =
440 
441  normFactor += G4Exp((*theVec)[i]);
442  }
443  G4double logNormFactor = std::log(normFactor);
444  //Normalize
445  for (size_t shellID=0;shellID<numberOfShells;shellID++)
446  {
447  G4PhysicsFreeVector* theVec =
449  G4PhysicsFreeVector* theFullVec =
451  G4double previousValue = (*theFullVec)[i]; //log(XS)
452  G4double logEnergy = theFullVec->GetLowEdgeEnergy(i);
453  //log(XS/normFactor) = log(XS) - log(normFactor)
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;
472  }
473  }
474  */
475 
476  return;
477 }
static const double cm2
Definition: G4SIunits.hh:119
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
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const double eV
Definition: G4SIunits.hh:212
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)