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