Geant4_10
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 69704 2013-05-13 09:06:12Z 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 
68 G4ElementTable G4Element::theElementTable;
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
71 
72 // Constructor to Generate an element from scratch
73 
75  G4double zeff, G4double aeff)
76  : fName(name), fSymbol(symbol)
77 {
78  G4int iz = G4lrint(zeff);
79  if (iz < 1) {
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) {
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  fNeff = aeff/(g/mole);
96  fAeff = aeff;
97 
98  if(fNeff < 1.0) fNeff = 1.0;
99 
100  if (fNeff < zeff) {
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 
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) {
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) {
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 < theIsotopeVector->size() ) {
164  // check same Z
165  if (fNumberOfIsotopes==0) { fZeff = G4double(iz); }
166  else if (G4double(iz) != fZeff) {
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 {
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 == theIsotopeVector->size() ) {
188  // Compute Neff, Aeff
189  G4double wtSum=0.0;
190 
191  G4double aeff = 0.0;
192  G4double neff = 0.0;
193  for (size_t i=0;i<fNumberOfIsotopes;i++) {
194  aeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetA();
195  neff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetN();
196  wtSum += fRelativeAbundanceVector[i];
197  }
198  aeff /= wtSum;
199  neff /= wtSum;
200  if(0.0 == fAeff) {
201  fAeff = aeff;
202  fNeff = neff;
203  }
204 
205  if(wtSum != 1.0) {
206  for(size_t i=0; i<fNumberOfIsotopes; ++i) {
207  fRelativeAbundanceVector[i] /= wtSum;
208  }
209  }
210 
211  fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
212  fAtomicShells = new G4double[fNbOfAtomicShells];
213  fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
214 
215  for ( G4int j = 0; j < fNbOfAtomicShells; j++ )
216  {
217  fAtomicShells[j] = G4AtomicShells::GetBindingEnergy(iz, j);
218  fNbOfShellElectrons[j] = G4AtomicShells::GetNumberOfElectrons(iz, j);
219  }
220  ComputeDerivedQuantities();
221 
222  }
223 }
224 
225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
226 
227 void G4Element::InitializePointers()
228 {
229  theIsotopeVector = 0;
230  fRelativeAbundanceVector = 0;
231  fAtomicShells = 0;
232  fNbOfShellElectrons = 0;
233  fIonisation = 0;
234  fNumberOfIsotopes = 0;
235  fNaturalAbundance = false;
236 
237  // add initialisation of all remaining members
238  fZeff = 0;
239  fNeff = 0;
240  fAeff = 0;
241  fNbOfAtomicShells = 0;
242  fIndexInTable = 0;
243  fCoulomb = 0.0;
244  fRadTsai = 0.0;
245 }
246 
247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
248 
249 // Fake default constructor - sets only member data and allocates memory
250 // for usage restricted to object persistency
251 
253  : fZeff(0), fNeff(0), fAeff(0)
254 {
255  InitializePointers();
256 }
257 
258 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
259 
261 {
262  // G4cout << "### Destruction of element " << fName << " started" <<G4endl;
263 
264  if (theIsotopeVector) { delete theIsotopeVector; }
265  if (fRelativeAbundanceVector) { delete [] fRelativeAbundanceVector; }
266  if (fAtomicShells) { delete [] fAtomicShells; }
267  if (fNbOfShellElectrons) { delete [] fNbOfShellElectrons; }
268  if (fIonisation) { delete fIonisation; }
269 
270  //remove this element from theElementTable
271  theElementTable[fIndexInTable] = 0;
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 
276 void G4Element::ComputeDerivedQuantities()
277 {
278  // some basic functions of the atomic number
279 
280  // Store in table
281  theElementTable.push_back(this);
282  fIndexInTable = theElementTable.size() - 1;
283 
284  // Radiation Length
285  ComputeCoulombFactor();
286  ComputeLradTsaiFactor();
287 
288  // parameters for energy loss by ionisation
289  if (fIonisation) { delete fIonisation; }
290  fIonisation = new G4IonisParamElm(fZeff);
291 }
292 
293 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
294 
295 void G4Element::ComputeCoulombFactor()
296 {
297  //
298  // Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
299 
300  static const G4double k1 = 0.0083 , k2 = 0.20206 ,k3 = 0.0020 , k4 = 0.0369 ;
301 
302  G4double az2 = (fine_structure_const*fZeff)*(fine_structure_const*fZeff);
303  G4double az4 = az2 * az2;
304 
305  fCoulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
306 }
307 
308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
309 
310 void G4Element::ComputeLradTsaiFactor()
311 {
312  //
313  // Compute Tsai's Expression for the Radiation Length
314  // (Phys Rev. D50 3-1 (1994) page 1254)
315 
316  static const G4double Lrad_light[] = {5.31 , 4.79 , 4.74 , 4.71} ;
317  static const G4double Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
318 
319  const G4double logZ3 = std::log(fZeff)/3.;
320 
321  G4double Lrad, Lprad;
322  G4int iz = (G4int)(fZeff+0.5) - 1 ;
323  if (iz <= 3) { Lrad = Lrad_light[iz] ; Lprad = Lprad_light[iz] ; }
324  else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}
325 
326  fRadTsai = 4*alpha_rcl2*fZeff*(fZeff*(Lrad-fCoulomb) + Lprad);
327 }
328 
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330 
331 void G4Element::AddNaturalIsotopes()
332 {
333  G4int Z = G4lrint(fZeff);
335  G4int n = nist->GetNumberOfNistIsotopes(Z);
336  G4int N0 = nist->GetNistFirstIsotopeN(Z);
337  fNumberOfIsotopes = 0;
338  for(G4int i=0; i<n; ++i) {
339  if(nist->GetIsotopeAbundance(Z, N0+i) > 0.0) { ++fNumberOfIsotopes; }
340  }
341  theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
342  fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
343  G4int idx = 0;
344  G4double xsum = 0.0;
345  for(G4int i=0; i<n; ++i) {
346  G4int N = N0 + i;
347  G4double x = nist->GetIsotopeAbundance(Z, N);
348  if(x > 0.0) {
349  std::ostringstream strm;
350  strm << fSymbol << N;
351  (*theIsotopeVector)[idx] = new G4Isotope(strm.str(),Z, N, 0.0, 0);
352  fRelativeAbundanceVector[idx] = x;
353  xsum += x;
354  ++idx;
355  }
356  }
357  if(xsum != 0.0 && xsum != 1.0) {
358  for(G4int i=0; i<idx; ++i) { fRelativeAbundanceVector[i] /= xsum; }
359  }
360  fNaturalAbundance = true;
361 }
362 
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
364 
366 {
367  if (i<0 || i>=fNbOfAtomicShells) {
369  ed << "Invalid argument " << i << " in for G4Element " << fName
370  << " with Z= " << fZeff
371  << " and Nshells= " << fNbOfAtomicShells;
372  G4Exception("G4Element::GetAtomicShell()", "mat016", FatalException, ed);
373  return 0.0;
374  }
375  return fAtomicShells[i];
376 }
377 
378 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
379 
381 {
382  if (i<0 || i>=fNbOfAtomicShells) {
384  ed << "Invalid argument " << i << " for G4Element " << fName
385  << " with Z= " << fZeff
386  << " and Nshells= " << fNbOfAtomicShells;
387  G4Exception("G4Element::GetNbOfShellElectrons()", "mat016", FatalException, ed);
388  return 0;
389  }
390  return fNbOfShellElectrons[i];
391 }
392 
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
394 
396 {
397  return &theElementTable;
398 }
399 
400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
401 
403 {
404  return theElementTable.size();
405 }
406 
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
408 
410 {
411  // search the element by its name
412  for (size_t J=0 ; J<theElementTable.size() ; J++)
413  {
414  if (theElementTable[J]->GetName() == elementName)
415  return theElementTable[J];
416  }
417 
418  // the element does not exist in the table
419  if (warning) {
420  G4cout << "\n---> warning from G4Element::GetElement(). The element: "
421  << elementName << " does not exist in the table. Return NULL pointer."
422  << G4endl;
423  }
424  return 0;
425 }
426 
427 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
428 
430 {
431  InitializePointers();
432  *this = right;
433 
434  // Store this new element in table and set the index
435  theElementTable.push_back(this);
436  fIndexInTable = theElementTable.size() - 1;
437 }
438 
439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440 
441 const G4Element& G4Element::operator=(const G4Element& right)
442 {
443  if (this != &right)
444  {
445  fName = right.fName;
446  fSymbol = right.fSymbol;
447  fZeff = right.fZeff;
448  fNeff = right.fNeff;
449  fAeff = right.fAeff;
450 
451  if (fAtomicShells) delete [] fAtomicShells;
452  fNbOfAtomicShells = right.fNbOfAtomicShells;
453  fAtomicShells = new G4double[fNbOfAtomicShells];
454 
455  if (fNbOfShellElectrons) delete [] fNbOfShellElectrons;
456  fNbOfAtomicShells = right.fNbOfAtomicShells;
457  fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
458 
459  for ( G4int i = 0; i < fNbOfAtomicShells; i++ )
460  {
461  fAtomicShells[i] = right.fAtomicShells[i];
462  fNbOfShellElectrons[i] = right.fNbOfShellElectrons[i];
463  }
464  if (theIsotopeVector) delete theIsotopeVector;
465  if (fRelativeAbundanceVector) delete [] fRelativeAbundanceVector;
466 
467  fNumberOfIsotopes = right.fNumberOfIsotopes;
468  if (fNumberOfIsotopes > 0)
469  {
470  theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
471  fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
472  for (size_t i=0;i<fNumberOfIsotopes;i++)
473  {
474  (*theIsotopeVector)[i] = (*right.theIsotopeVector)[i];
475  fRelativeAbundanceVector[i] = right.fRelativeAbundanceVector[i];
476  }
477  }
478  ComputeDerivedQuantities();
479  }
480  return *this;
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
484 
486 {
487  return (this == (G4Element*) &right);
488 }
489 
490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
491 
493 {
494  return (this != (G4Element*) &right);
495 }
496 
497 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
498 
499 std::ostream& operator<<(std::ostream& flux, G4Element* element)
500 {
501  std::ios::fmtflags mode = flux.flags();
502  flux.setf(std::ios::fixed,std::ios::floatfield);
503  G4long prec = flux.precision(3);
504 
505  flux
506  << " Element: " << element->fName << " (" << element->fSymbol << ")"
507  << " Z = " << std::setw(4) << std::setprecision(1) << element->fZeff
508  << " N = " << std::setw(5) << std::setprecision(1) << element->fNeff
509  << " A = " << std::setw(6) << std::setprecision(2)
510  << (element->fAeff)/(g/mole) << " g/mole";
511 
512  for (size_t i=0; i<element->fNumberOfIsotopes; i++)
513  flux
514  << "\n ---> " << (*(element->theIsotopeVector))[i]
515  << " abundance: " << std::setw(6) << std::setprecision(2)
516  << (element->fRelativeAbundanceVector[i])/perCent << " %";
517 
518  flux.precision(prec);
519  flux.setf(mode,std::ios::floatfield);
520  return flux;
521 }
522 
523 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
524 
525  std::ostream& operator<<(std::ostream& flux, G4Element& element)
526 {
527  flux << &element;
528  return flux;
529 }
530 
531 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
532 
533 std::ostream& operator<<(std::ostream& flux, G4ElementTable ElementTable)
534 {
535  //Dump info for all known elements
536  flux << "\n***** Table : Nb of elements = " << ElementTable.size()
537  << " *****\n" << G4endl;
538 
539  for (size_t i=0; i<ElementTable.size(); i++) flux << ElementTable[i]
540  << G4endl << G4endl;
541 
542  return flux;
543 }
544 
545 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4String symbol
Definition: TRTMaterials.hh:40
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:409
long G4long
Definition: G4Types.hh:80
const XML_Char * name
Definition: expat.h:151
tuple x
Definition: test.py:50
G4int GetNbOfShellElectrons(G4int index) const
Definition: G4Element.cc:380
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4int GetNistFirstIsotopeN(G4int Z) const
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
Float_t Z
Definition: plot.C:39
bool G4bool
Definition: G4Types.hh:79
static G4double GetBindingEnergy(G4int Z, G4int SubshellNb)
G4double iz
Definition: TRTMaterials.hh:39
void AddIsotope(G4Isotope *isotope, G4double RelativeAbundance)
Definition: G4Element.cc:151
G4Element(const G4String &name, const G4String &symbol, G4double Zeff, G4double Aeff)
Definition: G4Element.cc:74
G4double GetIsotopeAbundance(G4int Z, G4int N) const
G4int GetNumberOfNistIsotopes(G4int Z) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4int operator==(const G4Element &) const
Definition: G4Element.cc:485
G4int GetZ() const
Definition: G4Isotope.hh:91
int G4lrint(double ad)
Definition: templates.hh:163
G4int operator!=(const G4Element &) const
Definition: G4Element.cc:492
static G4int GetNumberOfElectrons(G4int Z, G4int SubshellNb)
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
virtual ~G4Element()
Definition: G4Element.cc:260
float perCent
Definition: hepunit.py:239
#define G4endl
Definition: G4ios.hh:61
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
G4double GetAtomicShell(G4int index) const
Definition: G4Element.cc:365
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
const G4String & GetName() const
Definition: G4Element.hh:127
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
float perMillion
Definition: hepunit.py:241
static G4int GetNumberOfShells(G4int Z)