Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4PenelopeCrossSection.cc
이 파일의 문서화 페이지로 가기
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 95950 2016-03-03 10:42:48Z 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 }
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
G4double GetLowEdgeEnergy(size_t binNumber) const
#define G4endl
Definition: G4ios.hh:61
G4double GetTotalCrossSection(G4double energy) const
Returns total cross section at the given energy.
Float_t normFactor
G4PenelopeCrossSection(size_t nOfEnergyPoints, size_t nOfShells=0)
G4double Value(G4double theEnergy, size_t &lastidx) const
void AddShellCrossSectionPoint(size_t binNumber, size_t shellID, G4double energy, G4double xs)
G4double GetShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (per molecule)
G4PhysicsTable * shellNormalizedCrossSections
double G4double
Definition: G4Types.hh:76
static constexpr double cm2
Definition: G4SIunits.hh:120
void AddCrossSectionPoint(size_t binNumber, G4double energy, G4double XH0, G4double XH1, G4double XH2, G4double XS0, G4double XS1, G4double XS2)
double energy
Definition: plottest35.C:25
static constexpr double eV
Definition: G4SIunits.hh:215
G4double G4ParticleHPJENDLHEData::G4double result
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
void push_back(G4PhysicsVector *)
G4double GetNormalizedShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (normalized to 1)
G4GLOB_DLL std::ostream G4cout
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
size_t GetVectorLength() const