source: trunk/source/materials/src/G4Element.cc @ 1058

Last change on this file since 1058 was 986, checked in by garnier, 15 years ago

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