Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc (Version 10.5)


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