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

Last change on this file since 822 was 822, checked in by garnier, 17 years ago

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