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: G4PenelopeGammaConversionModel.cc 99415 2016-09-21 09:05:43Z gcosmo $ 26 // 27 // 27 // Author: Luciano Pandola 28 // Author: Luciano Pandola 28 // 29 // 29 // History: 30 // History: 30 // -------- 31 // -------- 31 // 13 Jan 2010 L Pandola First implementa 32 // 13 Jan 2010 L Pandola First implementation (updated to Penelope08) 32 // 24 May 2011 L Pandola Renamed (make v2 33 // 24 May 2011 L Pandola Renamed (make v2008 as default Penelope) 33 // 18 Sep 2013 L Pandola Migration to MT 34 // 18 Sep 2013 L Pandola Migration to MT paradigm. Only master model deals with 34 // data and create 35 // data and creates shared tables 35 // 36 // 36 37 37 #include "G4PenelopeGammaConversionModel.hh" 38 #include "G4PenelopeGammaConversionModel.hh" 38 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh" 39 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh" 40 #include "G4ParticleDefinition.hh" 41 #include "G4ParticleDefinition.hh" 41 #include "G4MaterialCutsCouple.hh" 42 #include "G4MaterialCutsCouple.hh" 42 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh" 43 #include "G4DynamicParticle.hh" 44 #include "G4DynamicParticle.hh" 44 #include "G4Element.hh" 45 #include "G4Element.hh" 45 #include "G4Gamma.hh" 46 #include "G4Gamma.hh" 46 #include "G4Electron.hh" 47 #include "G4Electron.hh" 47 #include "G4Positron.hh" 48 #include "G4Positron.hh" 48 #include "G4PhysicsFreeVector.hh" 49 #include "G4PhysicsFreeVector.hh" 49 #include "G4MaterialCutsCouple.hh" 50 #include "G4MaterialCutsCouple.hh" 50 #include "G4AutoLock.hh" 51 #include "G4AutoLock.hh" 51 #include "G4Exp.hh" 52 #include "G4Exp.hh" 52 53 53 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 54 const G4int G4PenelopeGammaConversionModel::fM << 55 G4PhysicsFreeVector* G4PenelopeGammaConversion << 56 G4double G4PenelopeGammaConversionModel::fAtom << 57 1.2281e+02,7.3167e+01,6.9 << 58 6.4696e+01,6.1228e+01,5.7 << 59 5.0787e+01,4.7851e+01,4.6 << 60 4.4503e+01,4.3815e+01,4.3 << 61 4.1586e+01,4.0953e+01,4.0 << 62 3.9756e+01,3.9144e+01,3.8 << 63 3.7174e+01,3.6663e+01,3.5 << 64 3.4688e+01,3.4197e+01,3.3 << 65 3.3068e+01,3.2740e+01,3.2 << 66 3.1884e+01,3.1622e+01,3.1 << 67 3.0950e+01,3.0758e+01,3.0 << 68 3.0097e+01,2.9832e+01,2.9 << 69 2.9247e+01,2.9085e+01,2.8 << 70 2.8580e+01,2.8442e+01,2.8 << 71 2.7973e+01,2.7819e+01,2.7 << 72 2.7285e+01,2.7093e+01,2.6 << 73 2.6516e+01,2.6304e+01,2.6 << 74 2.5730e+01,2.5577e+01,2.5 << 75 2.5100e+01,2.4941e+01,2.4 << 76 2.4506e+01,2.4391e+01,2.4 << 77 2.4039e+01,2.3922e+01,2.3 << 78 2.3621e+01,2.3523e+01,2.3 << 79 2.3238e+01,2.3139e+01,2.3 << 80 2.2833e+01,2.2694e+01,2.2 << 81 2.2446e+01,2.2358e+01,2.2 << 82 << 83 //....oooOO0OOooo........oooOO0OOooo........oo << 84 55 85 G4PenelopeGammaConversionModel::G4PenelopeGamm 56 G4PenelopeGammaConversionModel::G4PenelopeGammaConversionModel(const G4ParticleDefinition* part, 86 const G4String& nam) 57 const G4String& nam) 87 :G4VEmModel(nam),fParticleChange(nullptr),fP << 58 :G4VEmModel(nam),fParticleChange(0),fParticle(0), 88 fEffectiveCharge(nullptr),fMaterialInvScree << 59 logAtomicCrossSection(0), 89 fScreeningFunction(nullptr),fIsInitialised( << 60 fEffectiveCharge(0),fMaterialInvScreeningRadius(0), >> 61 fScreeningFunction(0),isInitialised(false),fLocalTable(false) >> 62 90 { 63 { 91 fIntrinsicLowEnergyLimit = 2.0*electron_mass 64 fIntrinsicLowEnergyLimit = 2.0*electron_mass_c2; 92 fIntrinsicHighEnergyLimit = 100.0*GeV; 65 fIntrinsicHighEnergyLimit = 100.0*GeV; 93 fSmallEnergy = 1.1*MeV; 66 fSmallEnergy = 1.1*MeV; 94 67 >> 68 InitializeScreeningRadii(); >> 69 95 if (part) 70 if (part) 96 SetParticle(part); 71 SetParticle(part); 97 72 98 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 73 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); 99 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 74 SetHighEnergyLimit(fIntrinsicHighEnergyLimit); 100 // 75 // 101 fVerboseLevel= 0; << 76 verboseLevel= 0; 102 // Verbosity scale: 77 // Verbosity scale: 103 // 0 = nothing 78 // 0 = nothing 104 // 1 = warning for energy non-conservation 79 // 1 = warning for energy non-conservation 105 // 2 = details of energy budget 80 // 2 = details of energy budget 106 // 3 = calculation of cross sections, file o 81 // 3 = calculation of cross sections, file openings, sampling of atoms 107 // 4 = entering in methods 82 // 4 = entering in methods 108 } 83 } 109 84 110 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 111 86 112 G4PenelopeGammaConversionModel::~G4PenelopeGam 87 G4PenelopeGammaConversionModel::~G4PenelopeGammaConversionModel() 113 { 88 { 114 //Delete shared tables, they exist only in t 89 //Delete shared tables, they exist only in the master model 115 if (IsMaster() || fLocalTable) 90 if (IsMaster() || fLocalTable) 116 { 91 { 117 for(G4int i=0; i<=fMaxZ; ++i) << 92 if (logAtomicCrossSection) 118 { 93 { 119 if(fLogAtomicCrossSection[i]) { << 94 for (auto& item : (*logAtomicCrossSection)) 120 delete fLogAtomicCrossSection[i]; << 95 if (item.second) delete item.second; 121 fLogAtomicCrossSection[i] = nullptr; << 96 delete logAtomicCrossSection; 122 } << 123 } 97 } 124 if (fEffectiveCharge) 98 if (fEffectiveCharge) 125 delete fEffectiveCharge; 99 delete fEffectiveCharge; 126 if (fMaterialInvScreeningRadius) 100 if (fMaterialInvScreeningRadius) 127 delete fMaterialInvScreeningRadius; 101 delete fMaterialInvScreeningRadius; 128 if (fScreeningFunction) 102 if (fScreeningFunction) 129 delete fScreeningFunction; 103 delete fScreeningFunction; 130 } 104 } 131 } 105 } 132 106 >> 107 133 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 134 109 135 void G4PenelopeGammaConversionModel::Initialis 110 void G4PenelopeGammaConversionModel::Initialise(const G4ParticleDefinition* part, 136 const G4DataVector&) 111 const G4DataVector&) 137 { 112 { 138 if (fVerboseLevel > 3) << 113 if (verboseLevel > 3) 139 G4cout << "Calling G4PenelopeGammaConvers 114 G4cout << "Calling G4PenelopeGammaConversionModel::Initialise()" << G4endl; 140 115 141 SetParticle(part); 116 SetParticle(part); 142 117 143 //Only the master model creates/fills/destro 118 //Only the master model creates/fills/destroys the tables 144 if (IsMaster() && part == fParticle) 119 if (IsMaster() && part == fParticle) 145 { 120 { >> 121 // logAtomicCrossSection is created only once, since it is never cleared >> 122 if (!logAtomicCrossSection) >> 123 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>; >> 124 146 //delete old material data... 125 //delete old material data... 147 if (fEffectiveCharge) 126 if (fEffectiveCharge) 148 { 127 { 149 delete fEffectiveCharge; 128 delete fEffectiveCharge; 150 fEffectiveCharge = nullptr; 129 fEffectiveCharge = nullptr; 151 } 130 } 152 if (fMaterialInvScreeningRadius) 131 if (fMaterialInvScreeningRadius) 153 { 132 { 154 delete fMaterialInvScreeningRadius; 133 delete fMaterialInvScreeningRadius; 155 fMaterialInvScreeningRadius = nullptr; 134 fMaterialInvScreeningRadius = nullptr; 156 } 135 } 157 if (fScreeningFunction) 136 if (fScreeningFunction) 158 { 137 { 159 delete fScreeningFunction; 138 delete fScreeningFunction; 160 fScreeningFunction = nullptr; 139 fScreeningFunction = nullptr; 161 } 140 } 162 //and create new ones 141 //and create new ones 163 fEffectiveCharge = new std::map<const G4 142 fEffectiveCharge = new std::map<const G4Material*,G4double>; 164 fMaterialInvScreeningRadius = new std::m 143 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>; 165 fScreeningFunction = new std::map<const 144 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >; 166 145 167 G4ProductionCutsTable* theCoupleTable = 146 G4ProductionCutsTable* theCoupleTable = 168 G4ProductionCutsTable::GetProductionCutsTabl 147 G4ProductionCutsTable::GetProductionCutsTable(); 169 148 170 for (G4int i=0;i<(G4int)theCoupleTable-> << 149 for (size_t i=0;i<theCoupleTable->GetTableSize();i++) 171 { 150 { 172 const G4Material* material = 151 const G4Material* material = 173 theCoupleTable->GetMaterialCutsCouple(i) 152 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 174 const G4ElementVector* theElementVector = 153 const G4ElementVector* theElementVector = material->GetElementVector(); 175 154 176 for (std::size_t j=0;j<material->GetNumber << 155 for (size_t j=0;j<material->GetNumberOfElements();j++) 177 { 156 { 178 G4int iZ = theElementVector->at(j)->Ge << 157 G4int iZ = (G4int) theElementVector->at(j)->GetZ(); 179 //read data files only in the master 158 //read data files only in the master 180 if (iZ <= fMaxZ && !fLogAtomicCrossSe << 159 if (!logAtomicCrossSection->count(iZ)) 181 ReadDataFile(iZ); 160 ReadDataFile(iZ); 182 } 161 } 183 162 184 //check if material data are available 163 //check if material data are available 185 if (!fEffectiveCharge->count(material)) 164 if (!fEffectiveCharge->count(material)) 186 InitializeScreeningFunctions(material); 165 InitializeScreeningFunctions(material); 187 } 166 } 188 if (fVerboseLevel > 0) { << 167 >> 168 >> 169 if (verboseLevel > 0) { 189 G4cout << "Penelope Gamma Conversion model v 170 G4cout << "Penelope Gamma Conversion model v2008 is initialized " << G4endl 190 << "Energy range: " 171 << "Energy range: " 191 << LowEnergyLimit() / MeV << " MeV - 172 << LowEnergyLimit() / MeV << " MeV - " 192 << HighEnergyLimit() / GeV << " GeV" 173 << HighEnergyLimit() / GeV << " GeV" 193 << G4endl; 174 << G4endl; 194 } 175 } >> 176 195 } 177 } 196 if(fIsInitialised) return; << 178 >> 179 >> 180 if(isInitialised) return; 197 fParticleChange = GetParticleChangeForGamma( 181 fParticleChange = GetParticleChangeForGamma(); 198 fIsInitialised = true; << 182 isInitialised = true; 199 } 183 } 200 184 201 //....oooOO0OOooo........oooOO0OOooo........oo 185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 202 186 203 void G4PenelopeGammaConversionModel::Initialis 187 void G4PenelopeGammaConversionModel::InitialiseLocal(const G4ParticleDefinition* part, 204 G4VEmModel *masterModel) 188 G4VEmModel *masterModel) 205 { 189 { 206 if (fVerboseLevel > 3) << 190 if (verboseLevel > 3) 207 G4cout << "Calling G4PenelopeGammaConvers 191 G4cout << "Calling G4PenelopeGammaConversionModel::InitialiseLocal()" << G4endl; >> 192 208 // 193 // 209 //Check that particle matches: one might hav 194 //Check that particle matches: one might have multiple master models (e.g. 210 //for e+ and e-). 195 //for e+ and e-). 211 // 196 // 212 if (part == fParticle) 197 if (part == fParticle) 213 { 198 { 214 //Get the const table pointers from the 199 //Get the const table pointers from the master to the workers 215 const G4PenelopeGammaConversionModel* th 200 const G4PenelopeGammaConversionModel* theModel = 216 static_cast<G4PenelopeGammaConversionModel*> 201 static_cast<G4PenelopeGammaConversionModel*> (masterModel); 217 202 218 //Copy pointers to the data tables 203 //Copy pointers to the data tables 219 fEffectiveCharge = theModel->fEffectiveC 204 fEffectiveCharge = theModel->fEffectiveCharge; 220 fMaterialInvScreeningRadius = theModel-> 205 fMaterialInvScreeningRadius = theModel->fMaterialInvScreeningRadius; 221 fScreeningFunction = theModel->fScreenin << 206 fScreeningFunction = theModel->fScreeningFunction; 222 for(G4int i=0; i<=fMaxZ; ++i) << 207 logAtomicCrossSection = theModel->logAtomicCrossSection; 223 fLogAtomicCrossSection[i] = theModel->fLogAt << 224 208 225 //Same verbosity for all workers, as the 209 //Same verbosity for all workers, as the master 226 fVerboseLevel = theModel->fVerboseLevel; << 210 verboseLevel = theModel->verboseLevel; 227 } 211 } 228 212 229 return; 213 return; 230 } 214 } 231 215 232 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 233 namespace { G4Mutex PenelopeGammaConversionMo 217 namespace { G4Mutex PenelopeGammaConversionModelMutex = G4MUTEX_INITIALIZER; } 234 218 235 G4double G4PenelopeGammaConversionModel::Compu 219 G4double G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom( 236 const G4ParticleDefinition 220 const G4ParticleDefinition*, 237 G4double energy, 221 G4double energy, 238 G4double Z, G4double, 222 G4double Z, G4double, 239 G4double, G4double) 223 G4double, G4double) 240 { 224 { 241 // 225 // 242 // Penelope model v2008. 226 // Penelope model v2008. 243 // Cross section (including triplet producti 227 // Cross section (including triplet production) read from database and managed 244 // through the G4CrossSectionHandler utility 228 // through the G4CrossSectionHandler utility. Cross section data are from 245 // M.J. Berger and J.H. Hubbel (XCOM), Repor 229 // M.J. Berger and J.H. Hubbel (XCOM), Report NBSIR 887-3598 246 // 230 // 247 231 248 if (energy < fIntrinsicLowEnergyLimit) 232 if (energy < fIntrinsicLowEnergyLimit) 249 return 0; 233 return 0; 250 234 251 G4int iZ = G4int(Z); << 235 G4int iZ = (G4int) Z; 252 236 253 if (!fLogAtomicCrossSection[iZ]) << 237 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was >> 238 //not invoked >> 239 if (!logAtomicCrossSection) >> 240 { >> 241 //create a **thread-local** version of the table. Used only for G4EmCalculator and >> 242 //Unit Tests >> 243 fLocalTable = true; >> 244 logAtomicCrossSection = new std::map<G4int,G4PhysicsFreeVector*>; >> 245 } >> 246 //now it should be ok >> 247 if (!logAtomicCrossSection->count(iZ)) 254 { 248 { 255 //If we are here, it means that Initial 249 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 256 //not filled up. This can happen in a U 250 //not filled up. This can happen in a UnitTest or via G4EmCalculator 257 if (fVerboseLevel > 0) << 251 if (verboseLevel > 0) 258 { 252 { 259 //Issue a G4Exception (warning) only in v 253 //Issue a G4Exception (warning) only in verbose mode 260 G4ExceptionDescription ed; 254 G4ExceptionDescription ed; 261 ed << "Unable to retrieve the cross secti 255 ed << "Unable to retrieve the cross section table for Z=" << iZ << G4endl; 262 ed << "This can happen only in Unit Tests 256 ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl; 263 G4Exception("G4PenelopeGammaConversionMod 257 G4Exception("G4PenelopeGammaConversionModel::ComputeCrossSectionPerAtom()", 264 "em2018",JustWarning,ed); 258 "em2018",JustWarning,ed); 265 } 259 } 266 //protect file reading via autolock 260 //protect file reading via autolock 267 G4AutoLock lock(&PenelopeGammaConversio 261 G4AutoLock lock(&PenelopeGammaConversionModelMutex); 268 ReadDataFile(iZ); 262 ReadDataFile(iZ); 269 lock.unlock(); 263 lock.unlock(); 270 fLocalTable = true; << 271 } 264 } >> 265 272 G4double cs = 0; 266 G4double cs = 0; 273 G4double logene = G4Log(energy); << 267 G4double logene = std::log(energy); 274 G4PhysicsFreeVector* theVec = fLogAtomicCros << 268 >> 269 G4PhysicsFreeVector* theVec = logAtomicCrossSection->find(iZ)->second; >> 270 275 G4double logXS = theVec->Value(logene); 271 G4double logXS = theVec->Value(logene); 276 cs = G4Exp(logXS); 272 cs = G4Exp(logXS); 277 273 278 if (fVerboseLevel > 2) << 274 if (verboseLevel > 2) 279 G4cout << "Gamma conversion cross section 275 G4cout << "Gamma conversion cross section at " << energy/MeV << " MeV for Z=" << Z << 280 " = " << cs/barn << " barn" << G4endl; 276 " = " << cs/barn << " barn" << G4endl; 281 return cs; 277 return cs; 282 } 278 } 283 279 284 //....oooOO0OOooo........oooOO0OOooo........oo 280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 285 281 286 void 282 void 287 G4PenelopeGammaConversionModel::SampleSecondar 283 G4PenelopeGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 288 const G4MaterialCutsCouple* coup 284 const G4MaterialCutsCouple* couple, 289 const G4DynamicParticle* aDynami 285 const G4DynamicParticle* aDynamicGamma, 290 G4double, 286 G4double, 291 G4double) 287 G4double) 292 { 288 { 293 // 289 // 294 // Penelope model v2008. 290 // Penelope model v2008. 295 // Final state is sampled according to the B 291 // Final state is sampled according to the Bethe-Heitler model with Coulomb 296 // corrections, according to the semi-empiri 292 // corrections, according to the semi-empirical model of 297 // J. Baro' et al., Radiat. Phys. Chem. 44 293 // J. Baro' et al., Radiat. Phys. Chem. 44 (1994) 531. 298 // 294 // 299 // The model uses the high energy Coulomb co 295 // The model uses the high energy Coulomb correction from 300 // H. Davies et al., Phys. Rev. 93 (1954) 7 296 // H. Davies et al., Phys. Rev. 93 (1954) 788 301 // and atomic screening radii tabulated from 297 // and atomic screening radii tabulated from 302 // J.H. Hubbel et al., J. Phys. Chem. Ref. 298 // J.H. Hubbel et al., J. Phys. Chem. Ref. Data 9 (1980) 1023 303 // for Z= 1 to 92. 299 // for Z= 1 to 92. 304 // 300 // 305 if (fVerboseLevel > 3) << 301 if (verboseLevel > 3) 306 G4cout << "Calling SamplingSecondaries() o 302 G4cout << "Calling SamplingSecondaries() of G4PenelopeGammaConversionModel" << G4endl; 307 303 308 G4double photonEnergy = aDynamicGamma->GetKi 304 G4double photonEnergy = aDynamicGamma->GetKineticEnergy(); 309 305 310 // Always kill primary 306 // Always kill primary 311 fParticleChange->ProposeTrackStatus(fStopAnd 307 fParticleChange->ProposeTrackStatus(fStopAndKill); 312 fParticleChange->SetProposedKineticEnergy(0. 308 fParticleChange->SetProposedKineticEnergy(0.); 313 309 314 if (photonEnergy <= fIntrinsicLowEnergyLimit 310 if (photonEnergy <= fIntrinsicLowEnergyLimit) 315 { 311 { 316 fParticleChange->ProposeLocalEnergyDepos 312 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy); 317 return ; 313 return ; 318 } 314 } 319 315 320 G4ParticleMomentum photonDirection = aDynami 316 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection(); 321 const G4Material* mat = couple->GetMaterial( 317 const G4Material* mat = couple->GetMaterial(); 322 318 323 //Either Initialize() was not called, or we 319 //Either Initialize() was not called, or we are in a slave and InitializeLocal() was 324 //not invoked 320 //not invoked 325 if (!fEffectiveCharge) 321 if (!fEffectiveCharge) 326 { 322 { 327 //create a **thread-local** version of t 323 //create a **thread-local** version of the table. Used only for G4EmCalculator and 328 //Unit Tests 324 //Unit Tests 329 fLocalTable = true; 325 fLocalTable = true; 330 fEffectiveCharge = new std::map<const G4 326 fEffectiveCharge = new std::map<const G4Material*,G4double>; 331 fMaterialInvScreeningRadius = new std::m 327 fMaterialInvScreeningRadius = new std::map<const G4Material*,G4double>; 332 fScreeningFunction = new std::map<const 328 fScreeningFunction = new std::map<const G4Material*,std::pair<G4double,G4double> >; 333 } 329 } 334 330 335 if (!fEffectiveCharge->count(mat)) 331 if (!fEffectiveCharge->count(mat)) 336 { 332 { 337 //If we are here, it means that Initiali 333 //If we are here, it means that Initialize() was inkoved, but the MaterialTable was 338 //not filled up. This can happen in a Un 334 //not filled up. This can happen in a UnitTest or via G4EmCalculator 339 if (fVerboseLevel > 0) << 335 if (verboseLevel > 0) 340 { 336 { 341 //Issue a G4Exception (warning) only in ve 337 //Issue a G4Exception (warning) only in verbose mode 342 G4ExceptionDescription ed; 338 G4ExceptionDescription ed; 343 ed << "Unable to allocate the EffectiveCha 339 ed << "Unable to allocate the EffectiveCharge data for " << 344 mat->GetName() << G4endl; 340 mat->GetName() << G4endl; 345 ed << "This can happen only in Unit Tests" 341 ed << "This can happen only in Unit Tests" << G4endl; 346 G4Exception("G4PenelopeGammaConversionMode 342 G4Exception("G4PenelopeGammaConversionModel::SampleSecondaries()", 347 "em2019",JustWarning,ed); 343 "em2019",JustWarning,ed); 348 } 344 } 349 //protect file reading via autolock 345 //protect file reading via autolock 350 G4AutoLock lock(&PenelopeGammaConversion 346 G4AutoLock lock(&PenelopeGammaConversionModelMutex); 351 InitializeScreeningFunctions(mat); 347 InitializeScreeningFunctions(mat); 352 lock.unlock(); 348 lock.unlock(); 353 } 349 } 354 350 355 // eps is the fraction of the photon energy 351 // eps is the fraction of the photon energy assigned to e- (including rest mass) 356 G4double eps = 0; 352 G4double eps = 0; 357 G4double eki = electron_mass_c2/photonEnergy 353 G4double eki = electron_mass_c2/photonEnergy; 358 354 359 //Do it fast for photon energy < 1.1 MeV (cl 355 //Do it fast for photon energy < 1.1 MeV (close to threshold) 360 if (photonEnergy < fSmallEnergy) 356 if (photonEnergy < fSmallEnergy) 361 eps = eki + (1.0-2.0*eki)*G4UniformRand(); 357 eps = eki + (1.0-2.0*eki)*G4UniformRand(); 362 else 358 else 363 { 359 { 364 //Complete calculation 360 //Complete calculation 365 G4double effC = fEffectiveCharge->find(m 361 G4double effC = fEffectiveCharge->find(mat)->second; 366 G4double alz = effC*fine_structure_const 362 G4double alz = effC*fine_structure_const; 367 G4double T = std::sqrt(2.0*eki); 363 G4double T = std::sqrt(2.0*eki); 368 G4double F00=(-1.774-1.210e1*alz+1.118e1 364 G4double F00=(-1.774-1.210e1*alz+1.118e1*alz*alz)*T 369 +(8.523+7.326e1*alz-4.441e1*alz*alz)* 365 +(8.523+7.326e1*alz-4.441e1*alz*alz)*T*T 370 -(1.352e1+1.211e2*alz-9.641e1*alz*alz 366 -(1.352e1+1.211e2*alz-9.641e1*alz*alz)*T*T*T 371 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T 367 +(8.946+6.205e1*alz-6.341e1*alz*alz)*T*T*T*T; 372 368 373 G4double F0b = fScreeningFunction->find( 369 G4double F0b = fScreeningFunction->find(mat)->second.second; 374 G4double g0 = F0b + F00; 370 G4double g0 = F0b + F00; 375 G4double invRad = fMaterialInvScreeningR 371 G4double invRad = fMaterialInvScreeningRadius->find(mat)->second; 376 G4double bmin = 4.0*eki/invRad; 372 G4double bmin = 4.0*eki/invRad; 377 std::pair<G4double,G4double> scree = Ge 373 std::pair<G4double,G4double> scree = GetScreeningFunctions(bmin); 378 G4double g1 = scree.first; 374 G4double g1 = scree.first; 379 G4double g2 = scree.second; 375 G4double g2 = scree.second; 380 G4double g1min = g1+g0; 376 G4double g1min = g1+g0; 381 G4double g2min = g2+g0; 377 G4double g2min = g2+g0; 382 G4double xr = 0.5-eki; 378 G4double xr = 0.5-eki; 383 G4double a1 = 2.*g1min*xr*xr/3.; 379 G4double a1 = 2.*g1min*xr*xr/3.; 384 G4double p1 = a1/(a1+g2min); 380 G4double p1 = a1/(a1+g2min); 385 381 386 G4bool loopAgain = false; 382 G4bool loopAgain = false; 387 //Random sampling of eps 383 //Random sampling of eps 388 do{ 384 do{ 389 loopAgain = false; 385 loopAgain = false; 390 if (G4UniformRand() <= p1) 386 if (G4UniformRand() <= p1) 391 { 387 { 392 G4double ru2m1 = 2.0*G4UniformRand()-1. 388 G4double ru2m1 = 2.0*G4UniformRand()-1.0; 393 if (ru2m1 < 0) 389 if (ru2m1 < 0) 394 eps = 0.5-xr*std::pow(std::abs(ru2m1), 390 eps = 0.5-xr*std::pow(std::abs(ru2m1),1./3.); 395 else 391 else 396 eps = 0.5+xr*std::pow(ru2m1,1./3.); 392 eps = 0.5+xr*std::pow(ru2m1,1./3.); 397 G4double B = eki/(invRad*eps*(1.0-eps)); 393 G4double B = eki/(invRad*eps*(1.0-eps)); 398 scree = GetScreeningFunctions(B); 394 scree = GetScreeningFunctions(B); 399 g1 = scree.first; 395 g1 = scree.first; 400 g1 = std::max(g1+g0,0.); 396 g1 = std::max(g1+g0,0.); 401 if (G4UniformRand()*g1min > g1) 397 if (G4UniformRand()*g1min > g1) 402 loopAgain = true; 398 loopAgain = true; 403 } 399 } 404 else 400 else 405 { 401 { 406 eps = eki+2.0*xr*G4UniformRand(); 402 eps = eki+2.0*xr*G4UniformRand(); 407 G4double B = eki/(invRad*eps*(1.0-eps)); 403 G4double B = eki/(invRad*eps*(1.0-eps)); 408 scree = GetScreeningFunctions(B); 404 scree = GetScreeningFunctions(B); 409 g2 = scree.second; 405 g2 = scree.second; 410 g2 = std::max(g2+g0,0.); 406 g2 = std::max(g2+g0,0.); 411 if (G4UniformRand()*g2min > g2) 407 if (G4UniformRand()*g2min > g2) 412 loopAgain = true; 408 loopAgain = true; 413 } 409 } 414 }while(loopAgain); 410 }while(loopAgain); >> 411 415 } 412 } 416 if (fVerboseLevel > 4) << 413 if (verboseLevel > 4) 417 G4cout << "Sampled eps = " << eps << G4end 414 G4cout << "Sampled eps = " << eps << G4endl; 418 415 419 G4double electronTotEnergy = eps*photonEnerg 416 G4double electronTotEnergy = eps*photonEnergy; 420 G4double positronTotEnergy = (1.0-eps)*photo 417 G4double positronTotEnergy = (1.0-eps)*photonEnergy; 421 418 422 // Scattered electron (positron) angles. ( Z 419 // Scattered electron (positron) angles. ( Z - axis along the parent photon) 423 420 424 //electron kinematics 421 //electron kinematics 425 G4double electronKineEnergy = std::max(0.,el 422 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ; 426 G4double costheta_el = G4UniformRand()*2.0-1 423 G4double costheta_el = G4UniformRand()*2.0-1.0; 427 G4double kk = std::sqrt(electronKineEnergy*( 424 G4double kk = std::sqrt(electronKineEnergy*(electronKineEnergy+2.*electron_mass_c2)); 428 costheta_el = (costheta_el*electronTotEnergy 425 costheta_el = (costheta_el*electronTotEnergy+kk)/(electronTotEnergy+costheta_el*kk); 429 G4double phi_el = twopi * G4UniformRand() ; 426 G4double phi_el = twopi * G4UniformRand() ; 430 G4double dirX_el = std::sqrt(1.-costheta_el* 427 G4double dirX_el = std::sqrt(1.-costheta_el*costheta_el) * std::cos(phi_el); 431 G4double dirY_el = std::sqrt(1.-costheta_el* 428 G4double dirY_el = std::sqrt(1.-costheta_el*costheta_el) * std::sin(phi_el); 432 G4double dirZ_el = costheta_el; 429 G4double dirZ_el = costheta_el; 433 430 434 //positron kinematics 431 //positron kinematics 435 G4double positronKineEnergy = std::max(0.,po 432 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ; 436 G4double costheta_po = G4UniformRand()*2.0-1 433 G4double costheta_po = G4UniformRand()*2.0-1.0; 437 kk = std::sqrt(positronKineEnergy*(positronK 434 kk = std::sqrt(positronKineEnergy*(positronKineEnergy+2.*electron_mass_c2)); 438 costheta_po = (costheta_po*positronTotEnergy 435 costheta_po = (costheta_po*positronTotEnergy+kk)/(positronTotEnergy+costheta_po*kk); 439 G4double phi_po = twopi * G4UniformRand() ; 436 G4double phi_po = twopi * G4UniformRand() ; 440 G4double dirX_po = std::sqrt(1.-costheta_po* 437 G4double dirX_po = std::sqrt(1.-costheta_po*costheta_po) * std::cos(phi_po); 441 G4double dirY_po = std::sqrt(1.-costheta_po* 438 G4double dirY_po = std::sqrt(1.-costheta_po*costheta_po) * std::sin(phi_po); 442 G4double dirZ_po = costheta_po; 439 G4double dirZ_po = costheta_po; 443 440 444 // Kinematics of the created pair: 441 // Kinematics of the created pair: 445 // the electron and positron are assumed to 442 // the electron and positron are assumed to have a symetric angular 446 // distribution with respect to the Z axis a 443 // distribution with respect to the Z axis along the parent photon 447 G4double localEnergyDeposit = 0. ; 444 G4double localEnergyDeposit = 0. ; 448 445 449 if (electronKineEnergy > 0.0) 446 if (electronKineEnergy > 0.0) 450 { 447 { 451 G4ThreeVector electronDirection ( dirX_e 448 G4ThreeVector electronDirection ( dirX_el, dirY_el, dirZ_el); 452 electronDirection.rotateUz(photonDirecti 449 electronDirection.rotateUz(photonDirection); 453 G4DynamicParticle* electron = new G4Dyna 450 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 454 electronDirection, 451 electronDirection, 455 electronKineEnergy); 452 electronKineEnergy); 456 fvect->push_back(electron); 453 fvect->push_back(electron); 457 } 454 } 458 else 455 else 459 { 456 { 460 localEnergyDeposit += electronKineEnergy 457 localEnergyDeposit += electronKineEnergy; 461 electronKineEnergy = 0; 458 electronKineEnergy = 0; 462 } 459 } 463 460 464 //Generate the positron. Real particle in an 461 //Generate the positron. Real particle in any case, because it will annihilate. If below 465 //threshold, produce it at rest 462 //threshold, produce it at rest 466 if (positronKineEnergy < 0.0) 463 if (positronKineEnergy < 0.0) 467 { 464 { 468 localEnergyDeposit += positronKineEnergy 465 localEnergyDeposit += positronKineEnergy; 469 positronKineEnergy = 0; //produce it at 466 positronKineEnergy = 0; //produce it at rest 470 } 467 } 471 G4ThreeVector positronDirection(dirX_po,dirY 468 G4ThreeVector positronDirection(dirX_po,dirY_po,dirZ_po); 472 positronDirection.rotateUz(photonDirection); 469 positronDirection.rotateUz(photonDirection); 473 G4DynamicParticle* positron = new G4DynamicP 470 G4DynamicParticle* positron = new G4DynamicParticle(G4Positron::Positron(), 474 positronDirection, positronK 471 positronDirection, positronKineEnergy); 475 fvect->push_back(positron); 472 fvect->push_back(positron); 476 473 477 //Add rest of energy to the local energy dep 474 //Add rest of energy to the local energy deposit 478 fParticleChange->ProposeLocalEnergyDeposit(l 475 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit); 479 476 480 if (fVerboseLevel > 1) << 477 if (verboseLevel > 1) 481 { 478 { 482 G4cout << "----------------------------- 479 G4cout << "-----------------------------------------------------------" << G4endl; 483 G4cout << "Energy balance from G4Penelop 480 G4cout << "Energy balance from G4PenelopeGammaConversion" << G4endl; 484 G4cout << "Incoming photon energy: " << 481 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl; 485 G4cout << "----------------------------- 482 G4cout << "-----------------------------------------------------------" << G4endl; 486 if (electronKineEnergy) 483 if (electronKineEnergy) 487 G4cout << "Electron (explicitly produced) " << 484 G4cout << "Electron (explicitely produced) " << electronKineEnergy/keV << " keV" 488 << G4endl; 485 << G4endl; 489 if (positronKineEnergy) 486 if (positronKineEnergy) 490 G4cout << "Positron (not at rest) " << posit 487 G4cout << "Positron (not at rest) " << positronKineEnergy/keV << " keV" << G4endl; 491 G4cout << "Rest masses of e+/- " << 2.0* 488 G4cout << "Rest masses of e+/- " << 2.0*electron_mass_c2/keV << " keV" << G4endl; 492 if (localEnergyDeposit) 489 if (localEnergyDeposit) 493 G4cout << "Local energy deposit " << localEn 490 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl; 494 G4cout << "Total final state: " << (elec 491 G4cout << "Total final state: " << (electronKineEnergy+positronKineEnergy+ 495 localEnergyDeposit+2.0*electron_ma 492 localEnergyDeposit+2.0*electron_mass_c2)/keV << 496 " keV" << G4endl; 493 " keV" << G4endl; 497 G4cout << "----------------------------- 494 G4cout << "-----------------------------------------------------------" << G4endl; 498 } 495 } 499 if (fVerboseLevel > 0) << 496 if (verboseLevel > 0) 500 { 497 { 501 G4double energyDiff = std::fabs(electron 498 G4double energyDiff = std::fabs(electronKineEnergy+positronKineEnergy+ 502 localEnergyDeposit+2.0*electron_ 499 localEnergyDeposit+2.0*electron_mass_c2-photonEnergy); 503 if (energyDiff > 0.05*keV) 500 if (energyDiff > 0.05*keV) 504 G4cout << "Warning from G4PenelopeGammaConve 501 G4cout << "Warning from G4PenelopeGammaConversion: problem with energy conservation: " 505 << (electronKineEnergy+positronKineEn 502 << (electronKineEnergy+positronKineEnergy+ 506 localEnergyDeposit+2.0*electron_mass_c2 503 localEnergyDeposit+2.0*electron_mass_c2)/keV 507 << " keV (final) vs. " << photonEnerg 504 << " keV (final) vs. " << photonEnergy/keV << " keV (initial)" << G4endl; 508 } 505 } 509 } 506 } 510 507 511 //....oooOO0OOooo........oooOO0OOooo........oo 508 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 512 509 513 void G4PenelopeGammaConversionModel::ReadDataF 510 void G4PenelopeGammaConversionModel::ReadDataFile(const G4int Z) 514 { 511 { 515 if (!IsMaster()) 512 if (!IsMaster()) 516 //Should not be here! 513 //Should not be here! 517 G4Exception("G4PenelopeGammaConversionMode 514 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()", 518 "em0100",FatalException,"Worker thread in 515 "em0100",FatalException,"Worker thread in this method"); 519 516 520 if (fVerboseLevel > 2) << 517 if (verboseLevel > 2) 521 { 518 { 522 G4cout << "G4PenelopeGammaConversionMode 519 G4cout << "G4PenelopeGammaConversionModel::ReadDataFile()" << G4endl; 523 G4cout << "Going to read Gamma Conversio 520 G4cout << "Going to read Gamma Conversion data files for Z=" << Z << G4endl; 524 } 521 } 525 522 526 const char* path = G4FindDataDir("G4LEDATA << 523 char* path = getenv("G4LEDATA"); 527 if(!path) << 524 if (!path) 528 { 525 { 529 G4String excep = 526 G4String excep = 530 "G4PenelopeGammaConversionModel - G4LEDATA e 527 "G4PenelopeGammaConversionModel - G4LEDATA environment variable not set!"; 531 G4Exception("G4PenelopeGammaConversionMo 528 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()", 532 "em0006",FatalException,excep); 529 "em0006",FatalException,excep); 533 return; 530 return; 534 } 531 } 535 532 536 /* 533 /* 537 Read the cross section file 534 Read the cross section file 538 */ 535 */ 539 std::ostringstream ost; 536 std::ostringstream ost; 540 if (Z>9) 537 if (Z>9) 541 ost << path << "/penelope/pairproduction/p 538 ost << path << "/penelope/pairproduction/pdgpp" << Z << ".p08"; 542 else 539 else 543 ost << path << "/penelope/pairproduction/p 540 ost << path << "/penelope/pairproduction/pdgpp0" << Z << ".p08"; 544 std::ifstream file(ost.str().c_str()); 541 std::ifstream file(ost.str().c_str()); 545 if (!file.is_open()) 542 if (!file.is_open()) 546 { 543 { 547 G4String excep = "G4PenelopeGammaConvers 544 G4String excep = "G4PenelopeGammaConversionModel - data file " + 548 G4String(ost.str()) + " not found!"; 545 G4String(ost.str()) + " not found!"; 549 G4Exception("G4PenelopeGammaConversionMo 546 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()", 550 "em0003",FatalException,excep); 547 "em0003",FatalException,excep); 551 } 548 } 552 549 553 //I have to know in advance how many points 550 //I have to know in advance how many points are in the data list 554 //to initialize the G4PhysicsFreeVector() 551 //to initialize the G4PhysicsFreeVector() 555 std::size_t ndata=0; << 552 size_t ndata=0; 556 G4String line; 553 G4String line; 557 while( getline(file, line) ) 554 while( getline(file, line) ) 558 ndata++; 555 ndata++; 559 ndata -= 1; //remove one header line 556 ndata -= 1; //remove one header line >> 557 //G4cout << "Found: " << ndata << " lines" << G4endl; 560 558 561 file.clear(); 559 file.clear(); 562 file.close(); 560 file.close(); 563 file.open(ost.str().c_str()); 561 file.open(ost.str().c_str()); 564 G4int readZ =0; 562 G4int readZ =0; 565 file >> readZ; 563 file >> readZ; 566 564 567 if (fVerboseLevel > 3) << 565 if (verboseLevel > 3) 568 G4cout << "Element Z=" << Z << G4endl; 566 G4cout << "Element Z=" << Z << G4endl; 569 567 570 //check the right file is opened. 568 //check the right file is opened. 571 if (readZ != Z) 569 if (readZ != Z) 572 { 570 { 573 G4ExceptionDescription ed; 571 G4ExceptionDescription ed; 574 ed << "Corrupted data file for Z=" << Z 572 ed << "Corrupted data file for Z=" << Z << G4endl; 575 G4Exception("G4PenelopeGammaConversionMo 573 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()", 576 "em0005",FatalException,ed); 574 "em0005",FatalException,ed); 577 } 575 } 578 576 579 fLogAtomicCrossSection[Z] = new G4PhysicsFre << 577 G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(ndata); 580 G4double ene=0,xs=0; 578 G4double ene=0,xs=0; 581 for (std::size_t i=0;i<ndata;++i) << 579 for (size_t i=0;i<ndata;i++) 582 { 580 { 583 file >> ene >> xs; 581 file >> ene >> xs; 584 //dimensional quantities 582 //dimensional quantities 585 ene *= eV; 583 ene *= eV; 586 xs *= barn; 584 xs *= barn; 587 if (xs < 1e-40*cm2) //protection against 585 if (xs < 1e-40*cm2) //protection against log(0) 588 xs = 1e-40*cm2; 586 xs = 1e-40*cm2; 589 fLogAtomicCrossSection[Z]->PutValue(i,G4 << 587 theVec->PutValue(i,std::log(ene),std::log(xs)); 590 } 588 } 591 file.close(); 589 file.close(); 592 590 >> 591 if (!logAtomicCrossSection) >> 592 { >> 593 G4ExceptionDescription ed; >> 594 ed << "Problem with allocation of logAtomicCrossSection data table " << G4endl; >> 595 G4Exception("G4PenelopeGammaConversionModel::ReadDataFile()", >> 596 "em2020",FatalException,ed); >> 597 delete theVec; >> 598 return; >> 599 } >> 600 logAtomicCrossSection->insert(std::make_pair(Z,theVec)); >> 601 593 return; 602 return; >> 603 >> 604 } >> 605 >> 606 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 607 >> 608 void G4PenelopeGammaConversionModel::InitializeScreeningRadii() >> 609 { >> 610 G4double temp[99] = {1.2281e+02,7.3167e+01,6.9228e+01,6.7301e+01,6.4696e+01, >> 611 6.1228e+01,5.7524e+01,5.4033e+01,5.0787e+01,4.7851e+01,4.6373e+01, >> 612 4.5401e+01,4.4503e+01,4.3815e+01,4.3074e+01,4.2321e+01,4.1586e+01, >> 613 4.0953e+01,4.0524e+01,4.0256e+01,3.9756e+01,3.9144e+01,3.8462e+01, >> 614 3.7778e+01,3.7174e+01,3.6663e+01,3.5986e+01,3.5317e+01,3.4688e+01, >> 615 3.4197e+01,3.3786e+01,3.3422e+01,3.3068e+01,3.2740e+01,3.2438e+01, >> 616 3.2143e+01,3.1884e+01,3.1622e+01,3.1438e+01,3.1142e+01,3.0950e+01, >> 617 3.0758e+01,3.0561e+01,3.0285e+01,3.0097e+01,2.9832e+01,2.9581e+01, >> 618 2.9411e+01,2.9247e+01,2.9085e+01,2.8930e+01,2.8721e+01,2.8580e+01, >> 619 2.8442e+01,2.8312e+01,2.8139e+01,2.7973e+01,2.7819e+01,2.7675e+01, >> 620 2.7496e+01,2.7285e+01,2.7093e+01,2.6911e+01,2.6705e+01,2.6516e+01, >> 621 2.6304e+01,2.6108e+01,2.5929e+01,2.5730e+01,2.5577e+01,2.5403e+01, >> 622 2.5245e+01,2.5100e+01,2.4941e+01,2.4790e+01,2.4655e+01,2.4506e+01, >> 623 2.4391e+01,2.4262e+01,2.4145e+01,2.4039e+01,2.3922e+01,2.3813e+01, >> 624 2.3712e+01,2.3621e+01,2.3523e+01,2.3430e+01,2.3331e+01,2.3238e+01, >> 625 2.3139e+01,2.3048e+01,2.2967e+01,2.2833e+01,2.2694e+01,2.2624e+01, >> 626 2.2545e+01,2.2446e+01,2.2358e+01,2.2264e+01}; >> 627 >> 628 //copy temporary vector in class data member >> 629 for (G4int i=0;i<99;i++) >> 630 fAtomicScreeningRadius[i] = temp[i]; 594 } 631 } 595 632 596 //....oooOO0OOooo........oooOO0OOooo........oo 633 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 597 634 598 void G4PenelopeGammaConversionModel::Initializ 635 void G4PenelopeGammaConversionModel::InitializeScreeningFunctions(const G4Material* material) 599 { 636 { >> 637 /* >> 638 if (!IsMaster()) >> 639 //Should not be here! >> 640 G4Exception("G4PenelopeGammaConversionModel::InitializeScreeningFunctions()", >> 641 "em01001",FatalException,"Worker thread in this method"); >> 642 */ >> 643 600 // This is subroutine GPPa0 of Penelope 644 // This is subroutine GPPa0 of Penelope 601 // 645 // 602 // 1) calculate the effective Z for the purp 646 // 1) calculate the effective Z for the purpose 603 // 647 // 604 G4double zeff = 0; 648 G4double zeff = 0; 605 G4int intZ = 0; 649 G4int intZ = 0; 606 G4int nElements = (G4int)material->GetNumber << 650 G4int nElements = material->GetNumberOfElements(); 607 const G4ElementVector* elementVector = mater 651 const G4ElementVector* elementVector = material->GetElementVector(); 608 652 609 //avoid calculations if only one building el 653 //avoid calculations if only one building element! 610 if (nElements == 1) 654 if (nElements == 1) 611 { 655 { 612 zeff = (*elementVector)[0]->GetZ(); 656 zeff = (*elementVector)[0]->GetZ(); 613 intZ = (G4int) zeff; 657 intZ = (G4int) zeff; 614 } 658 } 615 else // many elements...let's do the calcula 659 else // many elements...let's do the calculation 616 { 660 { 617 const G4double* fractionVector = materia 661 const G4double* fractionVector = material->GetVecNbOfAtomsPerVolume(); 618 662 619 G4double atot = 0; 663 G4double atot = 0; 620 for (G4int i=0;i<nElements;i++) 664 for (G4int i=0;i<nElements;i++) 621 { 665 { 622 G4double Zelement = (*elementVector)[i]->G 666 G4double Zelement = (*elementVector)[i]->GetZ(); 623 G4double Aelement = (*elementVector)[i]->G 667 G4double Aelement = (*elementVector)[i]->GetAtomicMassAmu(); 624 atot += Aelement*fractionVector[i]; 668 atot += Aelement*fractionVector[i]; 625 zeff += Zelement*Aelement*fractionVector[i 669 zeff += Zelement*Aelement*fractionVector[i]; //average with the number of nuclei 626 } 670 } 627 atot /= material->GetTotNbOfAtomsPerVolu 671 atot /= material->GetTotNbOfAtomsPerVolume(); 628 zeff /= (material->GetTotNbOfAtomsPerVol 672 zeff /= (material->GetTotNbOfAtomsPerVolume()*atot); 629 673 630 intZ = (G4int) (zeff+0.25); 674 intZ = (G4int) (zeff+0.25); 631 if (intZ <= 0) 675 if (intZ <= 0) 632 intZ = 1; 676 intZ = 1; 633 if (intZ > fMaxZ) << 677 if (intZ > 99) 634 intZ = fMaxZ; << 678 intZ = 99; 635 } 679 } 636 680 637 if (fEffectiveCharge) 681 if (fEffectiveCharge) 638 fEffectiveCharge->insert(std::make_pair(ma 682 fEffectiveCharge->insert(std::make_pair(material,zeff)); 639 683 640 // 684 // 641 // 2) Calculate Coulomb Correction 685 // 2) Calculate Coulomb Correction 642 // 686 // 643 G4double alz = fine_structure_const*zeff; 687 G4double alz = fine_structure_const*zeff; 644 G4double alzSquared = alz*alz; 688 G4double alzSquared = alz*alz; 645 G4double fc = alzSquared*(0.202059-alzSquar 689 G4double fc = alzSquared*(0.202059-alzSquared* 646 (0.03693-alzSquared* 690 (0.03693-alzSquared* 647 (0.00835-alzSquared*(0.00201-alzSq 691 (0.00835-alzSquared*(0.00201-alzSquared* 648 (0.00049-alzSquared* 692 (0.00049-alzSquared* 649 (0.00012-alzSquared*0.00003))) 693 (0.00012-alzSquared*0.00003))))) 650 +1.0/(alzSquared+1.0)); 694 +1.0/(alzSquared+1.0)); 651 // 695 // 652 // 3) Screening functions and low-energy cor 696 // 3) Screening functions and low-energy corrections 653 // 697 // 654 G4double matRadius = 2.0/ fAtomicScreeningRa << 698 G4double matRadius = 2.0/ fAtomicScreeningRadius[intZ-1]; 655 if (fMaterialInvScreeningRadius) 699 if (fMaterialInvScreeningRadius) 656 fMaterialInvScreeningRadius->insert(std::m 700 fMaterialInvScreeningRadius->insert(std::make_pair(material,matRadius)); 657 701 658 std::pair<G4double,G4double> myPair(0,0); 702 std::pair<G4double,G4double> myPair(0,0); 659 G4double f0a = 4.0*G4Log(fAtomicScreeningRad << 703 G4double f0a = 4.0*std::log(fAtomicScreeningRadius[intZ-1]); 660 G4double f0b = f0a - 4.0*fc; 704 G4double f0b = f0a - 4.0*fc; 661 myPair.first = f0a; 705 myPair.first = f0a; 662 myPair.second = f0b; 706 myPair.second = f0b; 663 707 664 if (fScreeningFunction) 708 if (fScreeningFunction) 665 fScreeningFunction->insert(std::make_pair( 709 fScreeningFunction->insert(std::make_pair(material,myPair)); 666 710 667 if (fVerboseLevel > 2) << 711 if (verboseLevel > 2) 668 { 712 { 669 G4cout << "Average Z for material " << m 713 G4cout << "Average Z for material " << material->GetName() << " = " << 670 zeff << G4endl; 714 zeff << G4endl; 671 G4cout << "Effective radius for material 715 G4cout << "Effective radius for material " << material->GetName() << " = " << 672 fAtomicScreeningRadius[intZ] << " m_e*c/hbar << 716 fAtomicScreeningRadius[intZ-1] << " m_e*c/hbar --> BCB = " << 673 matRadius << G4endl; 717 matRadius << G4endl; 674 G4cout << "Screening parameters F0 for m 718 G4cout << "Screening parameters F0 for material " << material->GetName() << " = " << 675 f0a << "," << f0b << G4endl; 719 f0a << "," << f0b << G4endl; 676 } 720 } 677 return; 721 return; 678 } 722 } 679 723 680 //....oooOO0OOooo........oooOO0OOooo........oo 724 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 681 725 682 std::pair<G4double,G4double> 726 std::pair<G4double,G4double> 683 G4PenelopeGammaConversionModel::GetScreeningFu 727 G4PenelopeGammaConversionModel::GetScreeningFunctions(G4double B) 684 { 728 { 685 // This is subroutine SCHIFF of Penelope 729 // This is subroutine SCHIFF of Penelope 686 // 730 // 687 // Screening Functions F1(B) and F2(B) in th 731 // Screening Functions F1(B) and F2(B) in the Bethe-Heitler differential cross 688 // section for pair production 732 // section for pair production 689 // 733 // 690 std::pair<G4double,G4double> result(0.,0.); 734 std::pair<G4double,G4double> result(0.,0.); 691 G4double BSquared = B*B; 735 G4double BSquared = B*B; 692 G4double f1 = 2.0-2.0*G4Log(1.0+BSquared); << 736 G4double f1 = 2.0-2.0*std::log(1.0+BSquared); 693 G4double f2 = f1 - 6.66666666e-1; // (-2/3) 737 G4double f2 = f1 - 6.66666666e-1; // (-2/3) 694 if (B < 1.0e-10) 738 if (B < 1.0e-10) 695 f1 = f1-twopi*B; 739 f1 = f1-twopi*B; 696 else 740 else 697 { 741 { 698 G4double a0 = 4.0*B*std::atan(1./B); 742 G4double a0 = 4.0*B*std::atan(1./B); 699 f1 = f1 - a0; 743 f1 = f1 - a0; 700 f2 += 2.0*BSquared*(4.0-a0-3.0*G4Log((1. << 744 f2 += 2.0*BSquared*(4.0-a0-3.0*std::log((1.0+BSquared)/BSquared)); 701 } 745 } 702 G4double g1 = 0.5*(3.0*f1-f2); 746 G4double g1 = 0.5*(3.0*f1-f2); 703 G4double g2 = 0.25*(3.0*f1+f2); 747 G4double g2 = 0.25*(3.0*f1+f2); 704 748 705 result.first = g1; 749 result.first = g1; 706 result.second = g2; 750 result.second = g2; 707 751 708 return result; 752 return result; 709 } 753 } 710 754 711 //....oooOO0OOooo........oooOO0OOooo........oo 755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 712 756 713 void G4PenelopeGammaConversionModel::SetPartic 757 void G4PenelopeGammaConversionModel::SetParticle(const G4ParticleDefinition* p) 714 { 758 { 715 if(!fParticle) { 759 if(!fParticle) { 716 fParticle = p; 760 fParticle = p; 717 } 761 } 718 } 762 } 719 763