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