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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4eBremsstrahlungRelModel 32 // File name: G4eBremsstrahlungRelModel 33 // 33 // 34 // Author: Andreas Schaelicke 34 // Author: Andreas Schaelicke 35 // 35 // 36 // Creation date: 12.08.2008 36 // Creation date: 12.08.2008 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // 40 // 13.11.08 add SetLPMflag and SetLPMconsta 40 // 13.11.08 add SetLPMflag and SetLPMconstant methods 41 // 13.11.08 change default LPMconstant valu 41 // 13.11.08 change default LPMconstant value 42 // 13.10.10 add angular distributon interfa 42 // 13.10.10 add angular distributon interface (VI) 43 // 31.05.16 change LPMconstant such that it 43 // 31.05.16 change LPMconstant such that it gives suppression variable 's' 44 // that consistent to Migdal's one 44 // that consistent to Migdal's one; fix a small bug in 'logTS1' 45 // computation; better agreement w 45 // computation; better agreement with exp.(M.Novak) 46 // 15.07.18 improved LPM suppression functi 46 // 15.07.18 improved LPM suppression function approximation (no artificial 47 // steps), code cleanup and optimi 47 // steps), code cleanup and optimizations,more implementation and 48 // model related comments, consist 48 // model related comments, consistent variable naming (M.Novak) 49 // 49 // 50 // Main References: 50 // Main References: 51 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815; 51 // Y.-S.Tsai, Rev. Mod. Phys. 46 (1974) 815; Rev. Mod. Phys. 49 (1977) 421. 52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501. 52 // S.Klein, Rev. Mod. Phys. 71 (1999) 1501. 53 // T.Stanev et.al., Phys. Rev. D25 (1982) 129 53 // T.Stanev et.al., Phys. Rev. D25 (1982) 1291. 54 // M.L.Ter-Mikaelian, High-energy Electromagn 54 // M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media, 55 // Wiley, 1972. 55 // Wiley, 1972. 56 // 56 // 57 // ------------------------------------------- 57 // ------------------------------------------------------------------- 58 // 58 // 59 59 60 #include "G4eBremsstrahlungRelModel.hh" 60 #include "G4eBremsstrahlungRelModel.hh" 61 #include "G4PhysicalConstants.hh" 61 #include "G4PhysicalConstants.hh" 62 #include "G4SystemOfUnits.hh" 62 #include "G4SystemOfUnits.hh" 63 #include "G4Electron.hh" 63 #include "G4Electron.hh" 64 #include "G4Gamma.hh" 64 #include "G4Gamma.hh" 65 #include "Randomize.hh" 65 #include "Randomize.hh" 66 #include "G4Material.hh" 66 #include "G4Material.hh" 67 #include "G4Element.hh" 67 #include "G4Element.hh" 68 #include "G4ElementVector.hh" 68 #include "G4ElementVector.hh" 69 #include "G4ParticleChangeForLoss.hh" 69 #include "G4ParticleChangeForLoss.hh" 70 #include "G4ModifiedTsai.hh" 70 #include "G4ModifiedTsai.hh" 71 #include "G4Exp.hh" 71 #include "G4Exp.hh" 72 #include "G4Log.hh" 72 #include "G4Log.hh" 73 #include "G4Pow.hh" << 74 #include "G4EmParameters.hh" << 75 #include "G4AutoLock.hh" << 76 #include <thread> << 77 73 78 const G4int G4eBremsstrahlungRelModel::gMaxZet 74 const G4int G4eBremsstrahlungRelModel::gMaxZet = 120; 79 75 80 // constant DCS factor: 16\alpha r_0^2/3 76 // constant DCS factor: 16\alpha r_0^2/3 81 const G4double G4eBremsstrahlungRelModel::gBre 77 const G4double G4eBremsstrahlungRelModel::gBremFactor 82 = 16. * CLHEP::fine_structure_const * CLHEP: 78 = 16. * CLHEP::fine_structure_const * CLHEP::classic_electr_radius 83 * CLHEP::classic_electr_radius/3.; 79 * CLHEP::classic_electr_radius/3.; 84 80 85 // Migdal's constant: 4\pi r_0*electron_reduce 81 // Migdal's constant: 4\pi r_0*electron_reduced_compton_wavelength^2 86 const G4double G4eBremsstrahlungRelModel::gMig 82 const G4double G4eBremsstrahlungRelModel::gMigdalConstant 87 = 4. * CLHEP::pi * CLHEP::classic_electr_rad 83 = 4. * CLHEP::pi * CLHEP::classic_electr_radius 88 * CLHEP::electron_Compton_length * CLHEP:: 84 * CLHEP::electron_Compton_length * CLHEP::electron_Compton_length; 89 85 90 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c) 86 // LPM constant: \alpha(mc^2)^2/(4\pi*\hbar c) 91 const G4double G4eBremsstrahlungRelModel::gLPM 87 const G4double G4eBremsstrahlungRelModel::gLPMconstant 92 = CLHEP::fine_structure_const * CLHEP::elect 88 = CLHEP::fine_structure_const * CLHEP::electron_mass_c2 93 * CLHEP::electron_mass_c2 / (4. * CLHEP::p 89 * CLHEP::electron_mass_c2 / (4. * CLHEP::pi * CLHEP::hbarc); 94 90 95 // abscissas and weights of an 8 point Gauss-L 91 // abscissas and weights of an 8 point Gauss-Legendre quadrature 96 // for numerical integration on [0,1] 92 // for numerical integration on [0,1] 97 const G4double G4eBremsstrahlungRelModel::gXGL 93 const G4double G4eBremsstrahlungRelModel::gXGL[] = { 98 1.98550718e-02, 1.01666761e-01, 2.37233795e- 94 1.98550718e-02, 1.01666761e-01, 2.37233795e-01, 4.08282679e-01, 99 5.91717321e-01, 7.62766205e-01, 8.98333239e- 95 5.91717321e-01, 7.62766205e-01, 8.98333239e-01, 9.80144928e-01 100 }; 96 }; 101 const G4double G4eBremsstrahlungRelModel::gWGL 97 const G4double G4eBremsstrahlungRelModel::gWGL[] = { 102 5.06142681e-02, 1.11190517e-01, 1.56853323e- 98 5.06142681e-02, 1.11190517e-01, 1.56853323e-01, 1.81341892e-01, 103 1.81341892e-01, 1.56853323e-01, 1.11190517e- 99 1.81341892e-01, 1.56853323e-01, 1.11190517e-01, 5.06142681e-02 104 }; 100 }; 105 101 106 // elastic and inelatic radiation logarithms f 102 // elastic and inelatic radiation logarithms for light elements (where the 107 // Thomas-Fermi model doesn't work): computed 103 // Thomas-Fermi model doesn't work): computed by using Dirac-Fock model of atom. 108 const G4double G4eBremsstrahlungRelModel::gFel 104 const G4double G4eBremsstrahlungRelModel::gFelLowZet [] = { 109 0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 105 0.0, 5.3104, 4.7935, 4.7402, 4.7112, 4.6694, 4.6134, 4.5520 110 }; 106 }; 111 const G4double G4eBremsstrahlungRelModel::gFin 107 const G4double G4eBremsstrahlungRelModel::gFinelLowZet[] = { 112 0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 108 0.0, 5.9173, 5.6125, 5.5377, 5.4728, 5.4174, 5.3688, 5.3236 113 }; 109 }; 114 110 115 // LPM supression functions evaluated at initi 111 // LPM supression functions evaluated at initialisation time 116 std::shared_ptr<G4eBremsstrahlungRelModel::LPM << 112 G4eBremsstrahlungRelModel::LPMFuncs G4eBremsstrahlungRelModel::gLPMFuncs; 117 { << 118 // We have to use shared pointer for the LPM << 119 // by the G4eBremsstrahlungRelModel used in << 120 // model is owned (well deleted) by (at leas << 121 // a G4SeltzerBergerModel which is owned by << 122 // which owned by a G4ThreadLocalSingleton<G << 123 // which is a static global and thus deleted << 124 // is deleted. << 125 static auto _instance = std::make_shared<G4e << 126 return _instance; << 127 } << 128 113 129 // special data structure per element i.e. per 114 // special data structure per element i.e. per Z 130 std::shared_ptr<std::vector<G4eBremsstrahlungR << 115 std::vector<G4eBremsstrahlungRelModel::ElementData*> G4eBremsstrahlungRelModel::gElementData; 131 { << 132 // Same code comment as for gLPMFuncs. << 133 static auto _instance = std::make_shared<std << 134 return _instance; << 135 } << 136 << 137 static std::once_flag applyOnce; << 138 << 139 namespace << 140 { << 141 G4Mutex theBremRelMutex = G4MUTEX_INITIALIZE << 142 } << 143 116 144 G4eBremsstrahlungRelModel::G4eBremsstrahlungRe 117 G4eBremsstrahlungRelModel::G4eBremsstrahlungRelModel(const G4ParticleDefinition* p, 145 118 const G4String& nam) 146 : G4VEmModel(nam), fLPMFuncs(gLPMFuncs()), fEl << 119 : G4VEmModel(nam) 147 { 120 { 148 fGammaParticle = G4Gamma::Gamma(); 121 fGammaParticle = G4Gamma::Gamma(); 149 // 122 // 150 fLowestKinEnergy = 1.0*CLHEP::MeV; << 123 fLowestKinEnergy = 1.0*MeV; 151 SetLowEnergyLimit(fLowestKinEnergy); 124 SetLowEnergyLimit(fLowestKinEnergy); 152 // 125 // 153 fLPMEnergyThreshold = 1.e+39; 126 fLPMEnergyThreshold = 1.e+39; 154 fLPMEnergy = 0.; 127 fLPMEnergy = 0.; >> 128 >> 129 SetLPMFlag(true); >> 130 // 155 SetAngularDistribution(new G4ModifiedTsai()) 131 SetAngularDistribution(new G4ModifiedTsai()); 156 // 132 // 157 if (nullptr != p) { 133 if (nullptr != p) { 158 SetParticle(p); 134 SetParticle(p); 159 } 135 } 160 } 136 } 161 137 162 G4eBremsstrahlungRelModel::~G4eBremsstrahlungR 138 G4eBremsstrahlungRelModel::~G4eBremsstrahlungRelModel() 163 { 139 { 164 if (fIsInitializer) { << 140 if (IsMaster()) { 165 // clear ElementData container 141 // clear ElementData container 166 for (auto const & ptr : *fElementData) { d << 142 for (std::size_t iz = 0; iz < gElementData.size(); ++iz) { 167 fElementData->clear(); << 143 if (nullptr != gElementData[iz]) { >> 144 delete gElementData[iz]; >> 145 } >> 146 } >> 147 gElementData.clear(); 168 // clear LPMFunctions (if any) 148 // clear LPMFunctions (if any) 169 if (fLPMFuncs->fIsInitialized) { << 149 if (LPMFlag()) { 170 fLPMFuncs->fLPMFuncG.clear(); << 150 gLPMFuncs.fLPMFuncG.clear(); 171 fLPMFuncs->fLPMFuncPhi.clear(); << 151 gLPMFuncs.fLPMFuncPhi.clear(); 172 fLPMFuncs->fIsInitialized = false; << 152 gLPMFuncs.fIsInitialized = false; 173 } 153 } 174 } 154 } 175 } 155 } 176 156 177 void G4eBremsstrahlungRelModel::Initialise(con 157 void G4eBremsstrahlungRelModel::Initialise(const G4ParticleDefinition* p, 178 con 158 const G4DataVector& cuts) 179 { 159 { 180 // parameters in each thread << 160 if (nullptr != p) { 181 if (fPrimaryParticle != p) { << 182 SetParticle(p); 161 SetParticle(p); 183 } 162 } 184 fUseLPM = G4EmParameters::Instance()->LPM(); << 185 fCurrentIZ = 0; 163 fCurrentIZ = 0; 186 << 164 // init element data and precompute LPM functions (only if lpmflag is true) 187 // init static element data and precompute L << 188 std::call_once(applyOnce, [this]() { fIsInit << 189 << 190 // for all treads and derived classes << 191 if (fIsInitializer || fElementData->empty()) << 192 G4AutoLock l(&theBremRelMutex); << 193 if (fElementData->empty()) { << 194 fElementData->resize(gMaxZet+1, nullptr) << 195 } << 196 InitialiseElementData(); << 197 InitLPMFunctions(); << 198 l.unlock(); << 199 } << 200 << 201 // element selectors are initialized in the << 202 if (IsMaster()) { 165 if (IsMaster()) { 203 InitialiseElementSelectors(p, cuts); << 166 InitialiseElementData(); >> 167 if (LPMFlag()) { InitLPMFunctions(); } >> 168 if (LowEnergyLimit() < HighEnergyLimit()) { >> 169 InitialiseElementSelectors(p, cuts); >> 170 } 204 } 171 } 205 // initialisation in all threads << 206 if (nullptr == fParticleChange) { 172 if (nullptr == fParticleChange) { 207 fParticleChange = GetParticleChangeForLoss 173 fParticleChange = GetParticleChangeForLoss(); 208 } 174 } 209 if (GetTripletModel()) { 175 if (GetTripletModel()) { 210 GetTripletModel()->Initialise(p, cuts); 176 GetTripletModel()->Initialise(p, cuts); 211 fIsScatOffElectron = true; 177 fIsScatOffElectron = true; 212 } 178 } 213 } 179 } 214 180 215 void G4eBremsstrahlungRelModel::InitialiseLoca 181 void G4eBremsstrahlungRelModel::InitialiseLocal(const G4ParticleDefinition*, 216 182 G4VEmModel* masterModel) 217 { 183 { 218 SetElementSelectors(masterModel->GetElementS << 184 if (LowEnergyLimit() < HighEnergyLimit()) { >> 185 SetElementSelectors(masterModel->GetElementSelectors()); >> 186 } 219 } 187 } 220 188 221 void G4eBremsstrahlungRelModel::SetParticle(co 189 void G4eBremsstrahlungRelModel::SetParticle(const G4ParticleDefinition* p) 222 { 190 { 223 fPrimaryParticle = p; 191 fPrimaryParticle = p; 224 fPrimaryParticleMass = p->GetPDGMass(); 192 fPrimaryParticleMass = p->GetPDGMass(); 225 fIsElectron = (p==G4Electron::Elect 193 fIsElectron = (p==G4Electron::Electron()); 226 } 194 } 227 195 228 // Sets kinematical variables like E_kin, E_t 196 // Sets kinematical variables like E_kin, E_t and some material dependent 229 // variables like LPM energy and characteristi 197 // variables like LPM energy and characteristic photon energy k_p (more exactly 230 // k_p^2) for the Ter-Mikaelian suppression ef 198 // k_p^2) for the Ter-Mikaelian suppression effect. 231 void G4eBremsstrahlungRelModel::SetupForMateri 199 void G4eBremsstrahlungRelModel::SetupForMaterial(const G4ParticleDefinition*, 232 200 const G4Material* mat, 233 201 G4double kineticEnergy) 234 { 202 { 235 fDensityFactor = gMigdalConstant*mat->GetEle 203 fDensityFactor = gMigdalConstant*mat->GetElectronDensity(); 236 fLPMEnergy = gLPMconstant*mat->GetRadlen 204 fLPMEnergy = gLPMconstant*mat->GetRadlen(); 237 // threshold for LPM effect (i.e. below whic 205 // threshold for LPM effect (i.e. below which LPM hidden by density effect) 238 if (fUseLPM) { << 206 if (LPMFlag()) { 239 fLPMEnergyThreshold = std::sqrt(fDensityFa 207 fLPMEnergyThreshold = std::sqrt(fDensityFactor)*fLPMEnergy; 240 } else { 208 } else { 241 fLPMEnergyThreshold = 1.e+39; // i.e. do 209 fLPMEnergyThreshold = 1.e+39; // i.e. do not use LPM effect 242 } 210 } 243 // calculate threshold for density effect: k 211 // calculate threshold for density effect: k_p = sqrt(fDensityCorr) 244 fPrimaryKinEnergy = kineticEnergy; 212 fPrimaryKinEnergy = kineticEnergy; 245 fPrimaryTotalEnergy = kineticEnergy+fPrimary 213 fPrimaryTotalEnergy = kineticEnergy+fPrimaryParticleMass; 246 fDensityCorr = fDensityFactor*fPrimar 214 fDensityCorr = fDensityFactor*fPrimaryTotalEnergy*fPrimaryTotalEnergy; 247 // set activation flag for LPM effects in th 215 // set activation flag for LPM effects in the DCS 248 fIsLPMActive = (fPrimaryTotalEnergy>fLPMEner << 216 fIsLPMActive = (fPrimaryTotalEnergy>fLPMEnergyThreshold); 249 } 217 } 250 218 251 // minimum primary (e-/e+) energy at which dis 219 // minimum primary (e-/e+) energy at which discrete interaction is possible 252 G4double G4eBremsstrahlungRelModel::MinPrimary 220 G4double G4eBremsstrahlungRelModel::MinPrimaryEnergy(const G4Material*, 253 221 const G4ParticleDefinition*, 254 222 G4double cut) 255 { 223 { 256 return std::max(fLowestKinEnergy, cut); 224 return std::max(fLowestKinEnergy, cut); 257 } 225 } 258 226 259 // Computes the restricted dE/dx as the approp 227 // Computes the restricted dE/dx as the appropriate weight of the individual 260 // element contributions that are computed by 228 // element contributions that are computed by numerically integrating the DCS. 261 G4double 229 G4double 262 G4eBremsstrahlungRelModel::ComputeDEDXPerVolum 230 G4eBremsstrahlungRelModel::ComputeDEDXPerVolume(const G4Material* material, 263 231 const G4ParticleDefinition* p, 264 232 G4double kineticEnergy, 265 233 G4double cutEnergy) 266 { 234 { 267 G4double dedx = 0.0; 235 G4double dedx = 0.0; 268 if (nullptr == fPrimaryParticle) { 236 if (nullptr == fPrimaryParticle) { 269 SetParticle(p); 237 SetParticle(p); 270 } 238 } 271 if (kineticEnergy < LowEnergyLimit()) { 239 if (kineticEnergy < LowEnergyLimit()) { 272 return dedx; 240 return dedx; 273 } 241 } 274 // maximum value of the dE/dx integral (the 242 // maximum value of the dE/dx integral (the minimum is 0 of course) 275 G4double tmax = std::min(cutEnergy, kineticE 243 G4double tmax = std::min(cutEnergy, kineticEnergy); 276 if (tmax == 0.0) { 244 if (tmax == 0.0) { 277 return dedx; 245 return dedx; 278 } 246 } 279 // sets kinematical and material related var 247 // sets kinematical and material related variables 280 SetupForMaterial(fPrimaryParticle, material, 248 SetupForMaterial(fPrimaryParticle, material,kineticEnergy); 281 // get element compositions of the material 249 // get element compositions of the material 282 const G4ElementVector* theElemVector = mater 250 const G4ElementVector* theElemVector = material->GetElementVector(); 283 const G4double* theAtomNumDensVector = mater 251 const G4double* theAtomNumDensVector = material->GetAtomicNumDensityVector(); 284 const std::size_t numberOfElements = theElem 252 const std::size_t numberOfElements = theElemVector->size(); 285 // loop over the elements of the material an 253 // loop over the elements of the material and compute their contributions to 286 // the restricted dE/dx by numerical integra 254 // the restricted dE/dx by numerical integration of the dependent part of DCS 287 for (std::size_t ie = 0; ie < numberOfElemen 255 for (std::size_t ie = 0; ie < numberOfElements; ++ie) { 288 G4VEmModel::SetCurrentElement((*theElemVec 256 G4VEmModel::SetCurrentElement((*theElemVector)[ie]); 289 G4int zet = (*theElemVector)[ie]->GetZasIn 257 G4int zet = (*theElemVector)[ie]->GetZasInt(); 290 fCurrentIZ = std::min(zet, gMaxZet); 258 fCurrentIZ = std::min(zet, gMaxZet); 291 dedx += (zet*zet)*theAtomNumD 259 dedx += (zet*zet)*theAtomNumDensVector[ie]*ComputeBremLoss(tmax); 292 } 260 } 293 // apply the constant factor C/Z = 16\alpha 261 // apply the constant factor C/Z = 16\alpha r_0^2/3 294 dedx *= gBremFactor; 262 dedx *= gBremFactor; 295 return std::max(dedx,0.); 263 return std::max(dedx,0.); 296 } 264 } 297 265 298 // Computes the integral part of the restricte 266 // Computes the integral part of the restricted dE/dx contribution from a given 299 // element (Z) by numerically integrating the 267 // element (Z) by numerically integrating the k dependent part of the DCS between 300 // k_min=0 and k_max = tmax = min[gamma-cut, e 268 // k_min=0 and k_max = tmax = min[gamma-cut, electron-kinetic-eenrgy]. 301 // The numerical integration is done by dividi 269 // The numerical integration is done by dividing the integration range into 'n' 302 // subintervals and an 8 pint GL integral (on 270 // subintervals and an 8 pint GL integral (on [0,1]) is performed on each sub- 303 // inteval by tranforming k to alpha=k/E_t (E_ 271 // inteval by tranforming k to alpha=k/E_t (E_t is the total energy of the e-) 304 // and each sub-interavl is transformed to [0, 272 // and each sub-interavl is transformed to [0,1]. So the integrastion is done 305 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delt 273 // in xi(alpha) = xi(k) = [k/E_t-alpha_i]/delta where alpha_i=(i-1)*delta for 306 // the i = 1,2,..,n-th sub-interval so xi(k) i 274 // the i = 1,2,..,n-th sub-interval so xi(k) in [0,1] on each sub-intevals. 307 // This transformation from 'k' to 'xi(k)' res 275 // This transformation from 'k' to 'xi(k)' results in a multiplicative factor 308 // of E_t*delta at each step. 276 // of E_t*delta at each step. 309 // The restricted dE/dx = N int_{0}^{k_max} k* 277 // The restricted dE/dx = N int_{0}^{k_max} k*ds/dk dk. There are 2 DCS model 310 // one with LPM and one without LPM effects (s 278 // one with LPM and one without LPM effects (see them below). In both case not 311 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is co 279 // the ds/dk(Z,k) but ds/dk(Z,k)*[F*k/C] is computed since: 312 // (i) what we need here is ds/dk*k and not 280 // (i) what we need here is ds/dk*k and not k so this multiplication is done 313 // (ii) the Ter-Mikaelian suppression i.e. F 281 // (ii) the Ter-Mikaelian suppression i.e. F related factor is done here 314 // (iii) the constant factor C (includes Z^2 282 // (iii) the constant factor C (includes Z^2 as well)is accounted in the caller 315 G4double G4eBremsstrahlungRelModel::ComputeBre 283 G4double G4eBremsstrahlungRelModel::ComputeBremLoss(G4double tmax) 316 { 284 { 317 // number of intervals and integration step 285 // number of intervals and integration step 318 const G4double alphaMax = tmax/fPrimaryTotal 286 const G4double alphaMax = tmax/fPrimaryTotalEnergy; 319 const G4int nSub = (G4int)(20*alphaMa 287 const G4int nSub = (G4int)(20*alphaMax)+3; 320 const G4double delta = alphaMax/((G4doubl 288 const G4double delta = alphaMax/((G4double)nSub); 321 // set minimum value of the first sub-inteva 289 // set minimum value of the first sub-inteval 322 G4double alpha_i = 0.0; 290 G4double alpha_i = 0.0; 323 G4double dedxInteg = 0.0; 291 G4double dedxInteg = 0.0; 324 for (G4int l = 0; l < nSub; ++l) { 292 for (G4int l = 0; l < nSub; ++l) { 325 for (G4int igl = 0; igl < 8; ++igl) { 293 for (G4int igl = 0; igl < 8; ++igl) { 326 // compute the emitted photon energy k 294 // compute the emitted photon energy k 327 const G4double k = (alpha_i+gXGL[igl]* 295 const G4double k = (alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy; 328 // compute the DCS value at k (without t 296 // compute the DCS value at k (without the constant, the 1/k, 1/F factors) 329 const G4double dcs = fIsLPMActive 297 const G4double dcs = fIsLPMActive 330 ? ComputeRelDXSectio 298 ? ComputeRelDXSectionPerAtom(k) // DCS WITHOUT LPM 331 : ComputeDXSectionPe 299 : ComputeDXSectionPerAtom(k); // DCS WITH LPM 332 // account Ter-Mikaelian suppression: ti 300 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2 333 dedxInteg += gWGL[igl]*dcs/(1.0+fDensity 301 dedxInteg += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k)); 334 } 302 } 335 // update sub-interval minimum value 303 // update sub-interval minimum value 336 alpha_i += delta; 304 alpha_i += delta; 337 } 305 } 338 // apply corrections due to variable transfo 306 // apply corrections due to variable transformation i.e. E_t*delta 339 dedxInteg *= delta*fPrimaryTotalEnergy; 307 dedxInteg *= delta*fPrimaryTotalEnergy; 340 return std::max(dedxInteg,0.); 308 return std::max(dedxInteg,0.); 341 } 309 } 342 310 343 // Computes restrected atomic cross section by 311 // Computes restrected atomic cross section by numerically integrating the 344 // DCS between the proper kinematical limits a 312 // DCS between the proper kinematical limits accounting the gamma production cut 345 G4double G4eBremsstrahlungRelModel::ComputeCro 313 G4double G4eBremsstrahlungRelModel::ComputeCrossSectionPerAtom( 346 314 const G4ParticleDefinition* p, 347 315 G4double kineticEnergy, 348 316 G4double Z, 349 317 G4double, 350 318 G4double cut, 351 319 G4double maxEnergy) 352 { 320 { 353 G4double crossSection = 0.0; 321 G4double crossSection = 0.0; 354 if (nullptr == fPrimaryParticle) { 322 if (nullptr == fPrimaryParticle) { 355 SetParticle(p); 323 SetParticle(p); 356 } 324 } 357 if (kineticEnergy < LowEnergyLimit()) { 325 if (kineticEnergy < LowEnergyLimit()) { 358 return crossSection; 326 return crossSection; 359 } 327 } 360 // min/max kinetic energy limits of the DCS 328 // min/max kinetic energy limits of the DCS integration: 361 const G4double tmin = std::min(cut, kineticE 329 const G4double tmin = std::min(cut, kineticEnergy); 362 const G4double tmax = std::min(maxEnergy, ki 330 const G4double tmax = std::min(maxEnergy, kineticEnergy); 363 // zero restricted x-section if e- kinetic e 331 // zero restricted x-section if e- kinetic energy is below gamma cut 364 if (tmin >= tmax) { 332 if (tmin >= tmax) { 365 return crossSection; 333 return crossSection; 366 } 334 } 367 fCurrentIZ = std::min(G4lrint(Z), gMaxZet); 335 fCurrentIZ = std::min(G4lrint(Z), gMaxZet); 368 // integrate numerically (dependent part of) 336 // integrate numerically (dependent part of) the DCS between the kin. limits: 369 // a. integrate between tmin and kineticEner 337 // a. integrate between tmin and kineticEnergy of the e- 370 crossSection = ComputeXSectionPerAtom(tmin); 338 crossSection = ComputeXSectionPerAtom(tmin); 371 // allow partial integration: only if maxEne 339 // allow partial integration: only if maxEnergy < kineticEnergy 372 // b. integrate between tmax and kineticEner 340 // b. integrate between tmax and kineticEnergy (tmax=maxEnergy in this case) 373 // (so the result in this case is the integr 341 // (so the result in this case is the integral of DCS between tmin and 374 // maxEnergy) 342 // maxEnergy) 375 if (tmax < kineticEnergy) { 343 if (tmax < kineticEnergy) { 376 crossSection -= ComputeXSectionPerAtom(tma 344 crossSection -= ComputeXSectionPerAtom(tmax); 377 } 345 } 378 // multiply with the constant factors: 16\al 346 // multiply with the constant factors: 16\alpha r_0^2/3 Z^2 379 crossSection *= Z*Z*gBremFactor; 347 crossSection *= Z*Z*gBremFactor; 380 return std::max(crossSection, 0.); 348 return std::max(crossSection, 0.); 381 } 349 } 382 350 383 // Numerical integral of the (k dependent part 351 // Numerical integral of the (k dependent part of) DCS between k_min=tmin and 384 // k_max = E_k (where E_k is the kinetic energ 352 // k_max = E_k (where E_k is the kinetic energy of the e- and tmin is the 385 // minimum of energy of the emitted photon). 353 // minimum of energy of the emitted photon). The integration is done in the 386 // transformed alpha(k) = ln(k/E_t) variable ( 354 // transformed alpha(k) = ln(k/E_t) variable (with E_t being the total energy of 387 // the primary e-). The integration range is d 355 // the primary e-). The integration range is divided into n sub-intervals with 388 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n wid 356 // delta = [ln(k_min/E_t)-ln(k_max/E_t)]/n width each. An 8 point GL integral 389 // on [0,1] is applied on each sub-inteval so 357 // on [0,1] is applied on each sub-inteval so alpha is transformed to 390 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/del 358 // xi(alpha) = xi(k) = [ln(k/E_t)-alpha_i]/delta where alpha_i = ln(k_min/E_t) + 391 // (i-1)*delta for the i = 1,2,..,n-th sub-int 359 // (i-1)*delta for the i = 1,2,..,n-th sub-interval and xi(k) in [0,1] on each 392 // sub-intevals. From the transformed xi, k(xi 360 // sub-intevals. From the transformed xi, k(xi) = E_t exp[xi*delta+alpha_i]. 393 // Since the integration is done in variable x 361 // Since the integration is done in variable xi instead of k this 394 // transformation results in a multiplicative 362 // transformation results in a multiplicative factor of k*delta at each step. 395 // However, DCS differential in k is ~1/k so t 363 // However, DCS differential in k is ~1/k so the multiplicative factor is simple 396 // becomes delta and the 1/k factor is dropped 364 // becomes delta and the 1/k factor is dropped from the DCS computation. 397 // NOTE: 365 // NOTE: 398 // - LPM suppression is accounted above thre 366 // - LPM suppression is accounted above threshold e- energy (corresponidng 399 // flag is set in SetUpForMaterial() => 2 367 // flag is set in SetUpForMaterial() => 2 DCS with/without LPM 400 // - Ter-Mikaelian suppression is always acc 368 // - Ter-Mikaelian suppression is always accounted 401 G4double G4eBremsstrahlungRelModel::ComputeXSe 369 G4double G4eBremsstrahlungRelModel::ComputeXSectionPerAtom(G4double tmin) 402 { 370 { 403 G4double xSection = 0.0; 371 G4double xSection = 0.0; 404 const G4double alphaMin = G4Log(tmin/fPrimar 372 const G4double alphaMin = G4Log(tmin/fPrimaryTotalEnergy); 405 const G4double alphaMax = G4Log(fPrimaryKinE << 373 const G4double alphaMax = G4Log(fPrimaryKinEnergy/fPrimaryTotalEnergy); 406 const G4int nSub = std::max((G4int)(0.45*alp << 374 const G4int nSub = (G4int)(0.45*(alphaMax-alphaMin))+4; 407 const G4double delta = alphaMax/((G4double)n << 375 const G4double delta = (alphaMax-alphaMin)/((G4double)nSub); 408 // set minimum value of the first sub-inteva 376 // set minimum value of the first sub-inteval 409 G4double alpha_i = alphaMin; << 377 G4double alpha_i = alphaMin; 410 for (G4int l = 0; l < nSub; ++l) { 378 for (G4int l = 0; l < nSub; ++l) { 411 for (G4int igl = 0; igl < 8; ++igl) { 379 for (G4int igl = 0; igl < 8; ++igl) { 412 // compute the emitted photon energy k 380 // compute the emitted photon energy k 413 const G4double k = G4Exp(alpha_i+gXGL[ 381 const G4double k = G4Exp(alpha_i+gXGL[igl]*delta)*fPrimaryTotalEnergy; 414 // compute the DCS value at k (without t 382 // compute the DCS value at k (without the constant, the 1/k, 1/F factors) 415 const G4double dcs = fIsLPMActive 383 const G4double dcs = fIsLPMActive 416 ? ComputeRelDXSectio 384 ? ComputeRelDXSectionPerAtom(k) // DCS WITHOUT LPM 417 : ComputeDXSectionPe 385 : ComputeDXSectionPerAtom(k); // DCS WITH LPM 418 // account Ter-Mikaelian suppression: ti 386 // account Ter-Mikaelian suppression: times 1/F with F = 1+(k_p/k)^2 419 xSection += gWGL[igl]*dcs/(1.0+fDensityC 387 xSection += gWGL[igl]*dcs/(1.0+fDensityCorr/(k*k)); 420 } 388 } 421 // update sub-interval minimum value 389 // update sub-interval minimum value 422 alpha_i += delta; 390 alpha_i += delta; 423 } 391 } 424 // apply corrections due to variable transfo 392 // apply corrections due to variable transformation 425 xSection *= delta; 393 xSection *= delta; 426 // final check 394 // final check 427 return std::max(xSection, 0.); 395 return std::max(xSection, 0.); 428 } 396 } 429 397 430 // DCS WITH LPM EFFECT: complete screening apr 398 // DCS WITH LPM EFFECT: complete screening aprx. and includes LPM suppression 431 // ds/dk(Z,k) = C/[F*k]*{ Xi(s*F)*[y^2*G/4 +(1 399 // ds/dk(Z,k) = C/[F*k]*{ Xi(s*F)*[y^2*G/4 +(1-y+y^2/3)Phi]*[L_el-f_c+L_inel/Z] 432 // +(1-y)*[1+1/Z]/12} 400 // +(1-y)*[1+1/Z]/12} with C = 16\alpha r_0^2/3 Z^2 and 433 // Xi(s),G(s), Phi(s) are LPM suppression func 401 // Xi(s),G(s), Phi(s) are LPM suppression functions: 434 // 402 // 435 // LPM SUPPRESSION: The 's' is the suppression 403 // LPM SUPPRESSION: The 's' is the suppression variable and F = F(k,k_p) = 436 // 1+(k_p/k)^2 with k_p = hbar*w_p*E/(m*c^2) i 404 // 1+(k_p/k)^2 with k_p = hbar*w_p*E/(m*c^2) is a material (e- density) 437 // dependent constant. F accounts the Ter-Mika 405 // dependent constant. F accounts the Ter-Mikaelian suppression with a smooth 438 // transition in the emitted photon energy. Al 406 // transition in the emitted photon energy. Also, the LPM suppression functions 439 // goes to 0 when s goes to 0 and goes to 1 wh 407 // goes to 0 when s goes to 0 and goes to 1 when s is increasing (=1 at s=~2) 440 // So evaluating the LPM suppression functions 408 // So evaluating the LPM suppression functions at 'sF' instead of 's' ensures a 441 // smooth transition depending on the emitted 409 // smooth transition depending on the emitted photon energy 'k': LPM effect is 442 // smoothly turned off i.e. Xi(sF)=G(sF)=Phi(s 410 // smoothly turned off i.e. Xi(sF)=G(sF)=Phi(sF)=1 when k << k_p because F >> 1 443 // and sF ~ s when k >> k_p since F ~ 1 in tha 411 // and sF ~ s when k >> k_p since F ~ 1 in that case. 444 // HERE, ds/dk(Z,k)*[F*k/C] is computed since: 412 // HERE, ds/dk(Z,k)*[F*k/C] is computed since: 445 // (i) DCS ~ 1/k factor will disappear due 413 // (i) DCS ~ 1/k factor will disappear due to the variable transformation 446 // v(k)=ln(k/E_t) -> dk/dv=E_t*e^v=k -> 414 // v(k)=ln(k/E_t) -> dk/dv=E_t*e^v=k -> ds/dv= ds/dk*dk/dv=ds/dk*k so it 447 // would cnacell out the 1/k factor => 415 // would cnacell out the 1/k factor => 1/k don't included here 448 // (ii) the constant factor C and Z don't de 416 // (ii) the constant factor C and Z don't depend on 'k' => not included here 449 // (iii) the 1/F(k) factor is accounted in th 417 // (iii) the 1/F(k) factor is accounted in the callers: explicitly (cross sec- 450 // tion computation) or implicitly thro 418 // tion computation) or implicitly through further variable transformaton 451 // (in the final state sampling algorit 419 // (in the final state sampling algorithm) 452 // COMPLETE SCREENING: see more at the DCS wit 420 // COMPLETE SCREENING: see more at the DCS without LPM effect below. 453 G4double 421 G4double 454 G4eBremsstrahlungRelModel::ComputeRelDXSection 422 G4eBremsstrahlungRelModel::ComputeRelDXSectionPerAtom(G4double gammaEnergy) 455 { 423 { 456 G4double dxsec = 0.0; 424 G4double dxsec = 0.0; 457 if (gammaEnergy < 0.) { 425 if (gammaEnergy < 0.) { 458 return dxsec; 426 return dxsec; 459 } 427 } 460 const G4double y = gammaEnergy/fPrimaryT 428 const G4double y = gammaEnergy/fPrimaryTotalEnergy; 461 const G4double onemy = 1.-y; 429 const G4double onemy = 1.-y; 462 const G4double dum0 = 0.25*y*y; 430 const G4double dum0 = 0.25*y*y; 463 // evaluate LPM functions (combined with the 431 // evaluate LPM functions (combined with the Ter-Mikaelian effect) 464 G4double funcGS, funcPhiS, funcXiS; 432 G4double funcGS, funcPhiS, funcXiS; 465 ComputeLPMfunctions(funcXiS, funcGS, funcPhi 433 ComputeLPMfunctions(funcXiS, funcGS, funcPhiS, gammaEnergy); 466 const ElementData* elDat = (*fElementData)[f << 434 const ElementData* elDat = gElementData[fCurrentIZ]; 467 const G4double term1 = funcXiS*(dum0*fun 435 const G4double term1 = funcXiS*(dum0*funcGS+(onemy+2.0*dum0)*funcPhiS); 468 dxsec = term1*elDat->fZFactor1+onemy*elDat-> 436 dxsec = term1*elDat->fZFactor1+onemy*elDat->fZFactor2; 469 // 437 // 470 if (fIsScatOffElectron) { 438 if (fIsScatOffElectron) { 471 fSumTerm = dxsec; 439 fSumTerm = dxsec; 472 fNucTerm = term1*elDat->fZFactor11 + onemy 440 fNucTerm = term1*elDat->fZFactor11 + onemy/12.; 473 } 441 } 474 return std::max(dxsec,0.0); 442 return std::max(dxsec,0.0); 475 } 443 } 476 444 477 // DCS WITHOUT LPM EFFECT: DCS with sceening ( 445 // DCS WITHOUT LPM EFFECT: DCS with sceening (Z>5) and Coulomb cor. no LPM 478 // ds/dk(Z,k)=C/[F*k]*{(1-y+3*y^2/4)*[(0.25*ph 446 // ds/dk(Z,k)=C/[F*k]*{(1-y+3*y^2/4)*[(0.25*phi1(g)-ln(Z)/3-f_c)+(0.25*psi1(e) 479 // -2*ln(Z)/3)/Z]+ (1-y)*[(phi1(g)-phi2(g))+(p 447 // -2*ln(Z)/3)/Z]+ (1-y)*[(phi1(g)-phi2(g))+(psi1(e)-psi2(e))/Z]/8} 480 // where f_c(Z) is the Coulomb correction fact 448 // where f_c(Z) is the Coulomb correction factor and phi1(g),phi2(g) and psi1(e), 481 // psi2(e) are coherent and incoherent screeni 449 // psi2(e) are coherent and incoherent screening functions. In the Thomas-Fermi 482 // model of the atom, the screening functions 450 // model of the atom, the screening functions will have a form that do not 483 // depend on Z (not explicitly). These numeric 451 // depend on Z (not explicitly). These numerical screening functions can be 484 // approximated as Tsai Eqs. [3.38-3.41] with 452 // approximated as Tsai Eqs. [3.38-3.41] with the variables g=gamma and 485 // e=epsilon given by Tsai Eqs. [3.30 and 3.31 453 // e=epsilon given by Tsai Eqs. [3.30 and 3.31] (see more details at the method 486 // ComputeScreeningFunctions()). Note, that in 454 // ComputeScreeningFunctions()). Note, that in case of complete screening i.e. 487 // g = e = 0 => 0.25*phi1(0)-ln(Z)/3 = ln(184. 455 // g = e = 0 => 0.25*phi1(0)-ln(Z)/3 = ln(184.149/Z^(1/3)) = L_el and 488 // 0.25*psi1(0)-2*ln(Z)/3=ln(1193.923/Z^(2/3)) 456 // 0.25*psi1(0)-2*ln(Z)/3=ln(1193.923/Z^(2/3))=L_inel and phi1(0)-phi2(0) = 489 // psi1(0)-psi2(0) = 2/3 so the DCS in complet 457 // psi1(0)-psi2(0) = 2/3 so the DCS in complete screening => 490 // COMPLETE SCREENING: 458 // COMPLETE SCREENING: 491 // ds/dk(Z,k)=C/k*{(1-y+3*y^2/4)*[L_el-f_c+L_i 459 // ds/dk(Z,k)=C/k*{(1-y+3*y^2/4)*[L_el-f_c+L_inel/Z] + (1-y)*[1+1/Z]/12} that is 492 // used in case of DCS with LPM above (if all 460 // used in case of DCS with LPM above (if all the suprression functions are 493 // absent i.e. their value = 1). 461 // absent i.e. their value = 1). 494 // Since the Thomas-Fermi model of the atom is 462 // Since the Thomas-Fermi model of the atom is not accurate at low Z, the DCS in 495 // complete screening is used here at low Z(<5 463 // complete screening is used here at low Z(<5) with L_el(Z), L_inel(Z) values 496 // computed by using the Dirac-Fock model of t 464 // computed by using the Dirac-Fock model of the atom. 497 // NOTE: that the Ter-Mikaelian suppression is 465 // NOTE: that the Ter-Mikaelian suppression is accounted in the DCS through the 498 // 1/F factor but it is included in the caller 466 // 1/F factor but it is included in the caller and not considered here. 499 // HERE, ds/dk(Z,k)*[F*k/C] is computed exactl 467 // HERE, ds/dk(Z,k)*[F*k/C] is computed exactly like in the DCS with LPM case. 500 G4double 468 G4double 501 G4eBremsstrahlungRelModel::ComputeDXSectionPer 469 G4eBremsstrahlungRelModel::ComputeDXSectionPerAtom(G4double gammaEnergy) 502 { 470 { 503 G4double dxsec = 0.0; 471 G4double dxsec = 0.0; 504 if (gammaEnergy < 0.) { 472 if (gammaEnergy < 0.) { 505 return dxsec; 473 return dxsec; 506 } 474 } 507 const G4double y = gammaEnergy/fPrim 475 const G4double y = gammaEnergy/fPrimaryTotalEnergy; 508 const G4double onemy = 1.-y; 476 const G4double onemy = 1.-y; 509 const G4double dum0 = onemy+0.75*y*y; 477 const G4double dum0 = onemy+0.75*y*y; 510 const ElementData* elDat = (*fElementData)[f << 478 const ElementData* elDat = gElementData[fCurrentIZ]; 511 // use complete screening and L_el, L_inel f 479 // use complete screening and L_el, L_inel from Dirac-Fock model instead of TF 512 if (fCurrentIZ < 5 || fIsUseCompleteScreenin 480 if (fCurrentIZ < 5 || fIsUseCompleteScreening) { 513 dxsec = dum0*elDat->fZFactor1; 481 dxsec = dum0*elDat->fZFactor1; 514 dxsec += onemy*elDat->fZFactor2; 482 dxsec += onemy*elDat->fZFactor2; 515 if (fIsScatOffElectron) { 483 if (fIsScatOffElectron) { 516 fSumTerm = dxsec; 484 fSumTerm = dxsec; 517 fNucTerm = dum0*elDat->fZFactor11+onemy/ 485 fNucTerm = dum0*elDat->fZFactor11+onemy/12.; 518 } 486 } 519 } else { 487 } else { 520 // use Tsai's analytical approx. (Tsai Eqs 488 // use Tsai's analytical approx. (Tsai Eqs. [3.38-3.41]) to the 'universal' 521 // numerical screening functions computed 489 // numerical screening functions computed by using the TF model of the atom 522 const G4double invZ = 1./(G4double)fCur 490 const G4double invZ = 1./(G4double)fCurrentIZ; 523 const G4double Fz = elDat->fFz; 491 const G4double Fz = elDat->fFz; 524 const G4double logZ = elDat->fLogZ; 492 const G4double logZ = elDat->fLogZ; 525 const G4double dum1 = y/(fPrimaryTotalE 493 const G4double dum1 = y/(fPrimaryTotalEnergy-gammaEnergy); 526 const G4double gamma = dum1*elDat->fGamm 494 const G4double gamma = dum1*elDat->fGammaFactor; 527 const G4double epsilon = dum1*elDat->fEpsi 495 const G4double epsilon = dum1*elDat->fEpsilonFactor; 528 // evaluate the screening functions 496 // evaluate the screening functions 529 G4double phi1, phi1m2, psi1, psi1m2; 497 G4double phi1, phi1m2, psi1, psi1m2; 530 ComputeScreeningFunctions(phi1, phi1m2, ps 498 ComputeScreeningFunctions(phi1, phi1m2, psi1, psi1m2, gamma, epsilon); 531 dxsec = dum0*((0.25*phi1-Fz) + (0.25*psi1 499 dxsec = dum0*((0.25*phi1-Fz) + (0.25*psi1-2.*logZ/3.)*invZ); 532 dxsec += 0.125*onemy*(phi1m2 + psi1m2*invZ 500 dxsec += 0.125*onemy*(phi1m2 + psi1m2*invZ); 533 if (fIsScatOffElectron) { 501 if (fIsScatOffElectron) { 534 fSumTerm = dxsec; 502 fSumTerm = dxsec; 535 fNucTerm = dum0*(0.25*phi1-Fz) + 0.125*o 503 fNucTerm = dum0*(0.25*phi1-Fz) + 0.125*onemy*phi1m2; 536 } 504 } 537 } 505 } 538 return std::max(dxsec,0.0); 506 return std::max(dxsec,0.0); 539 } 507 } 540 508 541 // Coherent and incoherent screening function 509 // Coherent and incoherent screening function approximations (see Tsai 542 // Eqs.[3.38-3.41]). Tsai's analytical approxi 510 // Eqs.[3.38-3.41]). Tsai's analytical approximations to the numerical screening 543 // functions computed by using the Thomas-Ferm 511 // functions computed by using the Thomas-Fermi model of atom (Moliere's appro- 544 // ximation to the numerical TF screening func 512 // ximation to the numerical TF screening function). In the TF-model, these 545 // screening functions can be expressed in a ' 513 // screening functions can be expressed in a 'universal' i.e. Z (directly) inde- 546 // pendent variable (see Tsai Eqs. Eqs. [3.30 514 // pendent variable (see Tsai Eqs. Eqs. [3.30 and 3.31]). 547 void G4eBremsstrahlungRelModel::ComputeScreeni 515 void G4eBremsstrahlungRelModel::ComputeScreeningFunctions(G4double& phi1, 548 516 G4double& phi1m2, 549 517 G4double& psi1, 550 518 G4double& psi1m2, 551 519 const G4double gam, 552 520 const G4double eps) 553 { 521 { 554 const G4double gam2 = gam*gam; 522 const G4double gam2 = gam*gam; 555 phi1 = 16.863-2.0*G4Log(1.0+0.311877*gam2) 523 phi1 = 16.863-2.0*G4Log(1.0+0.311877*gam2)+2.4*G4Exp(-0.9*gam) 556 +1.6*G4Exp(-1.5*gam); 524 +1.6*G4Exp(-1.5*gam); 557 phi1m2 = 2.0/(3.0+19.5*gam+18.0*gam2); // 525 phi1m2 = 2.0/(3.0+19.5*gam+18.0*gam2); // phi1-phi2 558 const G4double eps2 = eps*eps; 526 const G4double eps2 = eps*eps; 559 psi1 = 24.34-2.0*G4Log(1.0+13.111641*eps2) 527 psi1 = 24.34-2.0*G4Log(1.0+13.111641*eps2)+2.8*G4Exp(-8.0*eps) 560 +1.2*G4Exp(-29.2*eps); 528 +1.2*G4Exp(-29.2*eps); 561 psi1m2 = 2.0/(3.0+120.0*eps+1200.0*eps2); // 529 psi1m2 = 2.0/(3.0+120.0*eps+1200.0*eps2); //psi1-psi2 562 } 530 } 563 531 564 void 532 void 565 G4eBremsstrahlungRelModel::SampleSecondaries(s 533 G4eBremsstrahlungRelModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, 566 c 534 const G4MaterialCutsCouple* couple, 567 c 535 const G4DynamicParticle* dp, 568 G 536 G4double cutEnergy, 569 G 537 G4double maxEnergy) 570 { 538 { 571 const G4double kineticEnergy = dp->GetKin 539 const G4double kineticEnergy = dp->GetKineticEnergy(); >> 540 // const G4double logKineticEnergy = dp->GetLogKineticEnergy(); 572 if (kineticEnergy < LowEnergyLimit()) { 541 if (kineticEnergy < LowEnergyLimit()) { 573 return; 542 return; 574 } 543 } 575 // min, max kinetic energy limits 544 // min, max kinetic energy limits 576 const G4double tmin = std::min(cutEnergy, ki 545 const G4double tmin = std::min(cutEnergy, kineticEnergy); 577 const G4double tmax = std::min(maxEnergy, ki 546 const G4double tmax = std::min(maxEnergy, kineticEnergy); 578 if (tmin >= tmax) { 547 if (tmin >= tmax) { 579 return; 548 return; 580 } 549 } 581 // 550 // 582 SetupForMaterial(fPrimaryParticle, couple->G 551 SetupForMaterial(fPrimaryParticle, couple->GetMaterial(), kineticEnergy); 583 const G4Element* elm = SelectTargetAtom(coup 552 const G4Element* elm = SelectTargetAtom(couple,fPrimaryParticle,kineticEnergy, 584 dp-> 553 dp->GetLogKineticEnergy(),tmin,tmax); 585 // 554 // 586 fCurrentIZ = elm->GetZasInt(); 555 fCurrentIZ = elm->GetZasInt(); 587 const ElementData* elDat = (*fElementData)[f << 556 const ElementData* elDat = gElementData[fCurrentIZ]; 588 const G4double funcMax = elDat->fZFactor1+el 557 const G4double funcMax = elDat->fZFactor1+elDat->fZFactor2; 589 // get the random engine 558 // get the random engine 590 G4double rndm[2]; 559 G4double rndm[2]; 591 CLHEP::HepRandomEngine* rndmEngine = G4Rando 560 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine(); 592 // min max of the transformed variable: x(k) 561 // min max of the transformed variable: x(k) = ln(k^2+k_p^2) that is in [ln(k_c^2+k_p^2), ln(E_k^2+k_p^2)] 593 const G4double xmin = G4Log(tmin*tmin+fDen 562 const G4double xmin = G4Log(tmin*tmin+fDensityCorr); 594 const G4double xrange = G4Log(tmax*tmax+fDen 563 const G4double xrange = G4Log(tmax*tmax+fDensityCorr)-xmin; 595 G4double gammaEnergy, funcVal; 564 G4double gammaEnergy, funcVal; 596 do { 565 do { 597 rndmEngine->flatArray(2, rndm); 566 rndmEngine->flatArray(2, rndm); 598 gammaEnergy = std::sqrt(std::max(G4Exp(xmi 567 gammaEnergy = std::sqrt(std::max(G4Exp(xmin+rndm[0]*xrange)-fDensityCorr, 0.0)); 599 funcVal = fIsLPMActive 568 funcVal = fIsLPMActive 600 ? ComputeRelDXSectionPerAtom( 569 ? ComputeRelDXSectionPerAtom(gammaEnergy) 601 : ComputeDXSectionPerAtom(gam 570 : ComputeDXSectionPerAtom(gammaEnergy); 602 // cross-check of proper function maximum 571 // cross-check of proper function maximum in the rejection 603 // if (funcVal > funcMax) { 572 // if (funcVal > funcMax) { 604 // G4cout << "### G4eBremsstrahlungRelMod 573 // G4cout << "### G4eBremsstrahlungRelModel Warning: Majoranta exceeded! " 605 // << funcVal << " > " << funcMax 574 // << funcVal << " > " << funcMax 606 // << " Egamma(MeV)= " << gammaEnergy 575 // << " Egamma(MeV)= " << gammaEnergy 607 // << " Ee(MeV)= " << kineticEnergy 576 // << " Ee(MeV)= " << kineticEnergy 608 // << " " << GetName() 577 // << " " << GetName() 609 // << G4endl; 578 // << G4endl; 610 // } 579 // } 611 // Loop checking, 03-Aug-2015, Vladimir Iv 580 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 612 } while (funcVal < funcMax*rndm[1]); 581 } while (funcVal < funcMax*rndm[1]); 613 // 582 // 614 // scattering off nucleus or off e- by tripl 583 // scattering off nucleus or off e- by triplet model 615 if (fIsScatOffElectron && rndmEngine->flat() 584 if (fIsScatOffElectron && rndmEngine->flat()*fSumTerm>fNucTerm) { 616 GetTripletModel()->SampleSecondaries(vdp, 585 GetTripletModel()->SampleSecondaries(vdp, couple, dp, cutEnergy, maxEnergy); 617 return; 586 return; 618 } 587 } 619 // 588 // 620 // angles of the emitted gamma. ( Z - axis a 589 // angles of the emitted gamma. ( Z - axis along the parent particle) 621 // use general interface 590 // use general interface 622 G4ThreeVector gamDir = 591 G4ThreeVector gamDir = 623 GetAngularDistribution()->SampleDirection( 592 GetAngularDistribution()->SampleDirection(dp,fPrimaryTotalEnergy-gammaEnergy, 624 593 fCurrentIZ, couple->GetMaterial()); 625 // create G4DynamicParticle object for the G 594 // create G4DynamicParticle object for the Gamma 626 auto gamma = new G4DynamicParticle(fGammaPar 595 auto gamma = new G4DynamicParticle(fGammaParticle, gamDir, gammaEnergy); 627 vdp->push_back(gamma); 596 vdp->push_back(gamma); 628 // compute post-interaction kinematics of pr 597 // compute post-interaction kinematics of primary e-/e+ based on 629 // energy-momentum conservation 598 // energy-momentum conservation 630 const G4double totMomentum = std::sqrt(kinet 599 const G4double totMomentum = std::sqrt(kineticEnergy*( 631 fPrimaryTotalEn 600 fPrimaryTotalEnergy + CLHEP::electron_mass_c2)); 632 G4ThreeVector dir = 601 G4ThreeVector dir = 633 (totMomentum*dp->GetMomentumDirec 602 (totMomentum*dp->GetMomentumDirection()-gammaEnergy*gamDir).unit(); 634 const G4double finalE = kineticEnergy-gamm 603 const G4double finalE = kineticEnergy-gammaEnergy; 635 // if secondary gamma energy is higher than 604 // if secondary gamma energy is higher than threshold(very high by default) 636 // then stop tracking the primary particle a 605 // then stop tracking the primary particle and create new secondary e-/e+ 637 // instead of the primary one 606 // instead of the primary one 638 if (gammaEnergy > SecondaryThreshold()) { 607 if (gammaEnergy > SecondaryThreshold()) { 639 fParticleChange->ProposeTrackStatus(fStopA 608 fParticleChange->ProposeTrackStatus(fStopAndKill); 640 fParticleChange->SetProposedKineticEnergy( 609 fParticleChange->SetProposedKineticEnergy(0.0); 641 auto el = new G4DynamicParticle( 610 auto el = new G4DynamicParticle( 642 const_cast<G4ParticleDefinition* 611 const_cast<G4ParticleDefinition*>(fPrimaryParticle), dir, finalE); 643 vdp->push_back(el); 612 vdp->push_back(el); 644 } else { // continue tracking the primary e- 613 } else { // continue tracking the primary e-/e+ otherwise 645 fParticleChange->SetProposedMomentumDirect 614 fParticleChange->SetProposedMomentumDirection(dir); 646 fParticleChange->SetProposedKineticEnergy( 615 fParticleChange->SetProposedKineticEnergy(finalE); 647 } 616 } 648 } 617 } 649 618 650 void G4eBremsstrahlungRelModel::InitialiseElem 619 void G4eBremsstrahlungRelModel::InitialiseElementData() 651 { 620 { >> 621 const G4int size = (G4int)gElementData.size(); >> 622 if (size < gMaxZet+1) { >> 623 gElementData.resize(gMaxZet+1, nullptr); >> 624 } 652 // create for all elements that are in the d 625 // create for all elements that are in the detector 653 auto elemTable = G4Element::GetElementTable( << 626 const G4ElementTable* elemTable = G4Element::GetElementTable(); 654 for (auto const & elem : *elemTable) { << 627 std::size_t numElems = (*elemTable).size(); 655 const G4double zet = elem->GetZ(); << 628 for (std::size_t ielem=0; ielem<numElems; ++ielem) { 656 const G4int izet = std::min(elem->GetZasIn << 629 const G4Element* elem = (*elemTable)[ielem]; 657 if (nullptr == (*fElementData)[izet]) { << 630 const G4double zet = elem->GetZ(); >> 631 const G4int izet = std::min(G4lrint(zet),gMaxZet); >> 632 if (!gElementData[izet]) { 658 auto elemData = new ElementData(); 633 auto elemData = new ElementData(); 659 const G4double fc = elem->GetfCoulomb(); 634 const G4double fc = elem->GetfCoulomb(); 660 G4double Fel = 1.; 635 G4double Fel = 1.; 661 G4double Finel = 1.; 636 G4double Finel = 1.; 662 elemData->fLogZ = G4Log(zet); 637 elemData->fLogZ = G4Log(zet); 663 elemData->fFz = elemData->fLogZ/3.+f 638 elemData->fFz = elemData->fLogZ/3.+fc; 664 if (izet < 5) { 639 if (izet < 5) { 665 Fel = gFelLowZet[izet]; 640 Fel = gFelLowZet[izet]; 666 Finel = gFinelLowZet[izet]; 641 Finel = gFinelLowZet[izet]; 667 } else { 642 } else { 668 Fel = G4Log(184.15) - elemData->f 643 Fel = G4Log(184.15) - elemData->fLogZ/3.; 669 Finel = G4Log(1194) - 2.*elemData->f 644 Finel = G4Log(1194) - 2.*elemData->fLogZ/3.; 670 } 645 } 671 const G4double z13 = G4Pow::GetInstance( << 646 const G4double z23 = std::pow(zet,2./3.); 672 const G4double z23 = z13*z13; << 647 const G4double z13 = std::pow(zet,1./3.); 673 elemData->fZFactor1 = (Fel-fc)+Fine 648 elemData->fZFactor1 = (Fel-fc)+Finel/zet; 674 elemData->fZFactor11 = (Fel-fc); // 649 elemData->fZFactor11 = (Fel-fc); // used only for the triplet 675 elemData->fZFactor2 = (1.+1./zet)/1 650 elemData->fZFactor2 = (1.+1./zet)/12.; 676 elemData->fVarS1 = z23/(184.15*1 651 elemData->fVarS1 = z23/(184.15*184.15); 677 elemData->fILVarS1Cond = 1./(G4Log(std 652 elemData->fILVarS1Cond = 1./(G4Log(std::sqrt(2.0)*elemData->fVarS1)); 678 elemData->fILVarS1 = 1./G4Log(elem 653 elemData->fILVarS1 = 1./G4Log(elemData->fVarS1); 679 elemData->fGammaFactor = 100.0*electro 654 elemData->fGammaFactor = 100.0*electron_mass_c2/z13; 680 elemData->fEpsilonFactor = 100.0*electro 655 elemData->fEpsilonFactor = 100.0*electron_mass_c2/z23; 681 (*fElementData)[izet] = elemData; << 656 gElementData[izet] = elemData; 682 } 657 } 683 } 658 } 684 } 659 } 685 660 686 void G4eBremsstrahlungRelModel::ComputeLPMfunc 661 void G4eBremsstrahlungRelModel::ComputeLPMfunctions(G4double& funcXiS, 687 662 G4double& funcGS, 688 663 G4double& funcPhiS, 689 664 const G4double egamma) 690 { 665 { 691 static const G4double sqrt2 = std::sqrt(2.); 666 static const G4double sqrt2 = std::sqrt(2.); 692 const G4double redegamma = egamma/fPrimar 667 const G4double redegamma = egamma/fPrimaryTotalEnergy; 693 const G4double varSprime = std::sqrt(0.12 668 const G4double varSprime = std::sqrt(0.125*redegamma*fLPMEnergy/ 694 ((1.0-redegamm 669 ((1.0-redegamma)*fPrimaryTotalEnergy)); 695 const ElementData* elDat = (*fElementData << 670 const ElementData* elDat = gElementData[fCurrentIZ]; 696 const G4double varS1 = elDat->fVarS1; 671 const G4double varS1 = elDat->fVarS1; 697 const G4double condition = sqrt2*varS1; 672 const G4double condition = sqrt2*varS1; 698 G4double funcXiSprime = 2.0; 673 G4double funcXiSprime = 2.0; 699 if (varSprime > 1.0) { 674 if (varSprime > 1.0) { 700 funcXiSprime = 1.0; 675 funcXiSprime = 1.0; 701 } else if (varSprime > condition) { 676 } else if (varSprime > condition) { 702 const G4double ilVarS1Cond = elDat->fILVar 677 const G4double ilVarS1Cond = elDat->fILVarS1Cond; 703 const G4double funcHSprime = G4Log(varSpri 678 const G4double funcHSprime = G4Log(varSprime)*ilVarS1Cond; 704 funcXiSprime = 1.0 + funcHSprime - 0.08*(1 679 funcXiSprime = 1.0 + funcHSprime - 0.08*(1.0-funcHSprime)*funcHSprime 705 *(2.0-fu 680 *(2.0-funcHSprime)*ilVarS1Cond; 706 } 681 } 707 const G4double varS = varSprime/std::sqrt 682 const G4double varS = varSprime/std::sqrt(funcXiSprime); 708 // - include dielectric suppression effect i 683 // - include dielectric suppression effect into s according to Migdal 709 const G4double varShat = varS*(1.0+fDensityC 684 const G4double varShat = varS*(1.0+fDensityCorr/(egamma*egamma)); 710 funcXiS = 2.0; 685 funcXiS = 2.0; 711 if (varShat > 1.0) { 686 if (varShat > 1.0) { 712 funcXiS = 1.0; 687 funcXiS = 1.0; 713 } else if (varShat > varS1) { 688 } else if (varShat > varS1) { 714 funcXiS = 1.0+G4Log(varShat)*elDat->fILVar 689 funcXiS = 1.0+G4Log(varShat)*elDat->fILVarS1; 715 } 690 } 716 GetLPMFunctions(funcGS, funcPhiS, varShat); 691 GetLPMFunctions(funcGS, funcPhiS, varShat); 717 //ComputeLPMGsPhis(funcGS, funcPhiS, varShat 692 //ComputeLPMGsPhis(funcGS, funcPhiS, varShat); 718 // 693 // 719 //MAKE SURE SUPPRESSION IS SMALLER THAN 1: d 694 //MAKE SURE SUPPRESSION IS SMALLER THAN 1: due to Migdal's approximation on xi 720 if (funcXiS*funcPhiS > 1. || varShat > 0.57) 695 if (funcXiS*funcPhiS > 1. || varShat > 0.57) { 721 funcXiS=1./funcPhiS; 696 funcXiS=1./funcPhiS; 722 } 697 } 723 } 698 } 724 699 725 void G4eBremsstrahlungRelModel::ComputeLPMGsPh 700 void G4eBremsstrahlungRelModel::ComputeLPMGsPhis(G4double& funcGS, 726 701 G4double& funcPhiS, 727 702 const G4double varShat) 728 { 703 { 729 if (varShat < 0.01) { 704 if (varShat < 0.01) { 730 funcPhiS = 6.0*varShat*(1.0-CLHEP::pi*varS 705 funcPhiS = 6.0*varShat*(1.0-CLHEP::pi*varShat); 731 funcGS = 12.0*varShat-2.0*funcPhiS; 706 funcGS = 12.0*varShat-2.0*funcPhiS; 732 } else { 707 } else { 733 const G4double varShat2 = varShat*varShat; 708 const G4double varShat2 = varShat*varShat; 734 const G4double varShat3 = varShat*varShat2 709 const G4double varShat3 = varShat*varShat2; 735 const G4double varShat4 = varShat2*varShat 710 const G4double varShat4 = varShat2*varShat2; 736 // use Stanev approximation: for \psi(s) a 711 // use Stanev approximation: for \psi(s) and compute G(s) 737 if (varShat < 0.415827) { 712 if (varShat < 0.415827) { 738 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+v 713 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi)) 739 + varShat3/(0.623+0.796*varSha 714 + varShat3/(0.623+0.796*varShat+0.658*varShat2)); 740 // 1-\exp \left\{-4s-\frac{8s^2}{1+3.936 715 // 1-\exp \left\{-4s-\frac{8s^2}{1+3.936s+4.97s^2-0.05s^3+7.5s^4} \right\} 741 const G4double funcPsiS = 1.0 - G4Exp(-4 716 const G4double funcPsiS = 1.0 - G4Exp(-4.0*varShat 742 - 8.0*varShat2/ 717 - 8.0*varShat2/(1.0+3.936*varShat+4.97*varShat2 743 - 0.05*varShat3 718 - 0.05*varShat3 + 7.5*varShat4)); 744 // G(s) = 3 \psi(s) - 2 \phi(s) 719 // G(s) = 3 \psi(s) - 2 \phi(s) 745 funcGS = 3.0*funcPsiS - 2.0*funcPhiS; 720 funcGS = 3.0*funcPsiS - 2.0*funcPhiS; 746 } else if (varShat<1.55) { 721 } else if (varShat<1.55) { 747 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+v 722 funcPhiS = 1.0-G4Exp(-6.0*varShat*(1.0+varShat*(3.0-CLHEP::pi)) 748 + varShat3/(0.623+0.796*varSha 723 + varShat3/(0.623+0.796*varShat+0.658*varShat2)); 749 const G4double dum0 = -0.160723 724 const G4double dum0 = -0.160723 + 3.755030*varShat 750 -1.798138*varShat 725 -1.798138*varShat2 + 0.672827*varShat3 751 -0.120772*varShat 726 -0.120772*varShat4; 752 funcGS = std::tanh(dum0); 727 funcGS = std::tanh(dum0); 753 } else { 728 } else { 754 funcPhiS = 1.0-0.011905/varShat4; 729 funcPhiS = 1.0-0.011905/varShat4; 755 if (varShat<1.9156) { 730 if (varShat<1.9156) { 756 const G4double dum0 = -0.160723 731 const G4double dum0 = -0.160723 + 3.755030*varShat 757 -1.798138*varSha 732 -1.798138*varShat2 + 0.672827*varShat3 758 -0.120772*varSha 733 -0.120772*varShat4; 759 funcGS = std::tanh(dum0); 734 funcGS = std::tanh(dum0); 760 } else { 735 } else { 761 funcGS = 1.0-0.023065/varShat4; 736 funcGS = 1.0-0.023065/varShat4; 762 } 737 } 763 } 738 } 764 } 739 } 765 } 740 } 766 741 767 // s goes up to 2 with ds = 0.01 to be the def 742 // s goes up to 2 with ds = 0.01 to be the default bining 768 void G4eBremsstrahlungRelModel::InitLPMFunctio 743 void G4eBremsstrahlungRelModel::InitLPMFunctions() 769 { 744 { 770 if (!fLPMFuncs->fIsInitialized) { << 745 if (!gLPMFuncs.fIsInitialized) { 771 const G4int num = fLPMFuncs->fSLimit*fLPMF << 746 const G4int num = gLPMFuncs.fSLimit*gLPMFuncs.fISDelta+1; 772 fLPMFuncs->fLPMFuncG.resize(num); << 747 gLPMFuncs.fLPMFuncG.resize(num); 773 fLPMFuncs->fLPMFuncPhi.resize(num); << 748 gLPMFuncs.fLPMFuncPhi.resize(num); 774 for (G4int i = 0; i < num; ++i) { 749 for (G4int i = 0; i < num; ++i) { 775 const G4double sval=i/fLPMFuncs->fISDelt << 750 const G4double sval=i/gLPMFuncs.fISDelta; 776 ComputeLPMGsPhis(fLPMFuncs->fLPMFuncG[i] << 751 ComputeLPMGsPhis(gLPMFuncs.fLPMFuncG[i],gLPMFuncs.fLPMFuncPhi[i],sval); 777 } 752 } 778 fLPMFuncs->fIsInitialized = true; << 753 gLPMFuncs.fIsInitialized = true; 779 } 754 } 780 } 755 } 781 756 782 void G4eBremsstrahlungRelModel::GetLPMFunction 757 void G4eBremsstrahlungRelModel::GetLPMFunctions(G4double& lpmGs, 783 758 G4double& lpmPhis, 784 759 const G4double sval) 785 { 760 { 786 if (sval < fLPMFuncs->fSLimit) { << 761 if (sval < gLPMFuncs.fSLimit) { 787 G4double val = sval*fLPMFuncs->fISDelt << 762 G4double val = sval*gLPMFuncs.fISDelta; 788 const G4int ilow = (G4int)val; 763 const G4int ilow = (G4int)val; 789 val -= ilow; 764 val -= ilow; 790 lpmGs = (fLPMFuncs->fLPMFuncG[ilow+1]-fL << 765 lpmGs = (gLPMFuncs.fLPMFuncG[ilow+1]-gLPMFuncs.fLPMFuncG[ilow])*val 791 + fLPMFuncs->fLPMFuncG[ilow]; << 766 + gLPMFuncs.fLPMFuncG[ilow]; 792 lpmPhis = (fLPMFuncs->fLPMFuncPhi[ilow+1]- << 767 lpmPhis = (gLPMFuncs.fLPMFuncPhi[ilow+1]-gLPMFuncs.fLPMFuncPhi[ilow])*val 793 + fLPMFuncs->fLPMFuncPhi[ilow]; << 768 + gLPMFuncs.fLPMFuncPhi[ilow]; 794 } else { 769 } else { 795 G4double ss = sval*sval; 770 G4double ss = sval*sval; 796 ss *= ss; 771 ss *= ss; 797 lpmPhis = 1.0-0.01190476/ss; 772 lpmPhis = 1.0-0.01190476/ss; 798 lpmGs = 1.0-0.0230655/ss; 773 lpmGs = 1.0-0.0230655/ss; 799 } 774 } 800 } 775 } 801 776 802 777