00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 #include <iomanip>
00060 #include <sstream>
00061
00062 #include "G4Element.hh"
00063 #include "G4AtomicShells.hh"
00064 #include "G4NistManager.hh"
00065 #include "G4PhysicalConstants.hh"
00066 #include "G4SystemOfUnits.hh"
00067
00068 G4ElementTable G4Element::theElementTable;
00069
00070
00071
00072
00073
00074 G4Element::G4Element(const G4String& name, const G4String& symbol,
00075 G4double zeff, G4double aeff)
00076 : fName(name), fSymbol(symbol)
00077 {
00078 G4int iz = (G4int)zeff;
00079 if (zeff<1.) {
00080 G4ExceptionDescription ed;
00081 ed << "Fail to create G4Element " << name
00082 << " Z= " << zeff << " < 1 !" << G4endl;
00083 G4Exception ("G4Element::G4Element()", "mat011", FatalException, ed);
00084 }
00085 if (std::abs(zeff - iz) > perMillion) {
00086 G4ExceptionDescription ed;
00087 ed << "G4Element Warning: " << name << " Z= " << zeff
00088 << " A= " << aeff/(g/mole) << G4endl;
00089 G4Exception("G4Element::G4Element()", "mat017", JustWarning, ed);
00090 }
00091
00092 InitializePointers();
00093
00094 fZeff = zeff;
00095 fNeff = aeff/(g/mole);
00096 fAeff = aeff;
00097
00098 if(fNeff < 1.0) fNeff = 1.0;
00099
00100 if (fNeff < zeff) {
00101 G4ExceptionDescription ed;
00102 ed << "Fail to create G4Element " << name
00103 << " with Z= " << zeff << " N= " << fNeff
00104 << " N < Z is not allowed" << G4endl;
00105 G4Exception ("G4Element::G4Element()", "mat012", FatalException, ed);
00106 }
00107
00108 fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
00109 fAtomicShells = new G4double[fNbOfAtomicShells];
00110 fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
00111
00112 AddNaturalIsotopes();
00113
00114 for (G4int i=0;i<fNbOfAtomicShells;i++)
00115 {
00116 fAtomicShells[i] = G4AtomicShells::GetBindingEnergy(iz, i);
00117 fNbOfShellElectrons[i] = G4AtomicShells::GetNumberOfElectrons(iz, i);
00118 }
00119 ComputeDerivedQuantities();
00120 }
00121
00122
00123
00124
00125
00126
00127 G4Element::G4Element(const G4String& name,
00128 const G4String& symbol, G4int nIsotopes)
00129 : fName(name),fSymbol(symbol)
00130 {
00131 InitializePointers();
00132
00133 size_t n = size_t(nIsotopes);
00134
00135 if(0 == nIsotopes) {
00136 AddNaturalIsotopes();
00137 } else {
00138 theIsotopeVector = new G4IsotopeVector(n,0);
00139 fRelativeAbundanceVector = new G4double[nIsotopes];
00140 }
00141 }
00142
00143
00144
00145
00146
00147 void G4Element::AddIsotope(G4Isotope* isotope, G4double abundance)
00148 {
00149 if (theIsotopeVector == 0) {
00150 G4ExceptionDescription ed;
00151 ed << "Fail to add Isotope to G4Element " << fName
00152 << " with Z= " << fZeff << " N= " << fNeff << G4endl;
00153 G4Exception ("G4Element::AddIsotope()", "mat013", FatalException, ed);
00154 return;
00155 }
00156 G4int iz = isotope->GetZ();
00157
00158
00159 if ( fNumberOfIsotopes < theIsotopeVector->size() ) {
00160
00161 if (fNumberOfIsotopes==0) { fZeff = G4double(iz); }
00162 else if (G4double(iz) != fZeff) {
00163 G4ExceptionDescription ed;
00164 ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
00165 << " with different Z= " << fZeff << fNeff
00166 << G4endl;
00167 G4Exception ("G4Element::AddIsotope()", "mat014", FatalException, ed);
00168 return;
00169 }
00170
00171 fRelativeAbundanceVector[fNumberOfIsotopes] = abundance;
00172 (*theIsotopeVector)[fNumberOfIsotopes] = isotope;
00173 ++fNumberOfIsotopes;
00174 isotope->increaseCountUse();
00175 } else {
00176 G4ExceptionDescription ed;
00177 ed << "Fail to add Isotope Z= " << iz << " to G4Element " << fName
00178 << " - more isotopes than declaired " << G4endl;
00179 G4Exception ("G4Element::AddIsotope()", "mat015", FatalException, ed);
00180 return;
00181 }
00182
00183
00184 if ( fNumberOfIsotopes == theIsotopeVector->size() ) {
00185
00186 G4double wtSum=0.0;
00187
00188 fNeff = fAeff = 0.0;
00189 for (size_t i=0;i<fNumberOfIsotopes;i++) {
00190 fAeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetA();
00191 fNeff += fRelativeAbundanceVector[i]*(*theIsotopeVector)[i]->GetN();
00192 wtSum += fRelativeAbundanceVector[i];
00193 }
00194 fAeff /= wtSum;
00195 fNeff /= wtSum;
00196
00197 if(wtSum != 1.0) {
00198 for(size_t i=0; i<fNumberOfIsotopes; ++i) {
00199 fRelativeAbundanceVector[i] /= wtSum;
00200 }
00201 }
00202
00203 fNbOfAtomicShells = G4AtomicShells::GetNumberOfShells(iz);
00204 fAtomicShells = new G4double[fNbOfAtomicShells];
00205 fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
00206
00207 for ( G4int j = 0; j < fNbOfAtomicShells; j++ )
00208 {
00209 fAtomicShells[j] = G4AtomicShells::GetBindingEnergy(iz, j);
00210 fNbOfShellElectrons[j] = G4AtomicShells::GetNumberOfElectrons(iz, j);
00211 }
00212 ComputeDerivedQuantities();
00213
00214 }
00215 }
00216
00217
00218
00219 void G4Element::InitializePointers()
00220 {
00221 theIsotopeVector = 0;
00222 fRelativeAbundanceVector = 0;
00223 fAtomicShells = 0;
00224 fNbOfShellElectrons = 0;
00225 fIonisation = 0;
00226 fNumberOfIsotopes = 0;
00227 fNaturalAbandances = false;
00228
00229
00230 fZeff = 0;
00231 fNeff = 0;
00232 fAeff = 0;
00233 fNbOfAtomicShells = 0;
00234 fCountUse = 0;
00235 fIndexZ = 0;
00236 fIndexInTable = 0;
00237 fNaturalAbandances = false;
00238 fCoulomb = 0.0;
00239 fRadTsai = 0.0;
00240 }
00241
00242
00243
00244
00245
00246
00247 G4Element::G4Element( __void__& )
00248 : fZeff(0), fNeff(0), fAeff(0)
00249 {
00250 InitializePointers();
00251 }
00252
00253
00254
00255 G4Element::~G4Element()
00256 {
00257
00258
00259 if (theIsotopeVector) { delete theIsotopeVector; }
00260 if (fRelativeAbundanceVector) { delete [] fRelativeAbundanceVector; }
00261 if (fAtomicShells) { delete [] fAtomicShells; }
00262 if (fNbOfShellElectrons) { delete [] fNbOfShellElectrons; }
00263 if (fIonisation) { delete fIonisation; }
00264
00265
00266 theElementTable[fIndexInTable] = 0;
00267 }
00268
00269
00270
00271 void G4Element::ComputeDerivedQuantities()
00272 {
00273
00274
00275
00276 theElementTable.push_back(this);
00277 fIndexInTable = theElementTable.size() - 1;
00278
00279
00280 fIndexZ = 0;
00281 for (size_t J=0 ; J<fIndexInTable ; J++) {
00282 if (theElementTable[J]->GetZ() == fZeff) { ++fIndexZ; }
00283 }
00284
00285 fCountUse = 0;
00286
00287
00288 ComputeCoulombFactor();
00289 ComputeLradTsaiFactor();
00290
00291
00292 if (fIonisation) { delete fIonisation; }
00293 fIonisation = new G4IonisParamElm(fZeff);
00294 }
00295
00296
00297
00298 void G4Element::ComputeCoulombFactor()
00299 {
00300
00301
00302
00303 const G4double k1 = 0.0083 , k2 = 0.20206 ,k3 = 0.0020 , k4 = 0.0369 ;
00304
00305 G4double az2 = (fine_structure_const*fZeff)*(fine_structure_const*fZeff);
00306 G4double az4 = az2 * az2;
00307
00308 fCoulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
00309 }
00310
00311
00312
00313 void G4Element::ComputeLradTsaiFactor()
00314 {
00315
00316
00317
00318
00319 const G4double Lrad_light[] = {5.31 , 4.79 , 4.74 , 4.71} ;
00320 const G4double Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
00321
00322 const G4double logZ3 = std::log(fZeff)/3.;
00323
00324 G4double Lrad, Lprad;
00325 G4int iz = (G4int)(fZeff+0.5) - 1 ;
00326 if (iz <= 3) { Lrad = Lrad_light[iz] ; Lprad = Lprad_light[iz] ; }
00327 else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}
00328
00329 fRadTsai = 4*alpha_rcl2*fZeff*(fZeff*(Lrad-fCoulomb) + Lprad);
00330 }
00331
00332
00333
00334 void G4Element::AddNaturalIsotopes()
00335 {
00336 G4int Z = G4lrint(fZeff);
00337 G4NistManager* nist = G4NistManager::Instance();
00338 G4int n = nist->GetNumberOfNistIsotopes(Z);
00339 G4int N0 = nist->GetNistFirstIsotopeN(Z);
00340 fNumberOfIsotopes = 0;
00341 for(G4int i=0; i<n; ++i) {
00342 if(nist->GetIsotopeAbundance(Z, N0+i) > 0.0) { ++fNumberOfIsotopes; }
00343 }
00344 theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
00345 fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
00346 G4int idx = 0;
00347 G4double xsum = 0.0;
00348 for(G4int i=0; i<n; ++i) {
00349 G4int N = N0 + i;
00350 G4double x = nist->GetIsotopeAbundance(Z, N);
00351 if(x > 0.0) {
00352 std::ostringstream strm;
00353 strm << fSymbol << N;
00354 (*theIsotopeVector)[idx] = new G4Isotope(strm.str(),Z, N, 0.0, 0);
00355 fRelativeAbundanceVector[idx] = x;
00356 xsum += x;
00357 ++idx;
00358 }
00359 }
00360 if(xsum != 0.0 && xsum != 1.0) {
00361 for(G4int i=0; i<idx; ++i) { fRelativeAbundanceVector[i] /= xsum; }
00362 }
00363 }
00364
00365
00366
00367 G4double G4Element::GetAtomicShell(G4int i) const
00368 {
00369 if (i<0 || i>=fNbOfAtomicShells) {
00370 G4ExceptionDescription ed;
00371 ed << "Invalid argument " << i << " in for G4Element " << fName
00372 << " with Z= " << fZeff
00373 << " and Nshells= " << fNbOfAtomicShells
00374 << G4endl;
00375 G4Exception("G4Element::GetAtomicShell()", "mat016", FatalException, ed);
00376 return 0.0;
00377 }
00378 return fAtomicShells[i];
00379 }
00380
00381
00382
00383 G4int G4Element::GetNbOfShellElectrons(G4int i) const
00384 {
00385 if (i<0 || i>=fNbOfAtomicShells) {
00386 G4ExceptionDescription ed;
00387 ed << "Invalid argument " << i << " for G4Element " << fName
00388 << " with Z= " << fZeff
00389 << " and Nshells= " << fNbOfAtomicShells
00390 << G4endl;
00391 G4Exception("G4Element::GetNbOfShellElectrons()", "mat016", FatalException, ed);
00392 return 0;
00393 }
00394 return fNbOfShellElectrons[i];
00395 }
00396
00397
00398
00399 const G4ElementTable* G4Element::GetElementTable()
00400 {
00401 return &theElementTable;
00402 }
00403
00404
00405
00406 size_t G4Element::GetNumberOfElements()
00407 {
00408 return theElementTable.size();
00409 }
00410
00411
00412
00413 G4Element* G4Element::GetElement(G4String elementName, G4bool warning)
00414 {
00415
00416 for (size_t J=0 ; J<theElementTable.size() ; J++)
00417 {
00418 if (theElementTable[J]->GetName() == elementName)
00419 return theElementTable[J];
00420 }
00421
00422
00423 if (warning) {
00424 G4cout << "\n---> warning from G4Element::GetElement(). The element: "
00425 << elementName << " does not exist in the table. Return NULL pointer."
00426 << G4endl;
00427 }
00428 return 0;
00429 }
00430
00431
00432
00433 G4Element::G4Element(G4Element& right)
00434 {
00435 InitializePointers();
00436 *this = right;
00437
00438
00439 theElementTable.push_back(this);
00440 fIndexInTable = theElementTable.size() - 1;
00441 }
00442
00443
00444
00445 const G4Element& G4Element::operator=(const G4Element& right)
00446 {
00447 if (this != &right)
00448 {
00449 fName = right.fName;
00450 fSymbol = right.fSymbol;
00451 fZeff = right.fZeff;
00452 fNeff = right.fNeff;
00453 fAeff = right.fAeff;
00454
00455 if (fAtomicShells) delete [] fAtomicShells;
00456 fNbOfAtomicShells = right.fNbOfAtomicShells;
00457 fAtomicShells = new G4double[fNbOfAtomicShells];
00458
00459 if (fNbOfShellElectrons) delete [] fNbOfShellElectrons;
00460 fNbOfAtomicShells = right.fNbOfAtomicShells;
00461 fNbOfShellElectrons = new G4int[fNbOfAtomicShells];
00462
00463 for ( G4int i = 0; i < fNbOfAtomicShells; i++ )
00464 {
00465 fAtomicShells[i] = right.fAtomicShells[i];
00466 fNbOfShellElectrons[i] = right.fNbOfShellElectrons[i];
00467 }
00468 if (theIsotopeVector) delete theIsotopeVector;
00469 if (fRelativeAbundanceVector) delete [] fRelativeAbundanceVector;
00470
00471 fNumberOfIsotopes = right.fNumberOfIsotopes;
00472 if (fNumberOfIsotopes > 0)
00473 {
00474 theIsotopeVector = new G4IsotopeVector(fNumberOfIsotopes,0);
00475 fRelativeAbundanceVector = new G4double[fNumberOfIsotopes];
00476 for (size_t i=0;i<fNumberOfIsotopes;i++)
00477 {
00478 (*theIsotopeVector)[i] = (*right.theIsotopeVector)[i];
00479 fRelativeAbundanceVector[i] = right.fRelativeAbundanceVector[i];
00480 }
00481 }
00482 ComputeDerivedQuantities();
00483 }
00484 return *this;
00485 }
00486
00487
00488
00489 G4int G4Element::operator==(const G4Element& right) const
00490 {
00491 return (this == (G4Element*) &right);
00492 }
00493
00494
00495
00496 G4int G4Element::operator!=(const G4Element& right) const
00497 {
00498 return (this != (G4Element*) &right);
00499 }
00500
00501
00502
00503 std::ostream& operator<<(std::ostream& flux, G4Element* element)
00504 {
00505 std::ios::fmtflags mode = flux.flags();
00506 flux.setf(std::ios::fixed,std::ios::floatfield);
00507 G4long prec = flux.precision(3);
00508
00509 flux
00510 << " Element: " << element->fName << " (" << element->fSymbol << ")"
00511 << " Z = " << std::setw(4) << std::setprecision(1) << element->fZeff
00512 << " N = " << std::setw(5) << std::setprecision(1) << element->fNeff
00513 << " A = " << std::setw(6) << std::setprecision(2)
00514 << (element->fAeff)/(g/mole) << " g/mole";
00515
00516 for (size_t i=0; i<element->fNumberOfIsotopes; i++)
00517 flux
00518 << "\n ---> " << (*(element->theIsotopeVector))[i]
00519 << " abundance: " << std::setw(6) << std::setprecision(2)
00520 << (element->fRelativeAbundanceVector[i])/perCent << " %";
00521
00522 flux.precision(prec);
00523 flux.setf(mode,std::ios::floatfield);
00524 return flux;
00525 }
00526
00527
00528
00529 std::ostream& operator<<(std::ostream& flux, G4Element& element)
00530 {
00531 flux << &element;
00532 return flux;
00533 }
00534
00535
00536
00537 std::ostream& operator<<(std::ostream& flux, G4ElementTable ElementTable)
00538 {
00539
00540 flux << "\n***** Table : Nb of elements = " << ElementTable.size()
00541 << " *****\n" << G4endl;
00542
00543 for (size_t i=0; i<ElementTable.size(); i++) flux << ElementTable[i]
00544 << G4endl << G4endl;
00545
00546 return flux;
00547 }
00548
00549