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

Last change on this file since 1355 was 1337, checked in by garnier, 15 years ago

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