Geant4  10.03
G4Element.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 //
27 // $Id: G4Element.cc 96794 2016-05-09 10:09:30Z gcosmo $
28 //
29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
30 
31 // 26-06-96: Code uses operators (+=, *=, ++, -> etc.) correctly, P. Urban
32 // 09-07-96: new data members added by L.Urban
33 // 17-01-97: aesthetic rearrangement, M.Maire
34 // 20-01-97: Compute Tsai's formula for the rad length, M.Maire
35 // 21-01-97: remove mixture flag, M.Maire
36 // 24-01-97: ComputeIonisationParameters().
37 // new data member: fTaul, M.Maire
38 // 29-01-97: Forbidden to create Element with Z<1 or N<Z, M.Maire
39 // 20-03-97: corrected initialization of pointers, M.Maire
40 // 28-04-98: atomic subshell binding energies stuff, V. Grichine
41 // 09-07-98: Ionisation parameters removed from the class, M.Maire
42 // 16-11-98: name Subshell -> Shell; GetBindingEnergy() (mma)
43 // 09-03-01: assignement operator revised (mma)
44 // 02-05-01: check identical Z in AddIsotope (marc)
45 // 03-05-01: flux.precision(prec) at begin/end of operator<<
46 // 13-09-01: suppression of the data member fIndexInTable
47 // 14-09-01: fCountUse: nb of materials which use this element
48 // 26-02-02: fIndexInTable renewed
49 // 30-03-05: warning in GetElement(elementName)
50 // 15-11-05: GetElement(elementName, G4bool warning=true)
51 // 17-10-06: Add fNaturalAbundances (V.Ivanchenko)
52 // 27-07-07: improve destructor (V.Ivanchenko)
53 // 18-10-07: move definition of material index to ComputeDerivedQuantities (VI)
54 // 25.10.11: new scheme for G4Exception (mma)
55 // 05-03-12: always create isotope vector (V.Ivanchenko)
56 
57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 
59 #include <iomanip>
60 #include <sstream>
61 
62 #include "G4Element.hh"
63 #include "G4AtomicShells.hh"
64 #include "G4NistManager.hh"
65 #include "G4PhysicalConstants.hh"
66 #include "G4SystemOfUnits.hh"
67 #include "G4Log.hh"
68 
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
73 // Constructor to Generate an element from scratch
74 
76  G4double zeff, G4double aeff)
77  : fName(name), fSymbol(symbol)
78 {
79  G4int iz = G4lrint(zeff);
80  if (iz < 1) {
82  ed << "Fail to create G4Element " << name
83  << " Z= " << zeff << " < 1 !";
84  G4Exception ("G4Element::G4Element()", "mat011", FatalException, ed);
85  }
86  if (std::abs(zeff - iz) > perMillion) {
88  ed << "G4Element Warning: " << name << " Z= " << zeff
89  << " A= " << aeff/(g/mole);
90  G4Exception("G4Element::G4Element()", "mat017", JustWarning, ed);
91  }
92 
94 
95  fZeff = zeff;
96  fAeff = aeff;
97  fNeff = fAeff/(g/mole);
98 
99  if(fNeff < 1.0) fNeff = 1.0;
100 
101  if (fNeff < zeff) {
103  ed << "Fail to create G4Element " << name
104  << " with Z= " << zeff << " N= " << fNeff
105  << " N < Z is not allowed" << G4endl;
106  G4Exception("G4Element::G4Element()", "mat012", FatalException, ed);
107  }
108 
112 
114 
115  for (G4int i=0;i<fNbOfAtomicShells;i++)
116  {
119  }
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
125 // Constructor to Generate element from a List of 'nIsotopes' isotopes, added
126 // via AddIsotope
127 
129  const G4String& symbol, G4int nIsotopes)
130  : fName(name),fSymbol(symbol)
131 {
133 
134  size_t n = size_t(nIsotopes);
135 
136  if(0 >= nIsotopes) {
138  ed << "Fail to create G4Element " << name
139  << " <" << symbol << "> with " << nIsotopes
140  << " isotopes";
141  G4Exception ("G4Element::G4Element()", "mat012", FatalException, ed);
142  } else {
144  fRelativeAbundanceVector = new G4double[nIsotopes];
145  }
146 }
147 
148 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
149 
150 // Add an isotope to the element
151 
152 void G4Element::AddIsotope(G4Isotope* isotope, G4double abundance)
153 {
154  if (theIsotopeVector == 0) {
156  ed << "Fail to add Isotope to G4Element " << fName
157  << " with Z= " << fZeff << " N= " << fNeff;
158  G4Exception ("G4Element::AddIsotope()", "mat013", FatalException, ed);
159  return;
160  }
161  G4int iz = isotope->GetZ();
162 
163  // filling ...
164  if ( fNumberOfIsotopes < (G4int)theIsotopeVector->size() ) {
165  // check same Z
166  if (fNumberOfIsotopes==0) { fZeff = G4double(iz); }
167  else if (G4double(iz) != fZeff) {
169  ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
170  << " with different Z= " << fZeff << fNeff;
171  G4Exception ("G4Element::AddIsotope()", "mat014", FatalException, ed);
172  return;
173  }
174  //Z ok
176  (*theIsotopeVector)[fNumberOfIsotopes] = isotope;
178 
179  } else {
181  ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
182  << " - more isotopes than declaired ";
183  G4Exception ("G4Element::AddIsotope()", "mat015", FatalException, ed);
184  return;
185  }
186 
187  // filled.
188  if ( fNumberOfIsotopes == (G4int)theIsotopeVector->size() ) {
189  G4double wtSum=0.0;
190  fAeff = 0.0;
191  for (G4int i=0; i<fNumberOfIsotopes; i++) {
192  fAeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetA();
193  wtSum += fRelativeAbundanceVector[i];
194  }
195  if(wtSum > 0.0) { fAeff /= wtSum; }
196  fNeff = fAeff/(g/mole);
197 
198  if(wtSum != 1.0) {
199  for(G4int i=0; i<fNumberOfIsotopes; ++i) {
200  fRelativeAbundanceVector[i] /= wtSum;
201  }
202  }
203 
207 
208  for ( G4int j = 0; j < fNbOfAtomicShells; j++ )
209  {
212  }
214  }
215 }
216 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 
220 {
221  theIsotopeVector = nullptr;
222  fRelativeAbundanceVector = nullptr;
223  fAtomicShells = nullptr;
224  fNbOfShellElectrons = nullptr;
225  fIonisation = nullptr;
226  fNumberOfIsotopes = 0;
227  fNaturalAbundance = false;
228 
229  // add initialisation of all remaining members
230  fZeff = 0;
231  fNeff = 0;
232  fAeff = 0;
233  fNbOfAtomicShells = 0;
234  fIndexInTable = 0;
235  fCoulomb = 0.0;
236  fRadTsai = 0.0;
237  fZ = 0;
238 }
239 
240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
241 
242 // Fake default constructor - sets only member data and allocates memory
243 // for usage restricted to object persistency
244 
246  : fZeff(0), fNeff(0), fAeff(0)
247 {
249 }
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252 
254 {
255  if (theIsotopeVector) { delete theIsotopeVector; }
257  if (fAtomicShells) { delete [] fAtomicShells; }
258  if (fNbOfShellElectrons) { delete [] fNbOfShellElectrons; }
259  if (fIonisation) { delete fIonisation; }
260 
261  //remove this element from theElementTable
263 }
264 
265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
266 
268 {
269  // some basic functions of the atomic number
270 
271  // Store in table
272  theElementTable.push_back(this);
273  fIndexInTable = theElementTable.size() - 1;
274 
275  // Radiation Length
278 
279  // parameters for energy loss by ionisation
280  if (fIonisation) { delete fIonisation; }
282  fZ = G4lrint(fZeff);
283 }
284 
285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
286 
288 {
289  //
290  // Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
291 
292  static const G4double k1 = 0.0083 , k2 = 0.20206 ,k3 = 0.0020 , k4 = 0.0369 ;
293 
294  G4double az2 = (fine_structure_const*fZeff)*(fine_structure_const*fZeff);
295  G4double az4 = az2 * az2;
296 
297  fCoulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
301 
303 {
304  //
305  // Compute Tsai's Expression for the Radiation Length
306  // (Phys Rev. D50 3-1 (1994) page 1254)
307 
308  static const G4double Lrad_light[] = {5.31 , 4.79 , 4.74 , 4.71} ;
309  static const G4double Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
310 
311  const G4double logZ3 = G4Log(fZeff)/3.;
312 
313  G4double Lrad, Lprad;
314  G4int iz = G4lrint(fZeff) - 1 ;
315  static const G4double log184 = G4Log(184.15);
316  static const G4double log1194 = G4Log(1194.);
317  if (iz <= 3) { Lrad = Lrad_light[iz] ; Lprad = Lprad_light[iz] ; }
318  else { Lrad = log184 - logZ3 ; Lprad = log1194 - 2*logZ3;}
319 
320  fRadTsai = 4*alpha_rcl2*fZeff*(fZeff*(Lrad-fCoulomb) + Lprad);
321 }
322 
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
324 
326 {
327  G4int Z = G4lrint(fZeff);
329  G4int n = nist->GetNumberOfNistIsotopes(Z);
330  G4int N0 = nist->GetNistFirstIsotopeN(Z);
331 
332  if("" == fSymbol) {
333  const std::vector<G4String> elmnames =
335  if(Z < (G4int)elmnames.size()) { fSymbol = elmnames[Z]; }
336  else { fSymbol = fName; }
337  }
338 
339  fNumberOfIsotopes = 0;
340  for(G4int i=0; i<n; ++i) {
341  if(nist->GetIsotopeAbundance(Z, N0+i) > 0.0) { ++fNumberOfIsotopes; }
342  }
343  theIsotopeVector = new G4IsotopeVector((unsigned int)fNumberOfIsotopes,0);
345  G4int idx = 0;
346  G4double xsum = 0.0;
347  for(G4int i=0; i<n; ++i) {
348  G4int N = N0 + i;
349  G4double x = nist->GetIsotopeAbundance(Z, N);
350  if(x > 0.0) {
351  std::ostringstream strm;
352  strm << fSymbol << N;
353  (*theIsotopeVector)[idx] = new G4Isotope(strm.str(), Z, N, 0.0, 0);
354  fRelativeAbundanceVector[idx] = x;
355  xsum += x;
356  ++idx;
357  }
358  }
359  if(xsum != 0.0 && xsum != 1.0) {
360  for(G4int i=0; i<idx; ++i) { fRelativeAbundanceVector[i] /= xsum; }
361  }
362  fNaturalAbundance = true;
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
366 
368 {
369  if (i<0 || i>=fNbOfAtomicShells) {
371  ed << "Invalid argument " << i << " in for G4Element " << fName
372  << " with Z= " << fZeff
373  << " and Nshells= " << fNbOfAtomicShells;
374  G4Exception("G4Element::GetAtomicShell()", "mat016", FatalException, ed);
375  return 0.0;
376  }
377  return fAtomicShells[i];
378 }
379 
380 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
381 
383 {
384  if (i<0 || i>=fNbOfAtomicShells) {
386  ed << "Invalid argument " << i << " for G4Element " << fName
387  << " with Z= " << fZeff
388  << " and Nshells= " << fNbOfAtomicShells;
389  G4Exception("G4Element::GetNbOfShellElectrons()", "mat016",
390  FatalException, ed);
391  return 0;
392  }
393  return fNbOfShellElectrons[i];
394 }
395 
396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
397 
399 {
400  return &theElementTable;
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 
406 {
407  return theElementTable.size();
408 }
409 
410 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
411 
413 {
414  // search the element by its name
415  for (size_t J=0; J<theElementTable.size(); ++J)
416  {
417  if (theElementTable[J]->GetName() == elementName)
418  return theElementTable[J];
419  }
420 
421  // the element does not exist in the table
422  if (warning) {
423  G4cout << "\n---> warning from G4Element::GetElement(). The element: "
424  << elementName << " does not exist in the table. Return NULL pointer."
425  << G4endl;
426  }
427  return nullptr;
428 }
429 
430 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
431 
432 std::ostream& operator<<(std::ostream& flux, const G4Element* element)
433 {
434  std::ios::fmtflags mode = flux.flags();
435  flux.setf(std::ios::fixed,std::ios::floatfield);
436  G4long prec = flux.precision(3);
437 
438  flux
439  << " Element: " << element->fName << " (" << element->fSymbol << ")"
440  << " Z = " << std::setw(4) << std::setprecision(1) << element->fZeff
441  << " N = " << std::setw(5) << std::setprecision(1)
442  << G4lrint(element->fNeff)
443  << " A = " << std::setw(6) << std::setprecision(3)
444  << (element->fAeff)/(g/mole) << " g/mole";
445 
446  for (G4int i=0; i<element->fNumberOfIsotopes; i++)
447  flux
448  << "\n ---> " << (*(element->theIsotopeVector))[i]
449  << " abundance: " << std::setw(6) << std::setprecision(3)
450  << (element->fRelativeAbundanceVector[i])/perCent << " %";
451 
452  flux.precision(prec);
453  flux.setf(mode,std::ios::floatfield);
454  return flux;
455 }
456 
457 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
458 
459  std::ostream& operator<<(std::ostream& flux, const G4Element& element)
460 {
461  flux << &element;
462  return flux;
463 }
464 
465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
466 
467 std::ostream& operator<<(std::ostream& flux, G4ElementTable ElementTable)
468 {
469  //Dump info for all known elements
470  flux << "\n***** Table : Nb of elements = " << ElementTable.size()
471  << " *****\n" << G4endl;
472 
473  for (size_t i=0; i<ElementTable.size(); i++) flux << ElementTable[i]
474  << G4endl << G4endl;
475 
476  return flux;
477 }
478 
479 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4double fRadTsai
Definition: G4Element.hh:258
static constexpr double perMillion
Definition: G4SIunits.hh:334
std::vector< G4Isotope * > G4IsotopeVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4String fName
Definition: G4AttUtils.hh:55
static G4Element * GetElement(G4String name, G4bool warning=true)
Definition: G4Element.cc:412
void ComputeCoulombFactor()
Definition: G4Element.cc:287
void AddNaturalIsotopes()
Definition: G4Element.cc:325
static G4ElementTable theElementTable
Definition: G4Element.hh:250
G4double fCoulomb
Definition: G4Element.hh:257
static constexpr double perCent
Definition: G4SIunits.hh:332
const std::vector< G4String > & GetNistElementNames() const
G4int * fNbOfShellElectrons
Definition: G4Element.hh:241
G4IsotopeVector * theIsotopeVector
Definition: G4Element.hh:245
long G4long
Definition: G4Types.hh:80
G4double * fAtomicShells
Definition: G4Element.hh:240
G4IonisParamElm * fIonisation
Definition: G4Element.hh:259
const char * name(G4int ptype)
G4int GetNbOfShellElectrons(G4int index) const
Definition: G4Element.cc:382
static constexpr double g
Definition: G4SIunits.hh:183
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4int GetNistFirstIsotopeN(G4int Z) const
G4bool fNaturalAbundance
Definition: G4Element.hh:252
static const double prec
Definition: RanecuEngine.cc:58
void ComputeLradTsaiFactor()
Definition: G4Element.cc:302
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
void InitializePointers()
Definition: G4Element.cc:219
bool G4bool
Definition: G4Types.hh:79
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
std::ostream & operator<<(std::ostream &flux, const G4Element *element)
Definition: G4Element.cc:432
void AddIsotope(G4Isotope *isotope, G4double RelativeAbundance)
Definition: G4Element.cc:152
G4Element(const G4String &name, const G4String &symbol, G4double Zeff, G4double Aeff)
Definition: G4Element.cc:75
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4int fZ
Definition: G4Element.hh:237
G4int GetNumberOfNistIsotopes(G4int Z) const
const G4int n
G4String fSymbol
Definition: G4Element.hh:233
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4int fNumberOfIsotopes
Definition: G4Element.hh:244
G4double fNeff
Definition: G4Element.hh:235
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4int GetZ() const
Definition: G4Isotope.hh:91
int G4lrint(double ad)
Definition: templates.hh:163
size_t fIndexInTable
Definition: G4Element.hh:251
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
G4String fName
Definition: G4Element.hh:232
virtual ~G4Element()
Definition: G4Element.cc:253
#define G4endl
Definition: G4ios.hh:61
G4int fNbOfAtomicShells
Definition: G4Element.hh:239
G4double fZeff
Definition: G4Element.hh:234
G4double * fRelativeAbundanceVector
Definition: G4Element.hh:246
G4double GetAtomicShell(G4int index) const
Definition: G4Element.cc:367
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
const G4String & GetName() const
Definition: G4Element.hh:127
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
static constexpr double mole
Definition: G4SIunits.hh:286
G4double fAeff
Definition: G4Element.hh:236
static G4int GetNumberOfShells(G4int Z)
void ComputeDerivedQuantities()
Definition: G4Element.cc:267