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 #include "G4PenelopeOscillatorManager.hh"
00051
00052 #include "globals.hh"
00053 #include "G4PhysicalConstants.hh"
00054 #include "G4SystemOfUnits.hh"
00055 #include "G4AtomicTransitionManager.hh"
00056 #include "G4AtomicShell.hh"
00057 #include "G4Material.hh"
00058
00059
00060
00061 G4PenelopeOscillatorManager::G4PenelopeOscillatorManager() :
00062 oscillatorStoreIonisation(0),oscillatorStoreCompton(0),atomicNumber(0),
00063 atomicMass(0),excitationEnergy(0),plasmaSquared(0),atomsPerMolecule(0),
00064 atomTablePerMolecule(0)
00065 {
00066 fReadElementData = false;
00067 for (G4int i=0;i<5;i++)
00068 {
00069 for (G4int j=0;j<2000;j++)
00070 elementData[i][j] = 0.;
00071 }
00072 verbosityLevel = 0;
00073 }
00074
00075
00076
00077 G4PenelopeOscillatorManager::~G4PenelopeOscillatorManager()
00078 {
00079 Clear();
00080 delete instance;
00081 }
00082
00083
00084
00085 G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::instance = 0;
00086
00087
00088
00089 G4PenelopeOscillatorManager* G4PenelopeOscillatorManager::GetOscillatorManager()
00090 {
00091 if (!instance)
00092 instance = new G4PenelopeOscillatorManager();
00093 return instance;
00094 }
00095
00096
00097
00098 void G4PenelopeOscillatorManager::Clear()
00099 {
00100 if (verbosityLevel > 1)
00101 G4cout << " G4PenelopeOscillatorManager::Clear() - Clean Oscillator Tables" << G4endl;
00102
00103
00104 std::map<const G4Material*,G4PenelopeOscillatorTable*>::iterator i;
00105 for (i=oscillatorStoreIonisation->begin();i != oscillatorStoreIonisation->end();i++)
00106 {
00107 G4PenelopeOscillatorTable* table = i->second;
00108 if (table)
00109 {
00110 for (size_t k=0;k<table->size();k++)
00111 {
00112 if ((*table)[k])
00113 delete ((*table)[k]);
00114 }
00115 delete table;
00116 }
00117 }
00118 delete oscillatorStoreIonisation;
00119
00120
00121 for (i=oscillatorStoreCompton->begin();i != oscillatorStoreCompton->end();i++)
00122 {
00123 G4PenelopeOscillatorTable* table = i->second;
00124 if (table)
00125 {
00126 for (size_t k=0;k<table->size();k++)
00127 {
00128 if ((*table)[k])
00129 delete ((*table)[k]);
00130 }
00131 delete table;
00132 }
00133 }
00134 delete oscillatorStoreCompton;
00135
00136 if (atomicMass) delete atomicMass;
00137 if (atomicNumber) delete atomicNumber;
00138 if (excitationEnergy) delete excitationEnergy;
00139 if (plasmaSquared) delete plasmaSquared;
00140 if (atomsPerMolecule) delete atomsPerMolecule;
00141 if (atomTablePerMolecule) delete atomTablePerMolecule;
00142 }
00143
00144
00145
00146 void G4PenelopeOscillatorManager::Dump(const G4Material* material)
00147 {
00148 G4PenelopeOscillatorTable* theTable = GetOscillatorTableIonisation(material);
00149 if (!theTable)
00150 {
00151 G4cout << " G4PenelopeOscillatorManager::Dump " << G4endl;
00152 G4cout << "Problem in retrieving the Ionisation Oscillator Table for " << material->GetName() << G4endl;
00153 return;
00154 }
00155 G4cout << "*********************************************************************" << G4endl;
00156 G4cout << " Penelope Oscillator Table Ionisation for " << material->GetName() << G4endl;
00157 G4cout << "*********************************************************************" << G4endl;
00158 G4cout << "The table contains " << theTable->size() << " oscillators " << G4endl;
00159 G4cout << "*********************************************************************" << G4endl;
00160 if (theTable->size() < 10)
00161 for (size_t k=0;k<theTable->size();k++)
00162 {
00163 G4cout << "Oscillator # " << k << " Z = " << (*theTable)[k]->GetParentZ() <<
00164 " Shell Flag = " << (*theTable)[k]->GetShellFlag() <<
00165 " Parent shell ID = " << (*theTable)[k]->GetParentShellID() << G4endl;
00166 G4cout << "Ionisation energy = " << (*theTable)[k]->GetIonisationEnergy()/eV << " eV" << G4endl;
00167 G4cout << "Occupation number = " << (*theTable)[k]->GetOscillatorStrength() << G4endl;
00168 G4cout << "Resonance energy = " << (*theTable)[k]->GetResonanceEnergy()/eV << " eV" << G4endl;
00169 G4cout << "Cufoff resonance energy = " <<
00170 (*theTable)[k]->GetCutoffRecoilResonantEnergy()/eV << " eV" << G4endl;
00171 G4cout << "*********************************************************************" << G4endl;
00172 }
00173 for (size_t k=0;k<theTable->size();k++)
00174 {
00175 G4cout << k << " " << (*theTable)[k]->GetOscillatorStrength() << " " <<
00176 (*theTable)[k]->GetIonisationEnergy()/eV << " " << (*theTable)[k]->GetResonanceEnergy()/eV << " " <<
00177 (*theTable)[k]->GetParentZ() << " " << (*theTable)[k]->GetShellFlag() << " " <<
00178 (*theTable)[k]->GetParentShellID() << G4endl;
00179 }
00180 G4cout << "*********************************************************************" << G4endl;
00181
00182
00183
00184 theTable = GetOscillatorTableCompton(material);
00185 if (!theTable)
00186 {
00187 G4cout << " G4PenelopeOscillatorManager::Dump " << G4endl;
00188 G4cout << "Problem in retrieving the Compton Oscillator Table for " << material->GetName() << G4endl;
00189 return;
00190 }
00191 G4cout << "*********************************************************************" << G4endl;
00192 G4cout << " Penelope Oscillator Table Compton for " << material->GetName() << G4endl;
00193 G4cout << "*********************************************************************" << G4endl;
00194 G4cout << "The table contains " << theTable->size() << " oscillators " << G4endl;
00195 G4cout << "*********************************************************************" << G4endl;
00196 if (theTable->size() < 10)
00197 for (size_t k=0;k<theTable->size();k++)
00198 {
00199 G4cout << "Oscillator # " << k << " Z = " << (*theTable)[k]->GetParentZ() <<
00200 " Shell Flag = " << (*theTable)[k]->GetShellFlag() <<
00201 " Parent shell ID = " << (*theTable)[k]->GetParentShellID() << G4endl;
00202 G4cout << "Compton index = " << (*theTable)[k]->GetHartreeFactor() << G4endl;
00203 G4cout << "Ionisation energy = " << (*theTable)[k]->GetIonisationEnergy()/eV << " eV" << G4endl;
00204 G4cout << "Occupation number = " << (*theTable)[k]->GetOscillatorStrength() << G4endl;
00205 G4cout << "*********************************************************************" << G4endl;
00206 }
00207 for (size_t k=0;k<theTable->size();k++)
00208 {
00209 G4cout << k << " " << (*theTable)[k]->GetOscillatorStrength() << " " <<
00210 (*theTable)[k]->GetIonisationEnergy()/eV << " " << (*theTable)[k]->GetHartreeFactor() << " " <<
00211 (*theTable)[k]->GetParentZ() << " " << (*theTable)[k]->GetShellFlag() << " " <<
00212 (*theTable)[k]->GetParentShellID() << G4endl;
00213 }
00214 G4cout << "*********************************************************************" << G4endl;
00215
00216 return;
00217 }
00218
00219
00220
00221 void G4PenelopeOscillatorManager::CheckForTablesCreated()
00222 {
00223
00224
00225 if (!oscillatorStoreIonisation)
00226 {
00227 oscillatorStoreIonisation = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
00228 if (!fReadElementData)
00229 ReadElementData();
00230 if (!oscillatorStoreIonisation)
00231
00232 G4Exception("G4PenelopeOscillatorManager::GetOscillatorTableIonisation()",
00233 "em2034",FatalException,
00234 "Problem in allocating the Oscillator Store for Ionisation");
00235 }
00236
00237 if (!oscillatorStoreCompton)
00238 {
00239 oscillatorStoreCompton = new std::map<const G4Material*,G4PenelopeOscillatorTable*>;
00240 if (!fReadElementData)
00241 ReadElementData();
00242 if (!oscillatorStoreCompton)
00243
00244 G4Exception("G4PenelopeOscillatorManager::GetOscillatorTableIonisation()",
00245 "em2034",FatalException,
00246 "Problem in allocating the Oscillator Store for Compton");
00247 }
00248
00249 if (!atomicNumber)
00250 atomicNumber = new std::map<const G4Material*,G4double>;
00251 if (!atomicMass)
00252 atomicMass = new std::map<const G4Material*,G4double>;
00253 if (!excitationEnergy)
00254 excitationEnergy = new std::map<const G4Material*,G4double>;
00255 if (!plasmaSquared)
00256 plasmaSquared = new std::map<const G4Material*,G4double>;
00257 if (!atomsPerMolecule)
00258 atomsPerMolecule = new std::map<const G4Material*,G4double>;
00259 if (!atomTablePerMolecule)
00260 atomTablePerMolecule = new std::map< std::pair<const G4Material*,G4int>, G4double>;
00261 }
00262
00263
00264
00265
00266 G4double G4PenelopeOscillatorManager::GetTotalZ(const G4Material* mat)
00267 {
00268
00269 CheckForTablesCreated();
00270
00271
00272 if (atomicNumber->count(mat))
00273 return atomicNumber->find(mat)->second;
00274
00275
00276 BuildOscillatorTable(mat);
00277
00278
00279 if (atomicNumber->count(mat))
00280 return atomicNumber->find(mat)->second;
00281 else
00282 {
00283 G4cout << "G4PenelopeOscillatorManager::GetTotalZ() " << G4endl;
00284 G4cout << "Impossible to retrieve the total Z for " << mat->GetName() << G4endl;
00285 return 0;
00286 }
00287 }
00288
00289
00290
00291 G4double G4PenelopeOscillatorManager::GetTotalA(const G4Material* mat)
00292 {
00293
00294 CheckForTablesCreated();
00295
00296
00297 if (atomicMass->count(mat))
00298 return atomicMass->find(mat)->second;
00299
00300
00301 BuildOscillatorTable(mat);
00302
00303
00304 if (atomicMass->count(mat))
00305 return atomicMass->find(mat)->second;
00306 else
00307 {
00308 G4cout << "G4PenelopeOscillatorManager::GetTotalA() " << G4endl;
00309 G4cout << "Impossible to retrieve the total A for " << mat->GetName() << G4endl;
00310 return 0;
00311 }
00312 }
00313
00314
00315
00316 G4PenelopeOscillatorTable* G4PenelopeOscillatorManager::GetOscillatorTableIonisation(const G4Material* mat)
00317 {
00318
00319 CheckForTablesCreated();
00320
00321
00322 if (oscillatorStoreIonisation->count(mat))
00323 {
00324
00325 return oscillatorStoreIonisation->find(mat)->second;
00326 }
00327
00328
00329 BuildOscillatorTable(mat);
00330
00331
00332 if (oscillatorStoreIonisation->count(mat))
00333 return oscillatorStoreIonisation->find(mat)->second;
00334 else
00335 {
00336 G4cout << "G4PenelopeOscillatorManager::GetOscillatorTableIonisation() " << G4endl;
00337 G4cout << "Impossible to create ionisation oscillator table for " << mat->GetName() << G4endl;
00338 return NULL;
00339 }
00340 }
00341
00342
00343
00344 G4PenelopeOscillator* G4PenelopeOscillatorManager::GetOscillatorIonisation(const G4Material* material,
00345 G4int index)
00346 {
00347 G4PenelopeOscillatorTable* theTable = GetOscillatorTableIonisation(material);
00348 if (((size_t)index) < theTable->size())
00349 return (*theTable)[index];
00350 else
00351 {
00352 G4cout << "WARNING: Ionisation table for material " << material->GetName() << " has " <<
00353 theTable->size() << " oscillators" << G4endl;
00354 G4cout << "Oscillator #" << index << " cannot be retrieved" << G4endl;
00355 G4cout << "Returning null pointer" << G4endl;
00356 return NULL;
00357 }
00358 }
00359
00360
00361
00362
00363 G4PenelopeOscillatorTable* G4PenelopeOscillatorManager::GetOscillatorTableCompton(const G4Material* mat)
00364 {
00365
00366 CheckForTablesCreated();
00367
00368
00369 if (oscillatorStoreCompton->count(mat))
00370 {
00371
00372 return oscillatorStoreCompton->find(mat)->second;
00373 }
00374
00375
00376 BuildOscillatorTable(mat);
00377
00378
00379 if (oscillatorStoreCompton->count(mat))
00380 return oscillatorStoreCompton->find(mat)->second;
00381 else
00382 {
00383 G4cout << "G4PenelopeOscillatorManager::GetOscillatorTableCompton() " << G4endl;
00384 G4cout << "Impossible to create Compton oscillator table for " << mat->GetName() << G4endl;
00385 return NULL;
00386 }
00387 }
00388
00389
00390
00391 G4PenelopeOscillator* G4PenelopeOscillatorManager::GetOscillatorCompton(const G4Material* material,
00392 G4int index)
00393 {
00394 G4PenelopeOscillatorTable* theTable = GetOscillatorTableCompton(material);
00395 if (((size_t)index) < theTable->size())
00396 return (*theTable)[index];
00397 else
00398 {
00399 G4cout << "WARNING: Compton table for material " << material->GetName() << " has " <<
00400 theTable->size() << " oscillators" << G4endl;
00401 G4cout << "Oscillator #" << index << " cannot be retrieved" << G4endl;
00402 G4cout << "Returning null pointer" << G4endl;
00403 return NULL;
00404 }
00405 }
00406
00407
00408
00409 void G4PenelopeOscillatorManager::BuildOscillatorTable(const G4Material* material)
00410 {
00411
00412
00413 G4double meanAtomExcitationEnergy[99] = {19.2*eV, 41.8*eV, 40.0*eV, 63.7*eV, 76.0*eV, 81.0*eV,
00414 82.0*eV, 95.0*eV,115.0*eV,137.0*eV,149.0*eV,156.0*eV,
00415 166.0*eV,
00416 173.0*eV,173.0*eV,180.0*eV,174.0*eV,188.0*eV,190.0*eV,191.0*eV,
00417 216.0*eV,233.0*eV,245.0*eV,257.0*eV,272.0*eV,286.0*eV,297.0*eV,
00418 311.0*eV,322.0*eV,330.0*eV,334.0*eV,350.0*eV,347.0*eV,348.0*eV,
00419 343.0*eV,352.0*eV,363.0*eV,366.0*eV,379.0*eV,393.0*eV,417.0*eV,
00420 424.0*eV,428.0*eV,441.0*eV,449.0*eV,470.0*eV,470.0*eV,469.0*eV,
00421 488.0*eV,488.0*eV,487.0*eV,485.0*eV,491.0*eV,482.0*eV,488.0*eV,
00422 491.0*eV,501.0*eV,523.0*eV,535.0*eV,546.0*eV,560.0*eV,574.0*eV,
00423 580.0*eV,591.0*eV,614.0*eV,628.0*eV,650.0*eV,658.0*eV,674.0*eV,
00424 684.0*eV,694.0*eV,705.0*eV,718.0*eV,727.0*eV,736.0*eV,746.0*eV,
00425 757.0*eV,790.0*eV,790.0*eV,800.0*eV,810.0*eV,823.0*eV,823.0*eV,
00426 830.0*eV,825.0*eV,794.0*eV,827.0*eV,826.0*eV,841.0*eV,847.0*eV,
00427 878.0*eV,890.0*eV,902.0*eV,921.0*eV,934.0*eV,939.0*eV,952.0*eV,
00428 966.0*eV,980.0*eV};
00429
00430 if (verbosityLevel > 0)
00431 G4cout << "Going to build Oscillator Table for " << material->GetName() << G4endl;
00432
00433 G4int nElements = material->GetNumberOfElements();
00434 const G4ElementVector* elementVector = material->GetElementVector();
00435
00436
00437
00438
00439 const G4double* fractionVector = material->GetFractionVector();
00440
00441
00442
00443
00444 G4double totalZ = 0;
00445 G4double totalMolecularWeight = 0;
00446 G4double meanExcitationEnergy = 0;
00447
00448 std::vector<G4double> *StechiometricFactors = new std::vector<G4double>;
00449
00450 for (G4int i=0;i<nElements;i++)
00451 {
00452
00453 G4double fraction = fractionVector[i];
00454 G4double atomicWeigth = (*elementVector)[i]->GetAtomicMassAmu();
00455 StechiometricFactors->push_back(fraction/atomicWeigth);
00456 }
00457
00458 G4double MaxStechiometricFactor = 0.;
00459 for (G4int i=0;i<nElements;i++)
00460 {
00461 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
00462 MaxStechiometricFactor = (*StechiometricFactors)[i];
00463 }
00464 if (MaxStechiometricFactor<1e-16)
00465 {
00466 G4ExceptionDescription ed;
00467 ed << "Problem with the mass composition of " << material->GetName() << G4endl;
00468 ed << "MaxStechiometricFactor = " << MaxStechiometricFactor << G4endl;
00469 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
00470 "em2035",FatalException,ed);
00471 }
00472
00473 for (G4int i=0;i<nElements;i++)
00474 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
00475
00476
00477 G4double theatomsPerMolecule = 0;
00478 for (G4int i=0;i<nElements;i++)
00479 theatomsPerMolecule += (*StechiometricFactors)[i];
00480 G4double moleculeDensity =
00481 material->GetTotNbOfAtomsPerVolume()/theatomsPerMolecule;
00482
00483
00484 if (verbosityLevel > 1)
00485 {
00486 for (size_t i=0;i<StechiometricFactors->size();i++)
00487 {
00488 G4cout << "Element " << (*elementVector)[i]->GetSymbol() << " (Z = " <<
00489 (*elementVector)[i]->GetZ() << ") --> " <<
00490 (*StechiometricFactors)[i] << " atoms/molecule " << G4endl;
00491 }
00492 }
00493
00494
00495 for (G4int i=0;i<nElements;i++)
00496 {
00497 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
00498 totalZ += iZ * (*StechiometricFactors)[i];
00499 totalMolecularWeight += (*elementVector)[i]->GetAtomicMassAmu() * (*StechiometricFactors)[i];
00500 meanExcitationEnergy += iZ*std::log(meanAtomExcitationEnergy[iZ-1])*(*StechiometricFactors)[i];
00501
00502
00503
00504
00505
00506
00507 std::pair<const G4Material*,G4int> theKey = std::make_pair(material,iZ);
00508 if (!atomTablePerMolecule->count(theKey))
00509 atomTablePerMolecule->insert(std::make_pair(theKey,(*StechiometricFactors)[i]));
00510 }
00511 meanExcitationEnergy = std::exp(meanExcitationEnergy/totalZ);
00512
00513 atomicNumber->insert(std::make_pair(material,totalZ));
00514 atomicMass->insert(std::make_pair(material,totalMolecularWeight));
00515 excitationEnergy->insert(std::make_pair(material,meanExcitationEnergy));
00516 atomsPerMolecule->insert(std::make_pair(material,theatomsPerMolecule));
00517
00518
00519 if (verbosityLevel > 1)
00520 {
00521 G4cout << "Calculated mean excitation energy for " << material->GetName() <<
00522 " = " << meanExcitationEnergy/eV << " eV" << G4endl;
00523 }
00524
00525 std::vector<G4PenelopeOscillator> *helper = new std::vector<G4PenelopeOscillator>;
00526
00527
00528
00529 G4PenelopeOscillator newOsc;
00530 newOsc.SetOscillatorStrength(0.);
00531 newOsc.SetIonisationEnergy(0*eV);
00532 newOsc.SetHartreeFactor(0);
00533 newOsc.SetParentZ(0);
00534 newOsc.SetShellFlag(30);
00535 newOsc.SetParentShellID(30);
00536 helper->push_back(newOsc);
00537
00538
00539 for (G4int k=0;k<nElements;k++)
00540 {
00541 G4double Z = (*elementVector)[k]->GetZ();
00542 G4bool finished = false;
00543 for (G4int i=0;i<2000 && !finished;i++)
00544 {
00545
00546
00547
00548
00549
00550
00551
00552 if (elementData[0][i] == Z)
00553 {
00554 G4int shellID = (G4int) elementData[1][i];
00555 G4double occup = elementData[2][i];
00556 if (shellID > 0)
00557 {
00558 if (std::fabs(occup) > 0)
00559 {
00560 G4PenelopeOscillator newOscLocal;
00561 newOscLocal.SetOscillatorStrength(std::fabs(occup)*(*StechiometricFactors)[k]);
00562 newOscLocal.SetIonisationEnergy(elementData[3][i]);
00563 newOscLocal.SetHartreeFactor(elementData[4][i]/fine_structure_const);
00564 newOscLocal.SetParentZ(elementData[0][i]);
00565
00566 newOscLocal.SetParentShellID((G4int)elementData[1][i]);
00567
00568
00569 if (elementData[0][i] > 6 && elementData[1][i] < 10)
00570 newOscLocal.SetShellFlag(((G4int)elementData[1][i]));
00571 else
00572 newOscLocal.SetShellFlag(30);
00573 helper->push_back(newOscLocal);
00574 if (occup < 0)
00575 {
00576 G4double ff = (*helper)[0].GetOscillatorStrength();
00577 ff += std::fabs(occup)*(*StechiometricFactors)[k];
00578 (*helper)[0].SetOscillatorStrength(ff);
00579 }
00580 }
00581 }
00582
00583 }
00584 if ( elementData[0][i] > Z)
00585 finished = true;
00586 }
00587 }
00588
00589 delete StechiometricFactors;
00590
00591
00592
00593
00594 std::sort(helper->begin(),helper->end());
00595
00596
00597 G4double RydbergEnergy = 13.60569*eV;
00598 G4double Omega = std::sqrt(4*pi*moleculeDensity*totalZ*Bohr_radius)*Bohr_radius*2.0*RydbergEnergy;
00599 G4double conductionStrength = (*helper)[0].GetOscillatorStrength();
00600 G4double plasmaEnergy = Omega*std::sqrt(conductionStrength/totalZ);
00601
00602 plasmaSquared->insert(std::make_pair(material,Omega*Omega));
00603
00604 G4bool isAConductor = false;
00605 G4int nullOsc = 0;
00606
00607 if (verbosityLevel > 1)
00608 {
00609 G4cout << "Estimated oscillator strenght and energy of plasmon: " <<
00610 conductionStrength << " and " << plasmaEnergy/eV << " eV" << G4endl;
00611 }
00612
00613 if (conductionStrength < 0.5 || plasmaEnergy<1.0*eV)
00614 {
00615 if (verbosityLevel >1 )
00616 G4cout << material->GetName() << " is an insulator " << G4endl;
00617
00618 helper->erase(helper->begin());
00619 }
00620 else
00621 {
00622 if (verbosityLevel >1 )
00623 G4cout << material->GetName() << " is a conductor " << G4endl;
00624 isAConductor = true;
00625
00626 G4double conductionStrengthCopy = conductionStrength;
00627 G4bool quit = false;
00628 for (size_t i = 1; i<helper->size() && !quit ;i++)
00629 {
00630 G4double oscStre = (*helper)[i].GetOscillatorStrength();
00631
00632 if (oscStre < conductionStrengthCopy)
00633 {
00634 conductionStrengthCopy = conductionStrengthCopy-oscStre;
00635 (*helper)[i].SetOscillatorStrength(0.);
00636 nullOsc++;
00637 }
00638 else
00639 {
00640 quit = true;
00641 (*helper)[i].SetOscillatorStrength(oscStre-conductionStrengthCopy);
00642 if (std::fabs((*helper)[i].GetOscillatorStrength()) < 1e-12)
00643 {
00644 conductionStrength += (*helper)[i].GetOscillatorStrength();
00645 (*helper)[i].SetOscillatorStrength(0.);
00646 nullOsc++;
00647 }
00648 }
00649 }
00650
00651
00652 (*helper)[0].SetOscillatorStrength(conductionStrength);
00653 (*helper)[0].SetIonisationEnergy(0.);
00654 (*helper)[0].SetResonanceEnergy(plasmaEnergy);
00655 G4double hartree = 0.75/std::sqrt(3.0*pi*pi*moleculeDensity*
00656 Bohr_radius*Bohr_radius*Bohr_radius*conductionStrength);
00657 (*helper)[0].SetHartreeFactor(hartree/fine_structure_const);
00658 }
00659
00660
00661 G4double sum = 0;
00662 for (size_t i=0;i<helper->size();i++)
00663 {
00664 sum += (*helper)[i].GetOscillatorStrength();
00665 }
00666 if (std::fabs(sum-totalZ) > (1e-6*totalZ))
00667 {
00668 G4ExceptionDescription ed;
00669 ed << "Inconsistent oscillator data for " << material->GetName() << G4endl;
00670 ed << sum << " " << totalZ << G4endl;
00671 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
00672 "em2036",FatalException,ed);
00673 }
00674 if (std::fabs(sum-totalZ) > (1e-12*totalZ))
00675 {
00676 G4double fact = totalZ/sum;
00677 for (size_t i=0;i<helper->size();i++)
00678 {
00679 G4double ff = (*helper)[i].GetOscillatorStrength()*fact;
00680 (*helper)[i].SetOscillatorStrength(ff);
00681 }
00682 }
00683
00684
00685 for (G4int k=0;k<nullOsc;k++)
00686 {
00687 G4bool exit=false;
00688 for (size_t i=0;i<helper->size() && !exit;i++)
00689 {
00690 if (std::fabs((*helper)[i].GetOscillatorStrength()) < 1e-12)
00691 {
00692 helper->erase(helper->begin()+i);
00693 exit = true;
00694 }
00695 }
00696 }
00697
00698
00699
00700 G4double adjustmentFactor = 0;
00701 if (helper->size() > 1)
00702 {
00703 G4double TST = totalZ*std::log(meanExcitationEnergy/eV);
00704 G4double AALow = 0.5;
00705 G4double AAHigh = 10.;
00706 do
00707 {
00708 adjustmentFactor = (AALow+AAHigh)*0.5;
00709 G4double sumLocal = 0;
00710 for (size_t i=0;i<helper->size();i++)
00711 {
00712 if (i == 0 && isAConductor)
00713 {
00714 G4double resEne = (*helper)[i].GetResonanceEnergy();
00715 sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
00716 }
00717 else
00718 {
00719 G4double ionEne = (*helper)[i].GetIonisationEnergy();
00720 G4double oscStre = (*helper)[i].GetOscillatorStrength();
00721 G4double WI2 = (adjustmentFactor*adjustmentFactor*ionEne*ionEne) +
00722 2./3.*(oscStre/totalZ)*Omega*Omega;
00723 G4double resEne = std::sqrt(WI2);
00724 (*helper)[i].SetResonanceEnergy(resEne);
00725 sumLocal += (*helper)[i].GetOscillatorStrength()*std::log(resEne/eV);
00726 }
00727 }
00728 if (sumLocal < TST)
00729 AALow = adjustmentFactor;
00730 else
00731 AAHigh = adjustmentFactor;
00732 }while((AAHigh-AALow)>(1e-14*adjustmentFactor));
00733 }
00734 else
00735 {
00736 G4double ionEne = (*helper)[0].GetIonisationEnergy();
00737 (*helper)[0].SetIonisationEnergy(std::fabs(ionEne));
00738 (*helper)[0].SetResonanceEnergy(meanExcitationEnergy);
00739 }
00740 if (verbosityLevel > 1)
00741 {
00742 G4cout << "Sternheimer's adjustment factor: " << adjustmentFactor << G4endl;
00743 }
00744
00745
00746 G4double xcheck = (*helper)[0].GetOscillatorStrength()*std::log((*helper)[0].GetResonanceEnergy());
00747 G4double TST = (*helper)[0].GetOscillatorStrength();
00748 for (size_t i=1;i<helper->size();i++)
00749 {
00750 xcheck += (*helper)[i].GetOscillatorStrength()*std::log((*helper)[i].GetResonanceEnergy());
00751 TST += (*helper)[i].GetOscillatorStrength();
00752 }
00753 if (std::fabs(TST-totalZ)>1e-8*totalZ)
00754 {
00755 G4ExceptionDescription ed;
00756 ed << "Inconsistent oscillator data " << G4endl;
00757 ed << TST << " " << totalZ << G4endl;
00758 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
00759 "em2036",FatalException,ed);
00760 }
00761 xcheck = std::exp(xcheck/totalZ);
00762 if (std::fabs(xcheck-meanExcitationEnergy) > 1e-8*meanExcitationEnergy)
00763 {
00764 G4ExceptionDescription ed;
00765 ed << "Error in Sterheimer factor calculation " << G4endl;
00766 ed << xcheck/eV << " " << meanExcitationEnergy/eV << G4endl;
00767 G4Exception("G4PenelopeOscillatorManager::BuildOscillatorTable()",
00768 "em2037",FatalException,ed);
00769 }
00770
00771
00772
00773
00774
00775 G4double Zmax = 0;
00776 for (G4int k=0;k<nElements;k++)
00777 {
00778 G4double Z = (*elementVector)[k]->GetZ();
00779 if (Z>Zmax) Zmax = Z;
00780 }
00781
00782 G4bool found = false;
00783 G4double cutEnergy = 50*eV;
00784 for (size_t i=0;i<helper->size() && !found;i++)
00785 {
00786 G4double Z = (*helper)[i].GetParentZ();
00787 G4int shID = (*helper)[i].GetParentShellID();
00788 if (shID == 10 && Z == Zmax)
00789 {
00790 found = true;
00791 if ((*helper)[i].GetIonisationEnergy() > cutEnergy)
00792 cutEnergy = (*helper)[i].GetIonisationEnergy();
00793 }
00794 }
00795
00796
00797 G4double lowEnergyLimitForFluorescence = 250*eV;
00798 cutEnergy = std::min(cutEnergy,lowEnergyLimitForFluorescence);
00799
00800 if (verbosityLevel > 1)
00801 G4cout << "Cutoff energy: " << cutEnergy/eV << " eV" << G4endl;
00802
00803
00804
00805
00806
00807 G4PenelopeOscillatorTable* theTable = new G4PenelopeOscillatorTable();
00808 G4PenelopeOscillatorResEnergyComparator comparator;
00809 std::sort(helper->begin(),helper->end(),comparator);
00810
00811
00812 for (size_t i=0;i<helper->size();i++)
00813 {
00814
00815 G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
00816 theTable->push_back(theOsc);
00817 }
00818
00819
00820
00821 G4double Rgroup = 1.05;
00822 size_t Nost = theTable->size();
00823
00824 size_t firstIndex = (isAConductor) ? 1 : 0;
00825 G4bool loopAgain = false;
00826 G4int removedLevels = 0;
00827 do
00828 {
00829 loopAgain = false;
00830 if (Nost>firstIndex+1)
00831 {
00832 removedLevels = 0;
00833 for (size_t i=firstIndex;i<theTable->size()-1;i++)
00834 {
00835 G4bool skipLoop = false;
00836 G4int shellFlag = (*theTable)[i]->GetShellFlag();
00837 G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
00838 G4double resEne = (*theTable)[i]->GetResonanceEnergy();
00839 G4double resEnePlus1 = (*theTable)[i+1]->GetResonanceEnergy();
00840 G4double oscStre = (*theTable)[i]->GetOscillatorStrength();
00841 G4double oscStrePlus1 = (*theTable)[i+1]->GetOscillatorStrength();
00842
00843 if (ionEne>cutEnergy)
00844 skipLoop = true;
00845 if (resEne<1.0*eV || resEnePlus1<1.0*eV)
00846 skipLoop = true;
00847 if (resEnePlus1 > Rgroup*resEne)
00848 skipLoop = true;
00849 if (!skipLoop)
00850 {
00851 G4double newRes = std::exp((oscStre*std::log(resEne)+
00852 oscStrePlus1*std::log(resEnePlus1))
00853 /(oscStre+oscStrePlus1));
00854 (*theTable)[i]->SetResonanceEnergy(newRes);
00855 G4double newIon = (oscStre*ionEne+
00856 oscStrePlus1*(*theTable)[i+1]->GetIonisationEnergy())/
00857 (oscStre+oscStrePlus1);
00858 (*theTable)[i]->SetIonisationEnergy(newIon);
00859 G4double newStre = oscStre+oscStrePlus1;
00860 (*theTable)[i]->SetOscillatorStrength(newStre);
00861 G4double newHartree = (oscStre*(*theTable)[i]->GetHartreeFactor()+
00862 oscStrePlus1*(*theTable)[i+1]->GetHartreeFactor())/
00863 (oscStre+oscStrePlus1);
00864 (*theTable)[i]->SetHartreeFactor(newHartree);
00865 if ((*theTable)[i]->GetParentZ() != (*theTable)[i+1]->GetParentZ())
00866 (*theTable)[i]->SetParentZ(0.);
00867 if (shellFlag < 10 || (*theTable)[i+1]->GetShellFlag() < 10)
00868 {
00869 G4int newFlag = std::min(shellFlag,(*theTable)[i+1]->GetShellFlag());
00870 (*theTable)[i]->SetShellFlag(newFlag);
00871 }
00872 else
00873 (*theTable)[i]->SetShellFlag(30);
00874
00875 (*theTable)[i]->SetParentShellID((*theTable)[i]->GetShellFlag());
00876
00877
00878 if (i<theTable->size()-2)
00879 {
00880 for (size_t ii=i+1;ii<theTable->size()-1;ii++)
00881 (*theTable)[ii] = (*theTable)[ii+1];
00882 }
00883
00884 theTable->erase(theTable->begin()+theTable->size()-1);
00885 removedLevels++;
00886 }
00887 }
00888 }
00889 if (removedLevels)
00890 {
00891 Nost -= removedLevels;
00892 loopAgain = true;
00893 }
00894 if (Rgroup < 1.414213 || Nost > 64)
00895 {
00896 Rgroup = Rgroup*Rgroup;
00897 loopAgain = true;
00898 }
00899 }while(loopAgain);
00900
00901 if (verbosityLevel > 1)
00902 {
00903 G4cout << "Final grouping factor for Ionisation: " << Rgroup << G4endl;
00904 }
00905
00906
00907 for (size_t i=0;i<theTable->size();i++)
00908 {
00909
00910 G4double ionEne = (*theTable)[i]->GetIonisationEnergy();
00911 if (ionEne < 1e-3*eV)
00912 {
00913 G4double resEne = (*theTable)[i]->GetResonanceEnergy();
00914 (*theTable)[i]->SetIonisationEnergy(0.*eV);
00915 (*theTable)[i]->SetCutoffRecoilResonantEnergy(resEne);
00916 }
00917 else
00918 (*theTable)[i]->SetCutoffRecoilResonantEnergy(ionEne);
00919 }
00920
00921
00922 oscillatorStoreIonisation->insert(std::make_pair(material,theTable));
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932 G4PenelopeOscillatorTable* theTableC = new G4PenelopeOscillatorTable();
00933
00934 std::sort(helper->begin(),helper->end());
00935
00936 for (size_t i=0;i<helper->size();i++)
00937 {
00938
00939 G4PenelopeOscillator* theOsc = new G4PenelopeOscillator((*helper)[i]);
00940 theTableC->push_back(theOsc);
00941 }
00942
00943
00944 Rgroup = 1.5;
00945 Nost = theTableC->size();
00946
00947 firstIndex = (isAConductor) ? 1 : 0;
00948 loopAgain = false;
00949 removedLevels = 0;
00950 do
00951 {
00952 loopAgain = false;
00953 if (Nost>firstIndex+1)
00954 {
00955 removedLevels = 0;
00956 for (size_t i=firstIndex;i<theTableC->size()-1;i++)
00957 {
00958 G4bool skipLoop = false;
00959
00960 G4double ionEne = (*theTableC)[i]->GetIonisationEnergy();
00961 G4double ionEnePlus1 = (*theTableC)[i+1]->GetIonisationEnergy();
00962 G4double oscStre = (*theTableC)[i]->GetOscillatorStrength();
00963 G4double oscStrePlus1 = (*theTableC)[i+1]->GetOscillatorStrength();
00964
00965 if (ionEne>cutEnergy)
00966 skipLoop = true;
00967 if (ionEne<1.0*eV || ionEnePlus1<1.0*eV)
00968 skipLoop = true;
00969 if (ionEnePlus1 > Rgroup*ionEne)
00970 skipLoop = true;
00971
00972 if (!skipLoop)
00973 {
00974 G4double newIon = (oscStre*ionEne+
00975 oscStrePlus1*ionEnePlus1)/
00976 (oscStre+oscStrePlus1);
00977 (*theTableC)[i]->SetIonisationEnergy(newIon);
00978 G4double newStre = oscStre+oscStrePlus1;
00979 (*theTableC)[i]->SetOscillatorStrength(newStre);
00980 G4double newHartree = (oscStre*(*theTableC)[i]->GetHartreeFactor()+
00981 oscStrePlus1*(*theTableC)[i+1]->GetHartreeFactor())/
00982 (oscStre+oscStrePlus1);
00983 (*theTableC)[i]->SetHartreeFactor(newHartree);
00984 if ((*theTableC)[i]->GetParentZ() != (*theTableC)[i+1]->GetParentZ())
00985 (*theTableC)[i]->SetParentZ(0.);
00986 (*theTableC)[i]->SetShellFlag(30);
00987 (*theTableC)[i]->SetParentShellID((*theTableC)[i]->GetShellFlag());
00988
00989 if (i<theTableC->size()-2)
00990 {
00991 for (size_t ii=i+1;ii<theTableC->size()-1;ii++)
00992 (*theTableC)[ii] = (*theTableC)[ii+1];
00993 }
00994 theTableC->erase(theTableC->begin()+theTableC->size()-1);
00995 removedLevels++;
00996 }
00997 }
00998 }
00999 if (removedLevels)
01000 {
01001 Nost -= removedLevels;
01002 loopAgain = true;
01003 }
01004 if (Rgroup < 2.0 || Nost > 64)
01005 {
01006 Rgroup = Rgroup*Rgroup;
01007 loopAgain = true;
01008 }
01009 }while(loopAgain);
01010
01011
01012 if (verbosityLevel > 1)
01013 {
01014 G4cout << "Final grouping factor for Compton: " << Rgroup << G4endl;
01015 }
01016
01017
01018 oscillatorStoreCompton->insert(std::make_pair(material,theTableC));
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
01045
01046
01047 delete helper;
01048 if (verbosityLevel > 1)
01049 Dump(material);
01050
01051 return;
01052 }
01053
01054
01055
01056 void G4PenelopeOscillatorManager::ReadElementData()
01057 {
01058 if (verbosityLevel > 0)
01059 {
01060 G4cout << "G4PenelopeOscillatorManager::ReadElementData()" << G4endl;
01061 G4cout << "Going to read Element Data" << G4endl;
01062 }
01063 char* path = getenv("G4LEDATA");
01064 if (!path)
01065 {
01066 G4String excep = "G4PenelopeOscillatorManager - G4LEDATA environment variable not set!";
01067 G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
01068 "em0006",FatalException,excep);
01069 return;
01070 }
01071 G4String pathString(path);
01072 G4String pathFile = pathString + "/penelope/pdatconf.p08";
01073 std::ifstream file(pathFile);
01074
01075 if (!file.is_open())
01076 {
01077 G4String excep = "G4PenelopeOscillatorManager - data file " + pathFile + " not found!";
01078 G4Exception("G4PenelopeOscillatorManager::ReadElementData()",
01079 "em0003",FatalException,excep);
01080 }
01081
01082 G4AtomicTransitionManager* theTransitionManager =
01083 G4AtomicTransitionManager::Instance();
01084
01085
01086 G4String theHeader;
01087 for (G4int iline=0;iline<22;iline++)
01088 getline(file,theHeader);
01089
01090 G4int Z=0;
01091 G4int shellCode = 0;
01092 G4String shellId = "NULL";
01093 G4int occupationNumber = 0;
01094 G4double ionisationEnergy = 0.0*eV;
01095 G4double hartreeProfile = 0.;
01096 G4int shellCounter = 0;
01097 G4int oldZ = -1;
01098 G4int numberOfShells = 0;
01099
01100 for (G4int i=0;!file.eof();i++)
01101 {
01102 file >> Z >> shellCode >> shellId >> occupationNumber >> ionisationEnergy >> hartreeProfile;
01103 if (Z>0 && i<2000)
01104 {
01105 elementData[0][i] = Z;
01106 elementData[1][i] = shellCode;
01107 elementData[2][i] = occupationNumber;
01108
01109 if (Z != oldZ)
01110 {
01111 shellCounter = 0;
01112 oldZ = Z;
01113 numberOfShells = theTransitionManager->NumberOfShells(Z);
01114 }
01115 G4double bindingEnergy = -1*eV;
01116 if (shellCounter<numberOfShells)
01117 {
01118 G4AtomicShell* shell = theTransitionManager->Shell(Z,shellCounter);
01119 bindingEnergy = shell->BindingEnergy();
01120 }
01121
01122
01123 elementData[3][i] = (bindingEnergy>100*eV) ? bindingEnergy : ionisationEnergy*eV;
01124
01125 elementData[4][i] = hartreeProfile;
01126 shellCounter++;
01127 }
01128 }
01129 file.close();
01130
01131 if (verbosityLevel > 1)
01132 {
01133 G4cout << "G4PenelopeOscillatorManager::ReadElementData(): Data file read" << G4endl;
01134 }
01135 fReadElementData = true;
01136 return;
01137
01138 }
01139
01140
01141 G4double G4PenelopeOscillatorManager::GetMeanExcitationEnergy(const G4Material* mat)
01142 {
01143
01144 CheckForTablesCreated();
01145
01146
01147 if (excitationEnergy->count(mat))
01148 return excitationEnergy->find(mat)->second;
01149
01150
01151 BuildOscillatorTable(mat);
01152
01153
01154 if (excitationEnergy->count(mat))
01155 return excitationEnergy->find(mat)->second;
01156 else
01157 {
01158 G4cout << "G4PenelopeOscillatorManager::GetMolecularExcitationEnergy() " << G4endl;
01159 G4cout << "Impossible to retrieve the excitation energy for " << mat->GetName() << G4endl;
01160 return 0;
01161 }
01162 }
01163
01164
01165 G4double G4PenelopeOscillatorManager::GetPlasmaEnergySquared(const G4Material* mat)
01166 {
01167
01168 CheckForTablesCreated();
01169
01170
01171 if (plasmaSquared->count(mat))
01172 return plasmaSquared->find(mat)->second;
01173
01174
01175 BuildOscillatorTable(mat);
01176
01177
01178 if (plasmaSquared->count(mat))
01179 return plasmaSquared->find(mat)->second;
01180 else
01181 {
01182 G4cout << "G4PenelopeOscillatorManager::GetPlasmaEnergySquared() " << G4endl;
01183 G4cout << "Impossible to retrieve the plasma energy for " << mat->GetName() << G4endl;
01184 return 0;
01185 }
01186 }
01187
01188
01189
01190 G4double G4PenelopeOscillatorManager::GetAtomsPerMolecule(const G4Material* mat)
01191 {
01192
01193 CheckForTablesCreated();
01194
01195
01196 if (atomsPerMolecule->count(mat))
01197 return atomsPerMolecule->find(mat)->second;
01198
01199
01200 BuildOscillatorTable(mat);
01201
01202
01203 if (atomsPerMolecule->count(mat))
01204 return atomsPerMolecule->find(mat)->second;
01205 else
01206 {
01207 G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
01208 G4cout << "Impossible to retrieve the number of atoms per molecule for "
01209 << mat->GetName() << G4endl;
01210 return 0;
01211 }
01212 }
01213
01214
01215
01216 G4double G4PenelopeOscillatorManager::GetNumberOfZAtomsPerMolecule(const G4Material* mat,G4int Z)
01217 {
01218
01219 CheckForTablesCreated();
01220
01221
01222 std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
01223 if (atomTablePerMolecule->count(theKey))
01224 return atomTablePerMolecule->find(theKey)->second;
01225
01226
01227 BuildOscillatorTable(mat);
01228
01229
01230 if (atomTablePerMolecule->count(theKey))
01231 return atomTablePerMolecule->find(theKey)->second;
01232 else
01233 {
01234 G4cout << "G4PenelopeOscillatorManager::GetAtomsPerMolecule() " << G4endl;
01235 G4cout << "Impossible to retrieve the number of atoms per molecule for Z = "
01236 << Z << " in material " << mat->GetName() << G4endl;
01237 return 0;
01238 }
01239 }