Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 67044 2013-01-30 08:50:06Z 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 fNaturalAbandances (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 
74 G4Element::G4Element(const G4String& name, const G4String& symbol,
75  G4double zeff, G4double aeff)
76  : fName(name), fSymbol(symbol)
77 {
78  G4int iz = (G4int)zeff;
79  if (zeff<1.) {
80  G4ExceptionDescription ed;
81  ed << "Fail to create G4Element " << name
82  << " Z= " << zeff << " < 1 !" << G4endl;
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) << G4endl;
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) {
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  AddNaturalIsotopes();
137  } else {
138  theIsotopeVector = new G4IsotopeVector(n,0);
139  fRelativeAbundanceVector = new G4double[nIsotopes];
140  }
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 
145 // Add an isotope to the element
146 
147 void G4Element::AddIsotope(G4Isotope* isotope, G4double abundance)
148 {
149  if (theIsotopeVector == 0) {
150  G4ExceptionDescription ed;
151  ed << "Fail to add Isotope to G4Element " << fName
152  << " with Z= " << fZeff << " N= " << fNeff << G4endl;
153  G4Exception ("G4Element::AddIsotope()", "mat013", FatalException, ed);
154  return;
155  }
156  G4int iz = isotope->GetZ();
157 
158  // filling ...
159  if ( fNumberOfIsotopes < theIsotopeVector->size() ) {
160  // check same Z
161  if (fNumberOfIsotopes==0) { fZeff = G4double(iz); }
162  else if (G4double(iz) != fZeff) {
163  G4ExceptionDescription ed;
164  ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
165  << " with different Z= " << fZeff << fNeff
166  << G4endl;
167  G4Exception ("G4Element::AddIsotope()", "mat014", FatalException, ed);
168  return;
169  }
170  //Z ok
171  fRelativeAbundanceVector[fNumberOfIsotopes] = abundance;
172  (*theIsotopeVector)[fNumberOfIsotopes] = isotope;
173  ++fNumberOfIsotopes;
174  isotope->increaseCountUse();
175  } else {
176  G4ExceptionDescription ed;
177  ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
178  << " - more isotopes than declaired " << G4endl;
179  G4Exception ("G4Element::AddIsotope()", "mat015", FatalException, ed);
180  return;
181  }
182 
183  // filled.
184  if ( fNumberOfIsotopes == theIsotopeVector->size() ) {
185  // Compute Neff, Aeff
186  G4double wtSum=0.0;
187 
188  fNeff = fAeff = 0.0;
189  for (size_t i=0;i<fNumberOfIsotopes;i++) {
190  fAeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetA();
191  fNeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetN();
192  wtSum += fRelativeAbundanceVector[i];
193  }
194  fAeff /= wtSum;
195  fNeff /= wtSum;
196 
197  if(wtSum != 1.0) {
198  for(size_t 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 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
218 
219 void G4Element::InitializePointers()
220 {
221  theIsotopeVector = 0;
222  fRelativeAbundanceVector = 0;
223  fAtomicShells = 0;
224  fNbOfShellElectrons = 0;
225  fIonisation = 0;
226  fNumberOfIsotopes = 0;
227  fNaturalAbandances = false;
228 
229  // add initialisation of all remaining members
230  fZeff = 0;
231  fNeff = 0;
232  fAeff = 0;
233  fNbOfAtomicShells = 0;
234  fCountUse = 0;
235  fIndexZ = 0;
236  fIndexInTable = 0;
237  fNaturalAbandances = false;
238  fCoulomb = 0.0;
239  fRadTsai = 0.0;
240 }
241 
242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
243 
244 // Fake default constructor - sets only member data and allocates memory
245 // for usage restricted to object persistency
246 
247 G4Element::G4Element( __void__& )
248  : fZeff(0), fNeff(0), fAeff(0)
249 {
250  InitializePointers();
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
255 G4Element::~G4Element()
256 {
257  // G4cout << "### Destruction of element " << fName << " started" <<G4endl;
258 
259  if (theIsotopeVector) { delete theIsotopeVector; }
260  if (fRelativeAbundanceVector) { delete [] fRelativeAbundanceVector; }
261  if (fAtomicShells) { delete [] fAtomicShells; }
262  if (fNbOfShellElectrons) { delete [] fNbOfShellElectrons; }
263  if (fIonisation) { delete fIonisation; }
264 
265  //remove this element from theElementTable
266  theElementTable[fIndexInTable] = 0;
267 }
268 
269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
270 
271 void G4Element::ComputeDerivedQuantities()
272 {
273  // some basic functions of the atomic number
274 
275  // Store in table
276  theElementTable.push_back(this);
277  fIndexInTable = theElementTable.size() - 1;
278 
279  // check if elements with same Z already exist
280  fIndexZ = 0;
281  for (size_t J=0 ; J<fIndexInTable ; J++) {
282  if (theElementTable[J]->GetZ() == fZeff) { ++fIndexZ; }
283  }
284  //nb of materials which use this element
285  fCountUse = 0;
286 
287  // Radiation Length
288  ComputeCoulombFactor();
289  ComputeLradTsaiFactor();
290 
291  // parameters for energy loss by ionisation
292  if (fIonisation) { delete fIonisation; }
293  fIonisation = new G4IonisParamElm(fZeff);
294 }
295 
296 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297 
298 void G4Element::ComputeCoulombFactor()
299 {
300  //
301  // Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
302 
303  const G4double k1 = 0.0083 , k2 = 0.20206 ,k3 = 0.0020 , k4 = 0.0369 ;
304 
305  G4double az2 = (fine_structure_const*fZeff)*(fine_structure_const*fZeff);
306  G4double az4 = az2 * az2;
307 
308  fCoulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
309 }
310 
311 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312 
313 void G4Element::ComputeLradTsaiFactor()
314 {
315  //
316  // Compute Tsai's Expression for the Radiation Length
317  // (Phys Rev. D50 3-1 (1994) page 1254)
318 
319  const G4double Lrad_light[] = {5.31 , 4.79 , 4.74 , 4.71} ;
320  const G4double Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
321 
322  const G4double logZ3 = std::log(fZeff)/3.;
323 
324  G4double Lrad, Lprad;
325  G4int iz = (G4int)(fZeff+0.5) - 1 ;
326  if (iz <= 3) { Lrad = Lrad_light[iz] ; Lprad = Lprad_light[iz] ; }
327  else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}
328 
329  fRadTsai = 4*alpha_rcl2*fZeff*(fZeff*(Lrad-fCoulomb) + Lprad);
330 }
331 
332 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
333 
334 void G4Element::AddNaturalIsotopes()
335 {
336  G4int Z = G4lrint(fZeff);
337  G4NistManager* nist = G4NistManager::Instance();
338  G4int n = nist->GetNumberOfNistIsotopes(Z);
339  G4int N0 = nist->GetNistFirstIsotopeN(Z);
340  fNumberOfIsotopes = 0;
341  for(G4int i=0; i<n; ++i) {
342  if(nist->GetIsotopeAbundance(Z, N0+i) > 0.0) { ++fNumberOfIsotopes; }
343  }
344  theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
345  fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
346  G4int idx = 0;
347  G4double xsum = 0.0;
348  for(G4int i=0; i<n; ++i) {
349  G4int N = N0 + i;
350  G4double x = nist->GetIsotopeAbundance(Z, N);
351  if(x > 0.0) {
352  std::ostringstream strm;
353  strm << fSymbol << N;
354  (*theIsotopeVector)[idx] = new G4Isotope(strm.str(),Z, N, 0.0, 0);
355  fRelativeAbundanceVector[idx] = x;
356  xsum += x;
357  ++idx;
358  }
359  }
360  if(xsum != 0.0 && xsum != 1.0) {
361  for(G4int i=0; i<idx; ++i) { fRelativeAbundanceVector[i] /= xsum; }
362  }
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
366 
367 G4double G4Element::GetAtomicShell(G4int i) const
368 {
369  if (i<0 || i>=fNbOfAtomicShells) {
370  G4ExceptionDescription ed;
371  ed << "Invalid argument " << i << " in for G4Element " << fName
372  << " with Z= " << fZeff
373  << " and Nshells= " << fNbOfAtomicShells
374  << G4endl;
375  G4Exception("G4Element::GetAtomicShell()", "mat016", FatalException, ed);
376  return 0.0;
377  }
378  return fAtomicShells[i];
379 }
380 
381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
382 
383 G4int G4Element::GetNbOfShellElectrons(G4int i) const
384 {
385  if (i<0 || i>=fNbOfAtomicShells) {
386  G4ExceptionDescription ed;
387  ed << "Invalid argument " << i << " for G4Element " << fName
388  << " with Z= " << fZeff
389  << " and Nshells= " << fNbOfAtomicShells
390  << G4endl;
391  G4Exception("G4Element::GetNbOfShellElectrons()", "mat016", FatalException, ed);
392  return 0;
393  }
394  return fNbOfShellElectrons[i];
395 }
396 
397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
398 
399 const G4ElementTable* G4Element::GetElementTable()
400 {
401  return &theElementTable;
402 }
403 
404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
405 
406 size_t G4Element::GetNumberOfElements()
407 {
408  return theElementTable.size();
409 }
410 
411 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
412 
413 G4Element* G4Element::GetElement(G4String elementName, G4bool warning)
414 {
415  // search the element by its name
416  for (size_t J=0 ; J<theElementTable.size() ; J++)
417  {
418  if (theElementTable[J]->GetName() == elementName)
419  return theElementTable[J];
420  }
421 
422  // the element does not exist in the table
423  if (warning) {
424  G4cout << "\n---> warning from G4Element::GetElement(). The element: "
425  << elementName << " does not exist in the table. Return NULL pointer."
426  << G4endl;
427  }
428  return 0;
429 }
430 
431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432 
433 G4Element::G4Element(G4Element& right)
434 {
435  InitializePointers();
436  *this = right;
437 
438  // Store this new element in table and set the index
439  theElementTable.push_back(this);
440  fIndexInTable = theElementTable.size() - 1;
441 }
442 
443 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
444 
445 const G4Element& G4Element::operator=(const G4Element& right)
446 {
447  if (this != &right)
448  {
449  fName = right.fName;
450  fSymbol = right.fSymbol;
451  fZeff = right.fZeff;
452  fNeff = right.fNeff;
453  fAeff = right.fAeff;
454 
455  if (fAtomicShells) delete [] fAtomicShells;
456  fNbOfAtomicShells = right.fNbOfAtomicShells;
457  fAtomicShells = new G4double[fNbOfAtomicShells];
458 
459  if (fNbOfShellElectrons) delete [] fNbOfShellElectrons;
460  fNbOfAtomicShells = right.fNbOfAtomicShells;
461  fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
462 
463  for ( G4int i = 0; i < fNbOfAtomicShells; i++ )
464  {
465  fAtomicShells[i] = right.fAtomicShells[i];
466  fNbOfShellElectrons[i] = right.fNbOfShellElectrons[i];
467  }
468  if (theIsotopeVector) delete theIsotopeVector;
469  if (fRelativeAbundanceVector) delete [] fRelativeAbundanceVector;
470 
471  fNumberOfIsotopes = right.fNumberOfIsotopes;
472  if (fNumberOfIsotopes > 0)
473  {
474  theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
475  fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
476  for (size_t i=0;i<fNumberOfIsotopes;i++)
477  {
478  (*theIsotopeVector)[i] = (*right.theIsotopeVector)[i];
479  fRelativeAbundanceVector[i] = right.fRelativeAbundanceVector[i];
480  }
481  }
482  ComputeDerivedQuantities();
483  }
484  return *this;
485 }
486 
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
488 
489 G4int G4Element::operator==(const G4Element& right) const
490 {
491  return (this == (G4Element*) &right);
492 }
493 
494 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
495 
496 G4int G4Element::operator!=(const G4Element& right) const
497 {
498  return (this != (G4Element*) &right);
499 }
500 
501 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
502 
503 std::ostream& operator<<(std::ostream& flux, G4Element* element)
504 {
505  std::ios::fmtflags mode = flux.flags();
506  flux.setf(std::ios::fixed,std::ios::floatfield);
507  G4long prec = flux.precision(3);
508 
509  flux
510  << " Element: " << element->fName << " (" << element->fSymbol << ")"
511  << " Z = " << std::setw(4) << std::setprecision(1) << element->fZeff
512  << " N = " << std::setw(5) << std::setprecision(1) << element->fNeff
513  << " A = " << std::setw(6) << std::setprecision(2)
514  << (element->fAeff)/(g/mole) << " g/mole";
515 
516  for (size_t i=0; i<element->fNumberOfIsotopes; i++)
517  flux
518  << "\n ---> " << (*(element->theIsotopeVector))[i]
519  << " abundance: " << std::setw(6) << std::setprecision(2)
520  << (element->fRelativeAbundanceVector[i])/perCent << " %";
521 
522  flux.precision(prec);
523  flux.setf(mode,std::ios::floatfield);
524  return flux;
525 }
526 
527 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
528 
529  std::ostream& operator<<(std::ostream& flux, G4Element& element)
530 {
531  flux << &element;
532  return flux;
533 }
534 
535 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
536 
537 std::ostream& operator<<(std::ostream& flux, G4ElementTable ElementTable)
538 {
539  //Dump info for all known elements
540  flux << "\n***** Table : Nb of elements = " << ElementTable.size()
541  << " *****\n" << G4endl;
542 
543  for (size_t i=0; i<ElementTable.size(); i++) flux << ElementTable[i]
544  << G4endl << G4endl;
545 
546  return flux;
547 }
548 
549 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Generated on Sat May 25 2013 14:33:20 for Geant4 by   doxygen 1.8.4