Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4PenelopePhotoElectricModel.cc 81067 2014-05-20 09:19:32Z gcosmo $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 08 Jan 2010 L Pandola First implementati 32 // 08 Jan 2010 L Pandola First implementation 32 // 01 Feb 2011 L Pandola Suppress fake ener << 33 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger 33 // is active. 34 // is active. 34 // Make sure that flu << 35 // Make sure that fluorescence/Auger is generated 35 // only if above thre 36 // only if above threshold 36 // 25 May 2011 L Pandola Renamed (make v200 37 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope) 37 // 10 Jun 2011 L Pandola Migrate atomic dee 38 // 10 Jun 2011 L Pandola Migrate atomic deexcitation interface 38 // 07 Oct 2011 L Pandola Bug fix (potential 39 // 07 Oct 2011 L Pandola Bug fix (potential violation of energy conservation) 39 // 27 Sep 2013 L Pandola Migrate to MT para << 40 // 27 Sep 2013 L Pandola Migrate to MT paradigm, only master model manages 40 // tables. 41 // tables. 41 // 02 Oct 2013 L Pandola Rewrite sampling a 42 // 02 Oct 2013 L Pandola Rewrite sampling algorithm of SelectRandomShell() 42 // to improve CPU per 43 // to improve CPU performances 43 // 44 // 44 45 45 #include "G4PenelopePhotoElectricModel.hh" 46 #include "G4PenelopePhotoElectricModel.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh" 48 #include "G4ParticleDefinition.hh" 49 #include "G4ParticleDefinition.hh" 49 #include "G4MaterialCutsCouple.hh" 50 #include "G4MaterialCutsCouple.hh" 50 #include "G4DynamicParticle.hh" 51 #include "G4DynamicParticle.hh" 51 #include "G4PhysicsTable.hh" 52 #include "G4PhysicsTable.hh" 52 #include "G4PhysicsFreeVector.hh" 53 #include "G4PhysicsFreeVector.hh" 53 #include "G4ElementTable.hh" 54 #include "G4ElementTable.hh" 54 #include "G4Element.hh" 55 #include "G4Element.hh" 55 #include "G4AtomicTransitionManager.hh" 56 #include "G4AtomicTransitionManager.hh" 56 #include "G4AtomicShell.hh" 57 #include "G4AtomicShell.hh" 57 #include "G4Gamma.hh" 58 #include "G4Gamma.hh" 58 #include "G4Electron.hh" 59 #include "G4Electron.hh" 59 #include "G4AutoLock.hh" 60 #include "G4AutoLock.hh" 60 #include "G4LossTableManager.hh" 61 #include "G4LossTableManager.hh" 61 #include "G4Exp.hh" << 62 << 63 //....oooOO0OOooo........oooOO0OOooo........oo << 64 << 65 const G4int G4PenelopePhotoElectricModel::fMax << 66 G4PhysicsTable* G4PenelopePhotoElectricModel:: << 67 62 68 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 64 70 G4PenelopePhotoElectricModel::G4PenelopePhotoE 65 G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition* part, 71 const G4String& nam) 66 const G4String& nam) 72 :G4VEmModel(nam),fParticleChange(nullptr),fP << 67 :G4VEmModel(nam),fParticleChange(0),fParticle(0), 73 fAtomDeexcitation(nullptr),fIsInitialised(f << 68 isInitialised(false),fAtomDeexcitation(0),logAtomicShellXS(0),fLocalTable(false) 74 { 69 { 75 fIntrinsicLowEnergyLimit = 100.0*eV; 70 fIntrinsicLowEnergyLimit = 100.0*eV; 76 fIntrinsicHighEnergyLimit = 100.0*GeV; 71 fIntrinsicHighEnergyLimit = 100.0*GeV; 77 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 72 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 73 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 79 // 74 // 80 75 81 if (part) 76 if (part) 82 SetParticle(part); 77 SetParticle(part); 83 78 84 fVerboseLevel= 0; << 79 verboseLevel= 0; 85 // Verbosity scale: 80 // Verbosity scale: 86 // 0 = nothing << 81 // 0 = nothing 87 // 1 = warning for energy non-conservation << 82 // 1 = warning for energy non-conservation 88 // 2 = details of energy budget 83 // 2 = details of energy budget 89 // 3 = calculation of cross sections, file o 84 // 3 = calculation of cross sections, file openings, sampling of atoms 90 // 4 = entering in methods 85 // 4 = entering in methods 91 86 92 //Mark this model as "applicable" for atomic 87 //Mark this model as "applicable" for atomic deexcitation 93 SetDeexcitationFlag(true); 88 SetDeexcitationFlag(true); 94 89 95 fTransitionManager = G4AtomicTransitionManag 90 fTransitionManager = G4AtomicTransitionManager::Instance(); 96 } 91 } 97 92 98 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 99 94 100 G4PenelopePhotoElectricModel::~G4PenelopePhoto 95 G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel() 101 { << 96 { 102 if (IsMaster() || fLocalTable) 97 if (IsMaster() || fLocalTable) 103 { 98 { 104 for(G4int i=0; i<=fMaxZ; ++i) << 99 std::map <G4int,G4PhysicsTable*>::iterator i; >> 100 if (logAtomicShellXS) 105 { 101 { 106 if(fLogAtomicShellXS[i]) { << 102 for (i=logAtomicShellXS->begin();i != logAtomicShellXS->end();i++) 107 fLogAtomicShellXS[i]->clearAndDestroy(); << 103 { 108 delete fLogAtomicShellXS[i]; << 104 G4PhysicsTable* tab = i->second; 109 fLogAtomicShellXS[i] = nullptr; << 105 //tab->clearAndDestroy(); 110 } << 106 delete tab; >> 107 } 111 } 108 } >> 109 delete logAtomicShellXS; 112 } 110 } 113 } 111 } 114 112 115 //....oooOO0OOooo........oooOO0OOooo........oo 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 114 117 void G4PenelopePhotoElectricModel::Initialise( 115 void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition* particle, 118 const G4DataVector& cuts) 116 const G4DataVector& cuts) 119 { 117 { 120 if (fVerboseLevel > 3) << 118 if (verboseLevel > 3) 121 G4cout << "Calling G4PenelopePhotoElectri 119 G4cout << "Calling G4PenelopePhotoElectricModel::Initialise()" << G4endl; 122 120 123 fAtomDeexcitation = G4LossTableManager::Inst 121 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 124 //Issue warning if the AtomicDeexcitation ha 122 //Issue warning if the AtomicDeexcitation has not been declared 125 if (!fAtomDeexcitation) 123 if (!fAtomDeexcitation) 126 { 124 { 127 G4cout << G4endl; 125 G4cout << G4endl; 128 G4cout << "WARNING from G4PenelopePhotoE 126 G4cout << "WARNING from G4PenelopePhotoElectricModel " << G4endl; 129 G4cout << "Atomic de-excitation module i 127 G4cout << "Atomic de-excitation module is not instantiated, so there will not be "; 130 G4cout << "any fluorescence/Auger emissi 128 G4cout << "any fluorescence/Auger emission." << G4endl; 131 G4cout << "Please make sure this is inte 129 G4cout << "Please make sure this is intended" << G4endl; 132 } 130 } 133 131 134 SetParticle(particle); 132 SetParticle(particle); 135 133 136 //Only the master model creates/fills/destro 134 //Only the master model creates/fills/destroys the tables 137 if (IsMaster() && particle == fParticle) << 135 if (IsMaster() && particle == fParticle) 138 { 136 { 139 G4ProductionCutsTable* theCoupleTable = << 140 G4ProductionCutsTable::GetProductionCutsTabl << 141 137 142 for (G4int i=0;i<(G4int)theCoupleTable-> << 138 // logAtomicShellXS is created only once, since it is never cleared >> 139 if (!logAtomicShellXS) >> 140 logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>; >> 141 >> 142 G4ProductionCutsTable* theCoupleTable = >> 143 G4ProductionCutsTable::GetProductionCutsTable(); >> 144 >> 145 for (size_t i=0;i<theCoupleTable->GetTableSize();i++) 143 { 146 { 144 const G4Material* material = << 147 const G4Material* material = 145 theCoupleTable->GetMaterialCutsCouple(i) 148 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 146 const G4ElementVector* theElementVector = 149 const G4ElementVector* theElementVector = material->GetElementVector(); 147 << 150 148 for (std::size_t j=0;j<material->GetNumber << 151 for (size_t j=0;j<material->GetNumberOfElements();j++) 149 { 152 { 150 G4int iZ = theElementVector->at(j)->Ge << 153 G4int iZ = (G4int) theElementVector->at(j)->GetZ(); 151 //read data files only in the master 154 //read data files only in the master 152 if (!fLogAtomicShellXS[iZ]) << 155 if (!logAtomicShellXS->count(iZ)) 153 ReadDataFile(iZ); << 156 ReadDataFile(iZ); 154 } 157 } 155 } 158 } >> 159 156 160 157 InitialiseElementSelectors(particle,cuts << 161 InitialiseElementSelectors(particle,cuts); 158 << 162 159 if (fVerboseLevel > 0) { << 163 if (verboseLevel > 0) { 160 G4cout << "Penelope Photo-Electric model v20 164 G4cout << "Penelope Photo-Electric model v2008 is initialized " << G4endl 161 << "Energy range: " 165 << "Energy range: " 162 << LowEnergyLimit() / MeV << " MeV - 166 << LowEnergyLimit() / MeV << " MeV - " 163 << HighEnergyLimit() / GeV << " GeV"; 167 << HighEnergyLimit() / GeV << " GeV"; 164 } 168 } 165 } 169 } 166 170 167 if(fIsInitialised) return; << 171 if(isInitialised) return; 168 fParticleChange = GetParticleChangeForGamma( 172 fParticleChange = GetParticleChangeForGamma(); 169 fIsInitialised = true; << 173 isInitialised = true; 170 174 171 } 175 } 172 176 173 void G4PenelopePhotoElectricModel::InitialiseL 177 void G4PenelopePhotoElectricModel::InitialiseLocal(const G4ParticleDefinition* part, 174 G4VEmModel *masterModel) 178 G4VEmModel *masterModel) 175 { 179 { 176 if (fVerboseLevel > 3) << 180 if (verboseLevel > 3) 177 G4cout << "Calling G4PenelopePhotoElectri 181 G4cout << "Calling G4PenelopePhotoElectricModel::InitialiseLocal()" << G4endl; >> 182 178 // 183 // 179 //Check that particle matches: one might hav << 184 //Check that particle matches: one might have multiple master models (e.g. 180 //for e+ and e-). 185 //for e+ and e-). 181 // 186 // 182 if (part == fParticle) 187 if (part == fParticle) 183 { 188 { 184 SetElementSelectors(masterModel->GetElem 189 SetElementSelectors(masterModel->GetElementSelectors()); 185 190 186 //Get the const table pointers from the 191 //Get the const table pointers from the master to the workers 187 const G4PenelopePhotoElectricModel* theM << 192 const G4PenelopePhotoElectricModel* theModel = 188 static_cast<G4PenelopePhotoElectricModel*> ( 193 static_cast<G4PenelopePhotoElectricModel*> (masterModel); 189 for(G4int i=0; i<=fMaxZ; ++i) << 194 190 fLogAtomicShellXS[i] = theModel->fLogAtomicS << 195 logAtomicShellXS = theModel->logAtomicShellXS; >> 196 191 //Same verbosity for all workers, as the 197 //Same verbosity for all workers, as the master 192 fVerboseLevel = theModel->fVerboseLevel; << 198 verboseLevel = theModel->verboseLevel; 193 } 199 } 194 200 195 return; 201 return; 196 } 202 } 197 203 198 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 199 namespace { G4Mutex PenelopePhotoElectricMode 205 namespace { G4Mutex PenelopePhotoElectricModelMutex = G4MUTEX_INITIALIZER; } 200 G4double G4PenelopePhotoElectricModel::Compute 206 G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom( 201 const G4ParticleDefinition*, 207 const G4ParticleDefinition*, 202 G4double energy, 208 G4double energy, 203 G4double Z, G4double, 209 G4double Z, G4double, 204 G4double, G4double) 210 G4double, G4double) 205 { 211 { 206 // 212 // 207 // Penelope model v2008 213 // Penelope model v2008 208 // 214 // 209 if (fVerboseLevel > 3) << 215 >> 216 if (verboseLevel > 3) 210 G4cout << "Calling ComputeCrossSectionPerA 217 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl; 211 218 212 G4int iZ = G4int(Z); << 219 G4int iZ = (G4int) Z; 213 220 214 if (!fLogAtomicShellXS[iZ]) << 221 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was >> 222 //not invoked >> 223 if (!logAtomicShellXS) >> 224 { >> 225 //create a **thread-local** version of the table. Used only for G4EmCalculator and >> 226 //Unit Tests >> 227 fLocalTable = true; >> 228 logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>; >> 229 } >> 230 >> 231 //now it should be ok >> 232 if (!logAtomicShellXS->count(iZ)) 215 { 233 { 216 //If we are here, it means that Initiali << 234 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 217 //not filled up. This can happen in a Un 235 //not filled up. This can happen in a UnitTest or via G4EmCalculator 218 if (fVerboseLevel > 0) << 236 if (verboseLevel > 0) 219 { << 237 { 220 //Issue a G4Exception (warning) only in ve 238 //Issue a G4Exception (warning) only in verbose mode 221 G4ExceptionDescription ed; 239 G4ExceptionDescription ed; 222 ed << "Unable to retrieve the shell cross 240 ed << "Unable to retrieve the shell cross section table for Z=" << iZ << G4endl; 223 ed << "This can happen only in Unit Tests 241 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl; 224 G4Exception("G4PenelopePhotoElectricModel: 242 G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()", 225 "em2038",JustWarning,ed); 243 "em2038",JustWarning,ed); 226 } 244 } 227 //protect file reading via autolock 245 //protect file reading via autolock 228 G4AutoLock lock(&PenelopePhotoElectricMo 246 G4AutoLock lock(&PenelopePhotoElectricModelMutex); 229 ReadDataFile(iZ); 247 ReadDataFile(iZ); 230 lock.unlock(); 248 lock.unlock(); >> 249 231 } 250 } 232 251 233 G4double cross = 0; 252 G4double cross = 0; 234 G4PhysicsTable* theTable = fLogAtomicShellX << 253 >> 254 G4PhysicsTable* theTable = logAtomicShellXS->find(iZ)->second; 235 G4PhysicsFreeVector* totalXSLog = (G4Physics 255 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0]; 236 256 237 if (!totalXSLog) 257 if (!totalXSLog) 238 { 258 { 239 G4Exception("G4PenelopePhotoElectricMod 259 G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()", 240 "em2039",FatalException, 260 "em2039",FatalException, 241 "Unable to retrieve the total cross sec << 261 "Unable to retrieve the total cross section table"); 242 return 0; 262 return 0; 243 } 263 } 244 G4double logene = G4Log(energy); << 264 G4double logene = std::log(energy); 245 G4double logXS = totalXSLog->Value(logene); 265 G4double logXS = totalXSLog->Value(logene); 246 cross = G4Exp(logXS); << 266 cross = std::exp(logXS); 247 << 267 248 if (fVerboseLevel > 2) << 268 if (verboseLevel > 2) 249 G4cout << "Photoelectric cross section at 269 G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z << 250 " = " << cross/barn << " barn" << G4endl 270 " = " << cross/barn << " barn" << G4endl; 251 return cross; 271 return cross; 252 } 272 } 253 273 254 //....oooOO0OOooo........oooOO0OOooo........oo 274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 255 275 256 void G4PenelopePhotoElectricModel::SampleSecon 276 void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 257 const G4MaterialCutsCouple* c 277 const G4MaterialCutsCouple* couple, 258 const G4DynamicParticle* aDyn 278 const G4DynamicParticle* aDynamicGamma, 259 G4double, 279 G4double, 260 G4double) 280 G4double) 261 { 281 { 262 // 282 // 263 // Photoelectric effect, Penelope model v200 283 // Photoelectric effect, Penelope model v2008 264 // 284 // 265 // The target atom and the target shell are << 285 // The target atom and the target shell are sampled according to the Livermore 266 // database << 286 // database 267 // D.E. Cullen et al., Report UCRL-50400 (1 287 // D.E. Cullen et al., Report UCRL-50400 (1989) 268 // The angular distribution of the electron << 288 // The angular distribution of the electron in the final state is sampled 269 // according to the Sauter distribution from << 289 // according to the Sauter distribution from 270 // F. Sauter, Ann. Phys. 11 (1931) 454 290 // F. Sauter, Ann. Phys. 11 (1931) 454 271 // The energy of the final electron is given << 291 // The energy of the final electron is given by the initial photon energy minus 272 // the binding energy. Fluorescence de-excit << 292 // the binding energy. Fluorescence de-excitation is subsequently produced 273 // (to fill the vacancy) according to the ge 293 // (to fill the vacancy) according to the general Geant4 G4DeexcitationManager: 274 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1 294 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1-1-9 (1997) 275 295 276 if (fVerboseLevel > 3) << 296 if (verboseLevel > 3) 277 G4cout << "Calling SamplingSecondaries() o 297 G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl; 278 298 279 G4double photonEnergy = aDynamicGamma->GetKi 299 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 280 300 281 // always kill primary 301 // always kill primary 282 fParticleChange->ProposeTrackStatus(fStopAnd 302 fParticleChange->ProposeTrackStatus(fStopAndKill); 283 fParticleChange->SetProposedKineticEnergy(0. 303 fParticleChange->SetProposedKineticEnergy(0.); 284 304 285 if (photonEnergy <= fIntrinsicLowEnergyLimit 305 if (photonEnergy <= fIntrinsicLowEnergyLimit) 286 { 306 { 287 fParticleChange->ProposeLocalEnergyDepos 307 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy); 288 return ; 308 return ; 289 } 309 } 290 310 291 G4ParticleMomentum photonDirection = aDynami 311 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); 292 312 293 // Select randomly one element in the curren 313 // Select randomly one element in the current material 294 if (fVerboseLevel > 2) << 314 if (verboseLevel > 2) 295 G4cout << "Going to select element in " << 315 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl; 296 316 297 // atom can be selected efficiently if eleme 317 // atom can be selected efficiently if element selectors are initialised 298 const G4Element* anElement = 318 const G4Element* anElement = 299 SelectRandomAtom(couple,G4Gamma::GammaDefi 319 SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy); 300 G4int Z = anElement->GetZasInt(); << 320 G4int Z = (G4int) anElement->GetZ(); 301 if (fVerboseLevel > 2) << 321 if (verboseLevel > 2) 302 G4cout << "Selected " << anElement->GetNam 322 G4cout << "Selected " << anElement->GetName() << G4endl; 303 << 323 304 // Select the ionised shell in the current a 324 // Select the ionised shell in the current atom according to shell cross sections 305 //shellIndex = 0 --> K shell 325 //shellIndex = 0 --> K shell 306 // 1-3 --> L shells 326 // 1-3 --> L shells 307 // 4-8 --> M shells 327 // 4-8 --> M shells 308 // 9 --> outer shells cumulative 328 // 9 --> outer shells cumulatively 309 // 329 // 310 std::size_t shellIndex = SelectRandomShell(Z << 330 size_t shellIndex = SelectRandomShell(Z,photonEnergy); 311 331 312 if (fVerboseLevel > 2) << 332 if (verboseLevel > 2) 313 G4cout << "Selected shell " << shellIndex 333 G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl; 314 334 315 // Retrieve the corresponding identifier and 335 // Retrieve the corresponding identifier and binding energy of the selected shell 316 const G4AtomicTransitionManager* transitionM 336 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance(); 317 337 318 //The number of shell cross section possibly << 338 //The number of shell cross section possibly reported in the Penelope database 319 //might be different from the number of shel 339 //might be different from the number of shells in the G4AtomicTransitionManager 320 //(namely, Penelope may contain more shell, 340 //(namely, Penelope may contain more shell, especially for very light elements). 321 //In order to avoid a warning message from t << 341 //In order to avoid a warning message from the G4AtomicTransitionManager, I 322 //add this protection. Results are anyway ch 342 //add this protection. Results are anyway changed, because when G4AtomicTransitionManager 323 //has a shellID>maxID, it sets the shellID t << 343 //has a shellID>maxID, it sets the shellID to the last valid shell. 324 std::size_t numberOfShells = (std::size_t) t << 344 size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z); 325 if (shellIndex >= numberOfShells) 345 if (shellIndex >= numberOfShells) 326 shellIndex = numberOfShells-1; 346 shellIndex = numberOfShells-1; 327 347 328 const G4AtomicShell* shell = fTransitionMana 348 const G4AtomicShell* shell = fTransitionManager->Shell(Z,shellIndex); 329 G4double bindingEnergy = shell->BindingEnerg 349 G4double bindingEnergy = shell->BindingEnergy(); >> 350 //G4int shellId = shell->ShellId(); 330 351 331 //Penelope considers only K, L and M shells. << 352 //Penelope considers only K, L and M shells. Cross sections of outer shells are 332 //not included in the Penelope database. If << 353 //not included in the Penelope database. If SelectRandomShell() returns 333 //shellIndex = 9, it means that an outer she << 354 //shellIndex = 9, it means that an outer shell was ionized. In this case the 334 //Penelope recipe is to set bindingEnergy = << 355 //Penelope recipe is to set bindingEnergy = 0 (the energy is entirely assigned 335 //to the electron) and to disregard fluoresc 356 //to the electron) and to disregard fluorescence. 336 if (shellIndex == 9) 357 if (shellIndex == 9) 337 bindingEnergy = 0.*eV; 358 bindingEnergy = 0.*eV; 338 359 >> 360 339 G4double localEnergyDeposit = 0.0; 361 G4double localEnergyDeposit = 0.0; 340 G4double cosTheta = 1.0; 362 G4double cosTheta = 1.0; 341 363 342 // Primary outcoming electron 364 // Primary outcoming electron 343 G4double eKineticEnergy = photonEnergy - bin 365 G4double eKineticEnergy = photonEnergy - bindingEnergy; 344 << 366 345 // There may be cases where the binding ener 367 // There may be cases where the binding energy of the selected shell is > photon energy 346 // In such cases do not generate secondaries 368 // In such cases do not generate secondaries 347 if (eKineticEnergy > 0.) 369 if (eKineticEnergy > 0.) 348 { 370 { 349 // The electron is created 371 // The electron is created 350 // Direction sampled from the Sauter dis 372 // Direction sampled from the Sauter distribution 351 cosTheta = SampleElectronDirection(eKine 373 cosTheta = SampleElectronDirection(eKineticEnergy); 352 G4double sinTheta = std::sqrt(1-cosTheta 374 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); 353 G4double phi = twopi * G4UniformRand() ; 375 G4double phi = twopi * G4UniformRand() ; 354 G4double dirx = sinTheta * std::cos(phi) 376 G4double dirx = sinTheta * std::cos(phi); 355 G4double diry = sinTheta * std::sin(phi) 377 G4double diry = sinTheta * std::sin(phi); 356 G4double dirz = cosTheta ; 378 G4double dirz = cosTheta ; 357 G4ThreeVector electronDirection(dirx,dir 379 G4ThreeVector electronDirection(dirx,diry,dirz); //electron direction 358 electronDirection.rotateUz(photonDirecti 380 electronDirection.rotateUz(photonDirection); 359 G4DynamicParticle* electron = new G4Dyna << 381 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 360 electronDirection, << 382 electronDirection, 361 eKineticEnergy); 383 eKineticEnergy); 362 fvect->push_back(electron); 384 fvect->push_back(electron); 363 } << 385 } 364 else << 386 else 365 bindingEnergy = photonEnergy; 387 bindingEnergy = photonEnergy; 366 388 >> 389 367 G4double energyInFluorescence = 0; //testing 390 G4double energyInFluorescence = 0; //testing purposes 368 G4double energyInAuger = 0; //testing purpos 391 G4double energyInAuger = 0; //testing purposes 369 << 392 370 //Now, take care of fluorescence, if require 393 //Now, take care of fluorescence, if required. According to the Penelope 371 //recipe, I have to skip fluoresence complet << 394 //recipe, I have to skip fluoresence completely if shellIndex == 9 372 //(= sampling of a shell outer than K,L,M) 395 //(= sampling of a shell outer than K,L,M) 373 if (fAtomDeexcitation && shellIndex<9) 396 if (fAtomDeexcitation && shellIndex<9) 374 { << 397 { 375 G4int index = couple->GetIndex(); 398 G4int index = couple->GetIndex(); 376 if (fAtomDeexcitation->CheckDeexcitation 399 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) 377 { << 400 { 378 std::size_t nBefore = fvect->size(); << 401 size_t nBefore = fvect->size(); 379 fAtomDeexcitation->GenerateParticles(fvect 402 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index); 380 std::size_t nAfter = fvect->size(); << 403 size_t nAfter = fvect->size(); 381 404 382 if (nAfter > nBefore) //actual production 405 if (nAfter > nBefore) //actual production of fluorescence 383 { 406 { 384 for (std::size_t j=nBefore;j<nAfter;++ << 407 for (size_t j=nBefore;j<nAfter;j++) //loop on products 385 { 408 { 386 G4double itsEnergy = ((*fvect)[j])->GetK 409 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy(); 387 if (itsEnergy < bindingEnergy) // valid << 410 bindingEnergy -= itsEnergy; 388 { << 411 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition()) 389 bindingEnergy -= itsEnergy; << 412 energyInFluorescence += itsEnergy; 390 if (((*fvect)[j])->GetParticleDefini << 413 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition()) 391 energyInFluorescence += itsEnergy; << 414 energyInAuger += itsEnergy; 392 else if (((*fvect)[j])->GetParticleD << 393 energyInAuger += itsEnergy; << 394 } << 395 else //invalid secondary: takes more tha << 396 { << 397 delete (*fvect)[j]; << 398 (*fvect)[j] = nullptr; << 399 } << 400 } 415 } 401 } 416 } 402 } 417 } 403 } 418 } 404 419 405 //Residual energy is deposited locally 420 //Residual energy is deposited locally 406 localEnergyDeposit += bindingEnergy; 421 localEnergyDeposit += bindingEnergy; 407 << 422 408 if (localEnergyDeposit < 0) //Should not be: << 423 if (localEnergyDeposit < 0) 409 { 424 { 410 G4Exception("G4PenelopePhotoElectricMode << 425 G4cout << "WARNING - " 411 "em2099",JustWarning,"WARNING: Negative << 426 << "G4PenelopePhotoElectricModel::SampleSecondaries() - Negative energy deposit" >> 427 << G4endl; 412 localEnergyDeposit = 0; 428 localEnergyDeposit = 0; 413 } 429 } 414 430 415 fParticleChange->ProposeLocalEnergyDeposit(l 431 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); 416 432 417 if (fVerboseLevel > 1) << 433 if (verboseLevel > 1) 418 { 434 { 419 G4cout << "----------------------------- 435 G4cout << "-----------------------------------------------------------" << G4endl; 420 G4cout << "Energy balance from G4Penelop 436 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl; 421 G4cout << "Selected shell: " << WriteTar << 437 G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " << 422 anElement->GetName() << G4endl; 438 anElement->GetName() << G4endl; 423 G4cout << "Incoming photon energy: " << 439 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl; 424 G4cout << "----------------------------- 440 G4cout << "-----------------------------------------------------------" << G4endl; 425 if (eKineticEnergy) 441 if (eKineticEnergy) 426 G4cout << "Outgoing electron " << eKineticEn 442 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl; 427 if (energyInFluorescence) 443 if (energyInFluorescence) 428 G4cout << "Fluorescence x-rays: " << energyI 444 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl; 429 if (energyInAuger) 445 if (energyInAuger) 430 G4cout << "Auger electrons: " << energyInAug 446 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl; 431 G4cout << "Local energy deposit " << loc 447 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 432 G4cout << "Total final state: " << << 448 G4cout << "Total final state: " << 433 (eKineticEnergy+energyInFluorescence+localEn << 449 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV << 434 " keV" << G4endl; 450 " keV" << G4endl; 435 G4cout << "----------------------------- 451 G4cout << "-----------------------------------------------------------" << G4endl; 436 } 452 } 437 if (fVerboseLevel > 0) << 453 if (verboseLevel > 0) 438 { 454 { 439 G4double energyDiff = << 455 G4double energyDiff = 440 std::fabs(eKineticEnergy+energyInFluorescenc 456 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy); 441 if (energyDiff > 0.05*keV) 457 if (energyDiff > 0.05*keV) 442 { 458 { 443 G4cout << "Warning from G4PenelopePhotoEle << 459 G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " << 444 (eKineticEnergy+energyInFluorescence+loc << 460 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV 445 << " keV (final) vs. " << << 461 << " keV (final) vs. " << 446 photonEnergy/keV << " keV (initial)" << 462 photonEnergy/keV << " keV (initial)" << G4endl; 447 G4cout << "------------------------------- 463 G4cout << "-----------------------------------------------------------" << G4endl; 448 G4cout << "Energy balance from G4PenelopeP 464 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl; 449 G4cout << "Selected shell: " << WriteTarge << 465 G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " << 450 anElement->GetName() << G4endl; 466 anElement->GetName() << G4endl; 451 G4cout << "Incoming photon energy: " << ph 467 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl; 452 G4cout << "------------------------------- 468 G4cout << "-----------------------------------------------------------" << G4endl; 453 if (eKineticEnergy) 469 if (eKineticEnergy) 454 G4cout << "Outgoing electron " << eKinet 470 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl; 455 if (energyInFluorescence) 471 if (energyInFluorescence) 456 G4cout << "Fluorescence x-rays: " << ene 472 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl; 457 if (energyInAuger) 473 if (energyInAuger) 458 G4cout << "Auger electrons: " << energyI 474 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl; 459 G4cout << "Local energy deposit " << local 475 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 460 G4cout << "Total final state: " << << 476 G4cout << "Total final state: " << 461 (eKineticEnergy+energyInFluorescence+loc << 477 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV << 462 " keV" << G4endl; 478 " keV" << G4endl; 463 G4cout << "------------------------------- 479 G4cout << "-----------------------------------------------------------" << G4endl; 464 } 480 } 465 } 481 } 466 } 482 } 467 483 468 //....oooOO0OOooo........oooOO0OOooo........oo 484 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 469 485 470 G4double G4PenelopePhotoElectricModel::SampleE 486 G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy) 471 { 487 { 472 G4double costheta = 1.0; 488 G4double costheta = 1.0; 473 if (energy>1*GeV) return costheta; 489 if (energy>1*GeV) return costheta; 474 << 490 475 //1) initialize energy-dependent variables 491 //1) initialize energy-dependent variables 476 // Variable naming according to Eq. (2.24) o 492 // Variable naming according to Eq. (2.24) of Penelope Manual 477 // (pag. 44) 493 // (pag. 44) 478 G4double gamma = 1.0 + energy/electron_mass_ 494 G4double gamma = 1.0 + energy/electron_mass_c2; 479 G4double gamma2 = gamma*gamma; 495 G4double gamma2 = gamma*gamma; 480 G4double beta = std::sqrt((gamma2-1.0)/gamma 496 G4double beta = std::sqrt((gamma2-1.0)/gamma2); 481 << 497 482 // ac corresponds to "A" of Eq. (2.31) 498 // ac corresponds to "A" of Eq. (2.31) 483 // 499 // 484 G4double ac = (1.0/beta) - 1.0; 500 G4double ac = (1.0/beta) - 1.0; 485 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(ga 501 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0); 486 G4double a2 = ac + 2.0; 502 G4double a2 = ac + 2.0; 487 G4double gtmax = 2.0*(a1 + 1.0/ac); 503 G4double gtmax = 2.0*(a1 + 1.0/ac); 488 << 504 489 G4double tsam = 0; 505 G4double tsam = 0; 490 G4double gtr = 0; 506 G4double gtr = 0; 491 507 492 //2) sampling. Eq. (2.31) of Penelope Manual 508 //2) sampling. Eq. (2.31) of Penelope Manual 493 // tsam = 1-std::cos(theta) 509 // tsam = 1-std::cos(theta) 494 // gtr = rejection function according to Eq. 510 // gtr = rejection function according to Eq. (2.28) 495 do{ 511 do{ 496 G4double rand = G4UniformRand(); 512 G4double rand = G4UniformRand(); 497 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(r 513 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand); 498 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam)); 514 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam)); 499 }while(G4UniformRand()*gtmax > gtr); 515 }while(G4UniformRand()*gtmax > gtr); 500 costheta = 1.0-tsam; 516 costheta = 1.0-tsam; >> 517 501 518 502 return costheta; 519 return costheta; 503 } 520 } 504 521 505 //....oooOO0OOooo........oooOO0OOooo........oo 522 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 506 523 507 void G4PenelopePhotoElectricModel::ReadDataFil 524 void G4PenelopePhotoElectricModel::ReadDataFile(G4int Z) 508 { 525 { 509 if (!IsMaster()) << 526 if (!IsMaster()) 510 //Should not be here! 527 //Should not be here! 511 G4Exception("G4PenelopePhotoElectricModel: 528 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()", 512 "em0100",FatalException,"Worker thread in << 529 "em0100",FatalException,"Worker thread in this method"); 513 530 514 if (fVerboseLevel > 2) << 531 if (verboseLevel > 2) 515 { 532 { 516 G4cout << "G4PenelopePhotoElectricModel: 533 G4cout << "G4PenelopePhotoElectricModel::ReadDataFile()" << G4endl; 517 G4cout << "Going to read PhotoElectric d 534 G4cout << "Going to read PhotoElectric data files for Z=" << Z << G4endl; 518 } 535 } 519 << 536 520 const char* path = G4FindDataDir("G4LEDATA << 537 char* path = getenv("G4LEDATA"); 521 if(!path) << 538 if (!path) 522 { 539 { 523 G4String excep = "G4PenelopePhotoElectri 540 G4String excep = "G4PenelopePhotoElectricModel - G4LEDATA environment variable not set!"; 524 G4Exception("G4PenelopePhotoElectricMode 541 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()", 525 "em0006",FatalException,excep); 542 "em0006",FatalException,excep); 526 return; 543 return; 527 } 544 } 528 << 545 529 /* 546 /* 530 Read the cross section file 547 Read the cross section file 531 */ 548 */ 532 std::ostringstream ost; 549 std::ostringstream ost; 533 if (Z>9) 550 if (Z>9) 534 ost << path << "/penelope/photoelectric/pd 551 ost << path << "/penelope/photoelectric/pdgph" << Z << ".p08"; 535 else 552 else 536 ost << path << "/penelope/photoelectric/pd 553 ost << path << "/penelope/photoelectric/pdgph0" << Z << ".p08"; 537 std::ifstream file(ost.str().c_str()); 554 std::ifstream file(ost.str().c_str()); 538 if (!file.is_open()) 555 if (!file.is_open()) 539 { 556 { 540 G4String excep = "G4PenelopePhotoElectri 557 G4String excep = "G4PenelopePhotoElectricModel - data file " + G4String(ost.str()) + " not found!"; 541 G4Exception("G4PenelopePhotoElectricMode 558 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()", 542 "em0003",FatalException,excep); 559 "em0003",FatalException,excep); 543 } 560 } 544 //I have to know in advance how many points 561 //I have to know in advance how many points are in the data list 545 //to initialize the G4PhysicsFreeVector() 562 //to initialize the G4PhysicsFreeVector() 546 std::size_t ndata=0; << 563 size_t ndata=0; 547 G4String line; 564 G4String line; 548 while( getline(file, line) ) 565 while( getline(file, line) ) 549 ndata++; 566 ndata++; 550 ndata -= 1; 567 ndata -= 1; 551 //G4cout << "Found: " << ndata << " lines" < 568 //G4cout << "Found: " << ndata << " lines" << G4endl; 552 569 553 file.clear(); 570 file.clear(); 554 file.close(); 571 file.close(); 555 file.open(ost.str().c_str()); 572 file.open(ost.str().c_str()); 556 573 557 G4int readZ =0; 574 G4int readZ =0; 558 std::size_t nShells= 0; << 575 size_t nShells= 0; 559 file >> readZ >> nShells; 576 file >> readZ >> nShells; 560 577 561 if (fVerboseLevel > 3) << 578 if (verboseLevel > 3) 562 G4cout << "Element Z=" << Z << " , nShells 579 G4cout << "Element Z=" << Z << " , nShells = " << nShells << G4endl; 563 580 564 //check the right file is opened. 581 //check the right file is opened. 565 if (readZ != Z || nShells <= 0 || nShells > 582 if (readZ != Z || nShells <= 0 || nShells > 50) //protect nShell against large values 566 { 583 { 567 G4ExceptionDescription ed; 584 G4ExceptionDescription ed; 568 ed << "Corrupted data file for Z=" << Z 585 ed << "Corrupted data file for Z=" << Z << G4endl; 569 G4Exception("G4PenelopePhotoElectricMode 586 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()", 570 "em0005",FatalException,ed); 587 "em0005",FatalException,ed); 571 return; 588 return; 572 } 589 } 573 G4PhysicsTable* thePhysicsTable = new G4Phys 590 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable(); 574 << 591 575 //the table has to contain nShell+1 G4Physic << 592 //the table has to contain nShell+1 G4PhysicsFreeVectors, 576 //(theTable)[0] --> total cross section 593 //(theTable)[0] --> total cross section 577 //(theTable)[ishell] --> cross section for s 594 //(theTable)[ishell] --> cross section for shell (ishell-1) 578 595 579 //reserve space for the vectors 596 //reserve space for the vectors 580 //everything is log-log 597 //everything is log-log 581 for (std::size_t i=0;i<nShells+1;++i) << 598 for (size_t i=0;i<nShells+1;i++) 582 thePhysicsTable->push_back(new G4PhysicsFr 599 thePhysicsTable->push_back(new G4PhysicsFreeVector(ndata)); 583 600 584 std::size_t k =0; << 601 size_t k =0; 585 for (k=0;k<ndata && !file.eof();++k) << 602 for (k=0;k<ndata && !file.eof();k++) 586 { 603 { 587 G4double energy = 0; 604 G4double energy = 0; 588 G4double aValue = 0; 605 G4double aValue = 0; 589 file >> energy ; 606 file >> energy ; 590 energy *= eV; 607 energy *= eV; 591 G4double logene = G4Log(energy); << 608 G4double logene = std::log(energy); 592 //loop on the columns 609 //loop on the columns 593 for (std::size_t i=0;i<nShells+1;++i) << 610 for (size_t i=0;i<nShells+1;i++) 594 { 611 { 595 file >> aValue; 612 file >> aValue; 596 aValue *= barn; 613 aValue *= barn; 597 G4PhysicsFreeVector* theVec = (G4PhysicsFr << 614 G4PhysicsFreeVector* theVec = (G4PhysicsFreeVector*) ((*thePhysicsTable)[i]); 598 if (aValue < 1e-40*cm2) //protection again 615 if (aValue < 1e-40*cm2) //protection against log(0) 599 aValue = 1e-40*cm2; 616 aValue = 1e-40*cm2; 600 theVec->PutValue(k,logene,G4Log(aValue)); << 617 theVec->PutValue(k,logene,std::log(aValue)); 601 } 618 } 602 } 619 } 603 620 604 if (fVerboseLevel > 2) << 621 if (verboseLevel > 2) 605 { 622 { 606 G4cout << "G4PenelopePhotoElectricModel: << 623 G4cout << "G4PenelopePhotoElectricModel: read " << k << " points for element Z = " 607 << Z << G4endl; 624 << Z << G4endl; 608 } 625 } 609 626 610 fLogAtomicShellXS[Z] = thePhysicsTable; << 627 logAtomicShellXS->insert(std::make_pair(Z,thePhysicsTable)); 611 << 628 612 file.close(); 629 file.close(); 613 return; 630 return; 614 } 631 } 615 632 616 //....oooOO0OOooo........oooOO0OOooo........oo 633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 617 634 618 std::size_t G4PenelopePhotoElectricModel::GetN << 635 size_t G4PenelopePhotoElectricModel::GetNumberOfShellXS(G4int Z) 619 { 636 { 620 if (!IsMaster()) 637 if (!IsMaster()) 621 //Should not be here! 638 //Should not be here! 622 G4Exception("G4PenelopePhotoElectricModel: 639 G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()", 623 "em0100",FatalException,"Worker thread in << 640 "em0100",FatalException,"Worker thread in this method"); 624 641 625 //read data files 642 //read data files 626 if (!fLogAtomicShellXS[Z]) << 643 if (!logAtomicShellXS->count(Z)) 627 ReadDataFile(Z); 644 ReadDataFile(Z); 628 //now it should be ok 645 //now it should be ok 629 if (!fLogAtomicShellXS[Z]) << 646 if (!logAtomicShellXS->count(Z)) 630 { 647 { 631 G4ExceptionDescription ed; 648 G4ExceptionDescription ed; 632 ed << "Cannot find shell cross section 649 ed << "Cannot find shell cross section data for Z=" << Z << G4endl; 633 G4Exception("G4PenelopePhotoElectricMod 650 G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()", 634 "em2038",FatalException,ed); 651 "em2038",FatalException,ed); 635 } 652 } 636 //one vector is allocated for the _total_ cr 653 //one vector is allocated for the _total_ cross section 637 std::size_t nEntries = fLogAtomicShellXS[Z]- << 654 size_t nEntries = logAtomicShellXS->find(Z)->second->entries(); 638 return (nEntries-1); 655 return (nEntries-1); 639 } 656 } 640 657 641 //....oooOO0OOooo........oooOO0OOooo........oo 658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 642 659 643 G4double G4PenelopePhotoElectricModel::GetShel << 660 G4double G4PenelopePhotoElectricModel::GetShellCrossSection(G4int Z,size_t shellID,G4double energy) 644 { 661 { 645 //this forces also the loading of the data 662 //this forces also the loading of the data 646 std::size_t entries = GetNumberOfShellXS(Z); << 663 size_t entries = GetNumberOfShellXS(Z); 647 664 648 if (shellID >= entries) 665 if (shellID >= entries) 649 { 666 { 650 G4cout << "Element Z=" << Z << " has dat 667 G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl; 651 G4cout << "so shellID should be from 0 t 668 G4cout << "so shellID should be from 0 to " << entries-1 << G4endl; 652 return 0; 669 return 0; 653 } 670 } 654 << 671 655 G4PhysicsTable* theTable = fLogAtomicShellX << 672 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second; 656 //[0] is the total XS, shellID is in the ele 673 //[0] is the total XS, shellID is in the element [shellID+1] 657 G4PhysicsFreeVector* totalXSLog = (G4Physics 674 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1]; 658 << 675 659 if (!totalXSLog) 676 if (!totalXSLog) 660 { 677 { 661 G4Exception("G4PenelopePhotoElectricMod 678 G4Exception("G4PenelopePhotoElectricModel::GetShellCrossSection()", 662 "em2039",FatalException, 679 "em2039",FatalException, 663 "Unable to retrieve the total cross sec 680 "Unable to retrieve the total cross section table"); 664 return 0; 681 return 0; 665 } 682 } 666 G4double logene = G4Log(energy); << 683 G4double logene = std::log(energy); 667 G4double logXS = totalXSLog->Value(logene); 684 G4double logXS = totalXSLog->Value(logene); 668 G4double cross = G4Exp(logXS); << 685 G4double cross = std::exp(logXS); 669 if (cross < 2e-40*cm2) cross = 0; 686 if (cross < 2e-40*cm2) cross = 0; 670 return cross; 687 return cross; 671 } 688 } 672 689 673 //....oooOO0OOooo........oooOO0OOooo........oo 690 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 674 691 675 G4String G4PenelopePhotoElectricModel::WriteTa << 692 G4String G4PenelopePhotoElectricModel::WriteTargetShell(size_t shellID) 676 { 693 { 677 G4String theShell = "outer shell"; 694 G4String theShell = "outer shell"; 678 if (shellID == 0) 695 if (shellID == 0) 679 theShell = "K"; 696 theShell = "K"; 680 else if (shellID == 1) 697 else if (shellID == 1) 681 theShell = "L1"; 698 theShell = "L1"; 682 else if (shellID == 2) 699 else if (shellID == 2) 683 theShell = "L2"; 700 theShell = "L2"; 684 else if (shellID == 3) 701 else if (shellID == 3) 685 theShell = "L3"; 702 theShell = "L3"; 686 else if (shellID == 4) 703 else if (shellID == 4) 687 theShell = "M1"; 704 theShell = "M1"; 688 else if (shellID == 5) 705 else if (shellID == 5) 689 theShell = "M2"; 706 theShell = "M2"; 690 else if (shellID == 6) 707 else if (shellID == 6) 691 theShell = "M3"; 708 theShell = "M3"; 692 else if (shellID == 7) 709 else if (shellID == 7) 693 theShell = "M4"; 710 theShell = "M4"; 694 else if (shellID == 8) 711 else if (shellID == 8) 695 theShell = "M5"; 712 theShell = "M5"; 696 << 713 697 return theShell; 714 return theShell; 698 } 715 } 699 716 700 //....oooOO0OOooo........oooOO0OOooo........oo 717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 701 718 702 void G4PenelopePhotoElectricModel::SetParticle 719 void G4PenelopePhotoElectricModel::SetParticle(const G4ParticleDefinition* p) 703 { 720 { 704 if(!fParticle) { 721 if(!fParticle) { 705 fParticle = p; << 722 fParticle = p; 706 } 723 } 707 } 724 } 708 725 709 //....oooOO0OOooo........oooOO0OOooo........oo 726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 710 727 711 std::size_t G4PenelopePhotoElectricModel::Sele << 728 size_t G4PenelopePhotoElectricModel::SelectRandomShell(G4int Z,G4double energy) 712 { 729 { 713 G4double logEnergy = G4Log(energy); << 730 G4double logEnergy = std::log(energy); 714 731 715 //Check if data have been read (it should be 732 //Check if data have been read (it should be!) 716 if (!fLogAtomicShellXS[Z]) << 733 if (!logAtomicShellXS->count(Z)) 717 { 734 { 718 G4ExceptionDescription ed; 735 G4ExceptionDescription ed; 719 ed << "Cannot find shell cross section 736 ed << "Cannot find shell cross section data for Z=" << Z << G4endl; 720 G4Exception("G4PenelopePhotoElectricMod 737 G4Exception("G4PenelopePhotoElectricModel::SelectRandomShell()", 721 "em2038",FatalException,ed); 738 "em2038",FatalException,ed); 722 } 739 } 723 740 724 G4PhysicsTable* theTable = fLogAtomicShellX << 741 // size_t shellIndex = 0; >> 742 >> 743 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second; 725 744 726 G4double sum = 0; 745 G4double sum = 0; >> 746 727 G4PhysicsFreeVector* totalXSLog = (G4Physics 747 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0]; 728 G4double logXS = totalXSLog->Value(logEnergy 748 G4double logXS = totalXSLog->Value(logEnergy); 729 G4double totalXS = G4Exp(logXS); << 749 G4double totalXS = std::exp(logXS); 730 << 750 731 //Notice: totalXS is the total cross section << 751 //Notice: totalXS is the total cross section and it does *not* correspond to 732 //the sum of partialXS's, since these includ 752 //the sum of partialXS's, since these include only K, L and M shells. 733 // 753 // 734 // Therefore, here one have to consider the << 754 // Therefore, here one have to consider the possibility of ionisation of 735 // an outer shell. Conventionally, it is ind 755 // an outer shell. Conventionally, it is indicated with id=10 in Penelope 736 // 756 // 737 G4double random = G4UniformRand()*totalXS; << 757 G4double random = G4UniformRand()*totalXS; 738 758 739 for (std::size_t k=1;k<theTable->entries();+ << 759 for (size_t k=1;k<theTable->entries();k++) 740 { 760 { 741 //Add one shell 761 //Add one shell 742 G4PhysicsFreeVector* partialXSLog = (G4P 762 G4PhysicsFreeVector* partialXSLog = (G4PhysicsFreeVector*) (*theTable)[k]; 743 G4double logXSLocal = partialXSLog->Valu 763 G4double logXSLocal = partialXSLog->Value(logEnergy); 744 G4double partialXS = G4Exp(logXSLocal); << 764 G4double partialXS = std::exp(logXSLocal); 745 sum += partialXS; 765 sum += partialXS; 746 if (random <= sum) << 766 if (random <= sum) 747 return k-1; << 767 return k-1; 748 } 768 } 749 //none of the shells K, L, M: return outer s 769 //none of the shells K, L, M: return outer shell 750 return 9; 770 return 9; 751 } 771 } 752 772