Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4WentzelOKandVIxSection.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/G4WentzelOKandVIxSection.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4WentzelOKandVIxSection.cc (Version 10.4)


  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: G4WentzelOKandVIxSection.cc 105734 2017-08-16 12:58:28Z gcosmo $
 26 //                                                 27 //
 27 // -------------------------------------------     28 // -------------------------------------------------------------------
 28 //                                                 29 //
 29 // GEANT4 Class file                               30 // GEANT4 Class file
 30 //                                                 31 //
 31 //                                                 32 //
 32 // File name:   G4WentzelOKandVIxSection           33 // File name:   G4WentzelOKandVIxSection
 33 //                                                 34 //
 34 // Author:      V.Ivanchenko                       35 // Author:      V.Ivanchenko 
 35 //                                                 36 //
 36 // Creation date: 09.04.2008 from G4MuMscModel     37 // Creation date: 09.04.2008 from G4MuMscModel
 37 //                                                 38 //
 38 // Modifications:                                  39 // Modifications:
 39 //                                                 40 //
 40 // -------------------------------------------     41 // -------------------------------------------------------------------
 41 //                                                 42 //
 42                                                    43 
 43 //....oooOO0OOooo........oooOO0OOooo........oo     44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 44 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 45                                                    46 
 46 #include "G4WentzelOKandVIxSection.hh"             47 #include "G4WentzelOKandVIxSection.hh"
 47 #include "G4ScreeningMottCrossSection.hh"      << 
 48 #include "G4PhysicalConstants.hh"                  48 #include "G4PhysicalConstants.hh"
 49 #include "G4SystemOfUnits.hh"                      49 #include "G4SystemOfUnits.hh"
 50 #include "Randomize.hh"                            50 #include "Randomize.hh"
 51 #include "G4Electron.hh"                           51 #include "G4Electron.hh"
 52 #include "G4Positron.hh"                           52 #include "G4Positron.hh"
 53 #include "G4Proton.hh"                             53 #include "G4Proton.hh"
 54 #include "G4EmParameters.hh"                       54 #include "G4EmParameters.hh"
 55 #include "G4Log.hh"                                55 #include "G4Log.hh"
 56 #include "G4Exp.hh"                                56 #include "G4Exp.hh"
 57 #include "G4AutoLock.hh"                       << 
 58                                                    57 
 59 //....oooOO0OOooo........oooOO0OOooo........oo     58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 60                                                    59 
 61 G4double G4WentzelOKandVIxSection::ScreenRSqua     60 G4double G4WentzelOKandVIxSection::ScreenRSquareElec[] = {0.0};
 62 G4double G4WentzelOKandVIxSection::ScreenRSqua     61 G4double G4WentzelOKandVIxSection::ScreenRSquare[]     = {0.0};
 63 G4double G4WentzelOKandVIxSection::FormFactor[     62 G4double G4WentzelOKandVIxSection::FormFactor[]        = {0.0};
 64                                                    63 
 65 namespace                                      <<  64 #ifdef G4MULTITHREADED
 66 {                                              <<  65 G4Mutex G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex = G4MUTEX_INITIALIZER;
 67   G4Mutex theWOKVIMutex = G4MUTEX_INITIALIZER; <<  66 #endif
 68 }                                              << 
 69                                                << 
 70 const G4double alpha2 = CLHEP::fine_structure_ << 
 71 const G4double factB1= 0.5*CLHEP::pi*CLHEP::fi << 
 72 const G4double numlimit = 0.1;                 << 
 73 const G4int nwarnlimit = 50;                   << 
 74                                                    67 
 75 using namespace std;                               68 using namespace std;
 76                                                    69 
 77 G4WentzelOKandVIxSection::G4WentzelOKandVIxSec     70 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection(G4bool comb) :
 78   temp(0.,0.,0.),                                  71   temp(0.,0.,0.),
 79   isCombined(comb)                             <<  72   numlimit(0.1),
                                                   >>  73   nwarnings(0),
                                                   >>  74   nwarnlimit(50),
                                                   >>  75   isCombined(comb),
                                                   >>  76   cosThetaMax(-1.0),
                                                   >>  77   alpha2(fine_structure_const*fine_structure_const)
 80 {                                                  78 {
 81   fNistManager = G4NistManager::Instance();        79   fNistManager = G4NistManager::Instance();
 82   fG4pow = G4Pow::GetInstance();                   80   fG4pow = G4Pow::GetInstance();
 83                                                << 
 84   theElectron = G4Electron::Electron();            81   theElectron = G4Electron::Electron();
 85   thePositron = G4Positron::Positron();            82   thePositron = G4Positron::Positron();
 86   theProton   = G4Proton::Proton();                83   theProton   = G4Proton::Proton();
                                                   >>  84   lowEnergyLimit = 1.0*eV;
                                                   >>  85   G4double p0 = electron_mass_c2*classic_electr_radius;
                                                   >>  86   coeff = twopi*p0*p0;
                                                   >>  87   particle = nullptr;
                                                   >>  88 
                                                   >>  89   fNucFormfactor = fExponentialNF;
                                                   >>  90 
                                                   >>  91   currentMaterial = nullptr;
                                                   >>  92   factB = factD = formfactA = screenZ = 0.0;
                                                   >>  93   cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
                                                   >>  94 
                                                   >>  95   factB1= 0.5*CLHEP::pi*fine_structure_const;
 87                                                    96 
 88   G4double p0 = CLHEP::electron_mass_c2*CLHEP: <<  97   tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
 89   coeff = CLHEP::twopi*p0*p0;                  <<  98   ecut = etag = DBL_MAX;
                                                   >>  99   targetZ = 0;
 90   targetMass = CLHEP::proton_mass_c2;             100   targetMass = CLHEP::proton_mass_c2;
 91 }                                                 101 }
 92                                                   102 
 93 //....oooOO0OOooo........oooOO0OOooo........oo    103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 94                                                   104 
 95 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSe    105 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
 96 {                                              << 106 {}
 97   delete fMottXSection;                        << 
 98 }                                              << 
 99                                                   107 
100 //....oooOO0OOooo........oooOO0OOooo........oo    108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101                                                   109 
102 void G4WentzelOKandVIxSection::Initialise(cons    110 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p, 
103                                           G4do    111                                           G4double cosThetaLim)
104 {                                                 112 {
105   SetupParticle(p);                               113   SetupParticle(p);
106   tkin = mom2 = momCM2 = 0.0;                     114   tkin = mom2 = momCM2 = 0.0;
107   ecut = etag = DBL_MAX;                          115   ecut = etag = DBL_MAX;
108   targetZ = 0;                                    116   targetZ = 0;
109                                                   117 
110   // cosThetaMax is below 1.0 only when MSC is    118   // cosThetaMax is below 1.0 only when MSC is combined with SS
111   if(isCombined) { cosThetaMax = cosThetaLim;     119   if(isCombined) { cosThetaMax = cosThetaLim; } 
112   G4EmParameters* param = G4EmParameters::Inst << 120   
113   G4double a = param->FactorForAngleLimit()*CL << 121   G4double a = G4EmParameters::Instance()->FactorForAngleLimit()
                                                   >> 122     *CLHEP::hbarc/CLHEP::fermi;
114   factorA2 = 0.5*a*a;                             123   factorA2 = 0.5*a*a;
115   currentMaterial = nullptr;                      124   currentMaterial = nullptr;
116                                                   125 
117   fNucFormfactor = param->NuclearFormfactorTyp << 126   fNucFormfactor = G4EmParameters::Instance()->NuclearFormfactorType();
118   if(0.0 == ScreenRSquare[0]) { InitialiseA();    127   if(0.0 == ScreenRSquare[0]) { InitialiseA(); }
119                                                << 
120   // Mott corrections always added             << 
121   if((p == theElectron || p == thePositron) && << 
122     fMottXSection = new G4ScreeningMottCrossSe << 
123     fMottXSection->Initialise(p, 1.0);         << 
124   }                                            << 
125   /*                                           << 
126   G4cout << "G4WentzelOKandVIxSection::Initial << 
127    << p->GetParticleName() << " cosThetaMax= " << 
128    << "  " << ScreenRSquare[0] << " coeff= " < << 
129   */                                           << 
130 }                                                 128 }
131                                                   129 
132 //....oooOO0OOooo........oooOO0OOooo........oo    130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133                                                   131 
134 void G4WentzelOKandVIxSection::InitialiseA()      132 void G4WentzelOKandVIxSection::InitialiseA()
135 {                                                 133 {
136   // Thomas-Fermi screening radii                 134   // Thomas-Fermi screening radii
137   // Formfactors from A.V. Butkevich et al., N    135   // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
138   if(0.0 != ScreenRSquare[0]) { return; }      << 136 #ifdef G4MULTITHREADED
139   G4AutoLock l(&theWOKVIMutex);                << 137   G4MUTEXLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
                                                   >> 138 #endif
140   if(0.0 == ScreenRSquare[0]) {                   139   if(0.0 == ScreenRSquare[0]) {
141     const G4double invmev2 = 1./(CLHEP::MeV*CL << 140     G4double a0 = electron_mass_c2/0.88534; 
142     G4double a0 = CLHEP::electron_mass_c2/0.88 << 141     G4double constn = 6.937e-6/(MeV*MeV);
143     G4double constn = 6.937e-6*invmev2;        << 
144     G4double fct = G4EmParameters::Instance()-    142     G4double fct = G4EmParameters::Instance()->ScreeningFactor();
145                                                   143 
146     G4double afact = 0.5*fct*alpha2*a0*a0;     << 144     G4double afact = fct*0.5*alpha2*a0*a0;
147     ScreenRSquare[0] = afact;                     145     ScreenRSquare[0] = afact;
148     ScreenRSquare[1] = afact;                     146     ScreenRSquare[1] = afact;
149     ScreenRSquareElec[1] = afact;                 147     ScreenRSquareElec[1] = afact; 
150     FormFactor[1] = 3.097e-6*invmev2;          << 148     FormFactor[1] = constn;
151                                                   149 
152     for(G4int j=2; j<100; ++j) {                  150     for(G4int j=2; j<100; ++j) {
153       G4double x = fG4pow->Z13(j);                151       G4double x = fG4pow->Z13(j);
154       ScreenRSquare[j] = afact*(1 + G4Exp(-j*j    152       ScreenRSquare[j] = afact*(1 + G4Exp(-j*j*0.001))*x*x;
155       ScreenRSquareElec[j] = afact*x*x;           153       ScreenRSquareElec[j] = afact*x*x;
156       x = fNistManager->GetA27(j);                154       x = fNistManager->GetA27(j);
157       FormFactor[j] = constn*x*x;                 155       FormFactor[j] = constn*x*x;
158     }                                             156     } 
159   }                                            << 157   }
160   l.unlock();                                  << 158 #ifdef G4MULTITHREADED
                                                   >> 159   G4MUTEXUNLOCK(&G4WentzelOKandVIxSection::WentzelOKandVIxSectionMutex);
                                                   >> 160 #endif
                                                   >> 161   
                                                   >> 162   //G4cout << "G4WentzelOKandVIxSection::Initialise  mass= " << mass
                                                   >> 163   //         << "  " << p->GetParticleName() 
                                                   >> 164   //     << "  cosThetaMax= " << cosThetaMax << G4endl; 
                                                   >> 165   
161 }                                                 166 }
162                                                   167 
163 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164                                                   169 
165 void G4WentzelOKandVIxSection::SetupParticle(c    170 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
166 {                                                 171 {
                                                   >> 172   /*
                                                   >> 173   G4cout << "G4WentzelOKandVIxSection::SetupParticle " << p 
                                                   >> 174    << "  " << particle << "  " << this << G4endl; 
                                                   >> 175   G4cout << this << "  " << p->GetParticleName() << G4endl; 
                                                   >> 176   */
167   particle = p;                                   177   particle = p;
168   mass = particle->GetPDGMass();                  178   mass = particle->GetPDGMass();
169   spin = particle->GetPDGSpin();                  179   spin = particle->GetPDGSpin();
170   if(0.0 != spin) { spin = 0.5; }                 180   if(0.0 != spin) { spin = 0.5; }
171   G4double q = std::abs(particle->GetPDGCharge    181   G4double q = std::abs(particle->GetPDGCharge()/eplus);
172   chargeSquare = q*q;                             182   chargeSquare = q*q;
173   charge3 = chargeSquare*q;                       183   charge3 = chargeSquare*q;
174   tkin = 0.0;                                     184   tkin = 0.0;
175   currentMaterial = nullptr;                      185   currentMaterial = nullptr;
176   targetZ = 0;                                    186   targetZ = 0;
177 }                                                 187 }
178                                                   188 
179 //....oooOO0OOooo........oooOO0OOooo........oo    189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180                                                   190 
181 G4double                                          191 G4double 
182 G4WentzelOKandVIxSection::SetupKinematic(G4dou    192 G4WentzelOKandVIxSection::SetupKinematic(G4double ekin, const G4Material* mat)
183 {                                                 193 {
184   if(ekin != tkin || mat != currentMaterial) {    194   if(ekin != tkin || mat != currentMaterial) { 
185     currentMaterial = mat;                        195     currentMaterial = mat;
186     tkin  = ekin;                                 196     tkin  = ekin;
187     mom2  = tkin*(tkin + 2.0*mass);               197     mom2  = tkin*(tkin + 2.0*mass);
188     invbeta2 = 1.0 +  mass*mass/mom2;             198     invbeta2 = 1.0 +  mass*mass/mom2;
189     factB = spin/invbeta2;                        199     factB = spin/invbeta2; 
190     cosTetMaxNuc = isCombined ?                   200     cosTetMaxNuc = isCombined ? 
191       std::max(cosThetaMax, 1.-factorA2*mat->G    201       std::max(cosThetaMax, 1.-factorA2*mat->GetIonisation()->GetInvA23()/mom2)
192       : cosThetaMax;                              202       : cosThetaMax;
193   }                                               203   } 
194   return cosTetMaxNuc;                            204   return cosTetMaxNuc;
195 }                                                 205 }
196                                                   206 
197 //....oooOO0OOooo........oooOO0OOooo........oo    207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198                                                   208   
199 G4double                                          209 G4double
200 G4WentzelOKandVIxSection::SetupTarget(G4int Z,    210 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
201 {                                                 211 {
202   G4double cosTetMaxNuc2 = cosTetMaxNuc;          212   G4double cosTetMaxNuc2 = cosTetMaxNuc;
203   if(Z != targetZ || tkin != etag) {              213   if(Z != targetZ || tkin != etag) {
204     etag    = tkin;                               214     etag    = tkin; 
205     targetZ = std::min(Z, 99);                    215     targetZ = std::min(Z, 99);
206     G4double massT = (1 == Z) ? CLHEP::proton_    216     G4double massT = (1 == Z) ? CLHEP::proton_mass_c2 :
207       fNistManager->GetAtomicMassAmu(Z)*CLHEP:    217       fNistManager->GetAtomicMassAmu(Z)*CLHEP::amu_c2;
208     SetTargetMass(massT);                         218     SetTargetMass(massT);
209                                                   219 
210     kinFactor = coeff*Z*chargeSquare*invbeta2/    220     kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
211     if(particle == theElectron && fMottXSectio << 
212       fMottFactor = (1.0 + 2.0e-4*Z*Z);        << 
213     }                                          << 
214                                                   221 
215     if(1 == Z) {                                  222     if(1 == Z) {
216       screenZ = ScreenRSquare[targetZ]/mom2;      223       screenZ = ScreenRSquare[targetZ]/mom2;
217     } else if(mass > MeV) {                       224     } else if(mass > MeV) {
218       screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z    225       screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
219         ScreenRSquare[targetZ]/mom2;              226         ScreenRSquare[targetZ]/mom2;
220     } else {                                      227     } else {
221       G4double tau = tkin/mass;                   228       G4double tau = tkin/mass;
222       screenZ = std::min(Z*1.13,(1.13 +3.76*Z*    229       screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
223           *invbeta2*alpha2*std::sqrt(tau/(tau     230           *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
224         ScreenRSquareElec[targetZ]/mom2;          231         ScreenRSquareElec[targetZ]/mom2;
225     }                                             232     }
226     if(targetZ == 1 && particle == theProton &    233     if(targetZ == 1 && particle == theProton && cosTetMaxNuc2 < 0.0) {
227       cosTetMaxNuc2 = 0.0;                        234       cosTetMaxNuc2 = 0.0;
228     }                                             235     }
229     formfactA = FormFactor[targetZ]*mom2;         236     formfactA = FormFactor[targetZ]*mom2;
230                                                   237 
231     cosTetMaxElec = 1.0;                          238     cosTetMaxElec = 1.0;
232     ComputeMaxElectronScattering(cut);            239     ComputeMaxElectronScattering(cut); 
233   }                                               240   }
234   //G4cout << "SetupTarget:  Z= " << targetZ < << 
235   //   << " fMottFactor= " << fMottFactor << " << 
236   return cosTetMaxNuc2;                           241   return cosTetMaxNuc2;
237 }                                                 242 } 
238                                                   243 
239 //....oooOO0OOooo........oooOO0OOooo........oo    244 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240                                                   245 
241 G4double                                          246 G4double 
242 G4WentzelOKandVIxSection::ComputeTransportCros    247 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
243 {                                                 248 {
244   G4double xSection = 0.0;                        249   G4double xSection = 0.0;
245   if(cosTMax >= 1.0) { return xSection; }         250   if(cosTMax >= 1.0) { return xSection; }
246                                                   251 
247   G4double costm = std::max(cosTMax,cosTetMaxE    252   G4double costm = std::max(cosTMax,cosTetMaxElec); 
248   G4double fb = screenZ*factB;                    253   G4double fb = screenZ*factB;
249                                                   254 
250   // scattering off electrons                     255   // scattering off electrons
251   if(costm < 1.0) {                               256   if(costm < 1.0) {
252     G4double x = (1.0 - costm)/screenZ;           257     G4double x = (1.0 - costm)/screenZ;
253     if(x < numlimit) {                            258     if(x < numlimit) { 
254       G4double x2 = 0.5*x*x;                      259       G4double x2 = 0.5*x*x;
255       xSection = x2*((1.0 - 1.3333333*x + 3*x2    260       xSection = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x)); 
256     } else {                                      261     } else { 
257       G4double x1= x/(1 + x);                     262       G4double x1= x/(1 + x);
258       G4double xlog = G4Log(1.0 + x);             263       G4double xlog = G4Log(1.0 + x);  
259       xSection = xlog - x1 - fb*(x + x1 - 2*xl    264       xSection = xlog - x1 - fb*(x + x1 - 2*xlog);
260     }                                             265     }
261                                                   266 
262     if(xSection < 0.0) {                          267     if(xSection < 0.0) {
263       ++nwarnings;                                268       ++nwarnings;
264       if(nwarnings < nwarnlimit) {                269       if(nwarnings < nwarnlimit) {
265         G4cout << "G4WentzelOKandVIxSection::C    270         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
266                << " scattering on e- <0"          271                << " scattering on e- <0"
267                << G4endl;                         272                << G4endl;
268         G4cout << "cross= " << xSection           273         G4cout << "cross= " << xSection
269                << " e(MeV)= " << tkin << " p(M    274                << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2) 
270                << " Z= " << targetZ << "  "       275                << " Z= " << targetZ << "  " 
271                << particle->GetParticleName()     276                << particle->GetParticleName() << G4endl;
272         G4cout << " 1-costm= " << 1.0-costm <<    277         G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ 
273                << " x= " << x << G4endl;          278                << " x= " << x << G4endl;
274       }                                           279       }
275       xSection = 0.0;                             280       xSection = 0.0;
276     }                                             281     }
277   }                                               282   }
278   /*                                              283   /*  
279       G4cout << "G4WentzelOKandVIxSection::Com    284       G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom: \n"
280       << " Z= " << targetZ                        285       << " Z= " << targetZ 
281       << " e(MeV)= " << tkin/MeV << " XSel= "     286       << " e(MeV)= " << tkin/MeV << " XSel= " << xSection  
282       << " zmaxE= " << (1.0 - cosTetMaxElec)/s    287       << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
283       << " zmaxN= " << (1.0 - cosThetaMax)/scr    288       << " zmaxN= " << (1.0 - cosThetaMax)/screenZ 
284       << " 1-costm= " << 1.0 - cosThetaMax <<     289       << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
285   */                                              290   */
286   // scattering off nucleus                       291   // scattering off nucleus
287   if(cosTMax < 1.0) {                             292   if(cosTMax < 1.0) {
288     G4double x = (1.0 - cosTMax)/screenZ;         293     G4double x = (1.0 - cosTMax)/screenZ;
289     G4double y;                                   294     G4double y;
290     if(x < numlimit) {                            295     if(x < numlimit) { 
291       G4double x2 = 0.5*x*x;                      296       G4double x2 = 0.5*x*x;
292       y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*    297       y = x2*((1.0 - 1.3333333*x + 3*x2) - fb*x*(0.6666667 - x)); 
293     } else {                                      298     } else { 
294       G4double x1= x/(1 + x);                     299       G4double x1= x/(1 + x);
295       G4double xlog = G4Log(1.0 + x);             300       G4double xlog = G4Log(1.0 + x);  
296       y = xlog - x1 - fb*(x + x1 - 2*xlog);       301       y = xlog - x1 - fb*(x + x1 - 2*xlog); 
297     }                                             302     }
298                                                   303 
299     if(y < 0.0) {                                 304     if(y < 0.0) {
300       ++nwarnings;                                305       ++nwarnings;
301       if(nwarnings < nwarnlimit) {                306       if(nwarnings < nwarnlimit) {
302         G4cout << "G4WentzelOKandVIxSection::C    307         G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom"
303                << " scattering on nucleus <0"     308                << " scattering on nucleus <0"
304                << G4endl;                         309                << G4endl;
305         G4cout << "y= " << y                      310         G4cout << "y= " << y 
306                << " e(MeV)= " << tkin << " Z=     311                << " e(MeV)= " << tkin << " Z= " << targetZ << "  " 
307                << particle->GetParticleName()     312                << particle->GetParticleName() << G4endl;
308         G4cout << " formfactA= " << formfactA     313         G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ 
309                << " x= " << x <<G4endl;           314                << " x= " << x <<G4endl;
310       }                                           315       }
311       y = 0.0;                                    316       y = 0.0;
312     }                                             317     }
313     xSection += y*targetZ;                        318     xSection += y*targetZ; 
314   }                                               319   }
315   xSection *= kinFactor;                          320   xSection *= kinFactor;
316                                                   321  
317   /*                                              322   /* 
318   G4cout << "Z= " << targetZ << " XStot= " <<     323   G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
319          << " screenZ= " << screenZ << " formF    324          << " screenZ= " << screenZ << " formF= " << formfactA 
320          << " for " << particle->GetParticleNa    325          << " for " << particle->GetParticleName() 
321    << " m= " << mass << " 1/v= " << sqrt(invbe    326    << " m= " << mass << " 1/v= " << sqrt(invbeta2) 
322    << " p= " << sqrt(mom2)                        327    << " p= " << sqrt(mom2)
323          << " x= " << x << G4endl;                328          << " x= " << x << G4endl;
324   */                                              329   */
325   return xSection;                                330   return xSection; 
326 }                                                 331 }
327                                                   332 
328 //....oooOO0OOooo........oooOO0OOooo........oo    333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329                                                   334 
330 G4ThreeVector&                                    335 G4ThreeVector&
331 G4WentzelOKandVIxSection::SampleSingleScatteri    336 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
332                                                   337                                                  G4double cosTMax,
333                                                   338                                                  G4double elecRatio)
334 {                                                 339 {
335   temp.set(0.0,0.0,1.0);                          340   temp.set(0.0,0.0,1.0);
336   CLHEP::HepRandomEngine* rndmEngineMod = G4Ra    341   CLHEP::HepRandomEngine* rndmEngineMod = G4Random::getTheEngine();
337                                                   342  
338   G4double formf = formfactA;                     343   G4double formf = formfactA;
339   G4double cost1 = cosTMin;                       344   G4double cost1 = cosTMin;
340   G4double cost2 = cosTMax;                       345   G4double cost2 = cosTMax;
341   if(elecRatio > 0.0) {                           346   if(elecRatio > 0.0) {
342     if(rndmEngineMod->flat() <= elecRatio) {      347     if(rndmEngineMod->flat() <= elecRatio) {
343       formf = 0.0;                                348       formf = 0.0;
344       cost1 = std::max(cost1,cosTetMaxElec);      349       cost1 = std::max(cost1,cosTetMaxElec);
345       cost2 = std::max(cost2,cosTetMaxElec);      350       cost2 = std::max(cost2,cosTetMaxElec);
346     }                                             351     }
347   }                                               352   }
348   if(cost1 > cost2) {                             353   if(cost1 > cost2) {
349     G4double w1 = 1. - cost1;                  << 
350     G4double w2 = 1. - cost2;                  << 
351     G4double w3 = rndmEngineMod->flat()*(w2 -  << 
352     G4double z1 = ((w2 - w3)*screenZ + w1*w2)/ << 
353     G4double fm = 1.0;                         << 
354                                                   354 
                                                   >> 355     G4double w1 = 1. - cost1 + screenZ;
                                                   >> 356     G4double w2 = 1. - cost2 + screenZ;
                                                   >> 357     G4double z1 = w1*w2/(w1 + rndmEngineMod->flat()*(w2 - w1)) - screenZ;
                                                   >> 358 
                                                   >> 359     G4double fm = 1.0;
355     if(fNucFormfactor == fExponentialNF) {        360     if(fNucFormfactor == fExponentialNF) {
356       fm += formf*z1;                             361       fm += formf*z1;
357       fm = 1.0/(fm*fm);                           362       fm = 1.0/(fm*fm);
358     } else if(fNucFormfactor == fGaussianNF) {    363     } else if(fNucFormfactor == fGaussianNF) {
359       fm = G4Exp(-2*formf*z1);                    364       fm = G4Exp(-2*formf*z1);
360     } else if(fNucFormfactor == fFlatNF) {        365     } else if(fNucFormfactor == fFlatNF) {
361       static const G4double ccoef = 0.00508/CL << 366       static const G4double ccoef = 0.00508/MeV;
362       G4double x = std::sqrt(2.*mom2*z1)*ccoef    367       G4double x = std::sqrt(2.*mom2*z1)*ccoef*2.;
363       fm = FlatFormfactor(x);                     368       fm = FlatFormfactor(x);
364       fm *= FlatFormfactor(x*0.6*fG4pow->A13(f << 369       fm *= FlatFormfactor(x*0.6
                                                   >> 370       *fG4pow->A13(fNistManager->GetAtomicMassAmu(targetZ)));
365     }                                             371     }
366     // G4cout << " fm=" << fm << "  " << fMott << 372 
367     G4double grej;                             << 373     G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2. - z1))
368     if(nullptr != fMottXSection) {             << 374       *fm*fm/(1.0 + z1*factD);
369       fMottXSection->SetupKinematic(tkin, targ << 375 
370       grej = fMottXSection->RatioMottRutherfor << 376     if(rndmEngineMod->flat() <= grej ) {
371     } else {                                   << 
372       grej = (1. - z1*factB + factB1*targetZ*s << 
373       *fm/(1.0 + z1*factD);                    << 
374     }                                          << 
375     if(fMottFactor*rndmEngineMod->flat() <= gr << 
376       // exclude "false" scattering due to for    377       // exclude "false" scattering due to formfactor and spin effect
377       G4double cost = 1.0 - z1;                   378       G4double cost = 1.0 - z1;
378       if(cost > 1.0)       { cost = 1.0; }        379       if(cost > 1.0)       { cost = 1.0; }
379       else if(cost < -1.0) { cost =-1.0; }        380       else if(cost < -1.0) { cost =-1.0; }
380       G4double sint = sqrt((1.0 - cost)*(1.0 +    381       G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
381       //G4cout << "sint= " << sint << G4endl;     382       //G4cout << "sint= " << sint << G4endl;
382       G4double phi  = twopi*rndmEngineMod->fla    383       G4double phi  = twopi*rndmEngineMod->flat();
383       temp.set(sint*cos(phi),sint*sin(phi),cos    384       temp.set(sint*cos(phi),sint*sin(phi),cost);
384     }                                             385     }
385   }                                               386   }
386   return temp;                                    387   return temp;
387 }                                                 388 }
388                                                   389 
389 //....oooOO0OOooo........oooOO0OOooo........oo    390 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390                                                   391 
391 void                                              392 void 
392 G4WentzelOKandVIxSection::ComputeMaxElectronSc    393 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
393 {                                                 394 {
394   if(mass > MeV) {                                395   if(mass > MeV) {
395     G4double ratio = electron_mass_c2/mass;       396     G4double ratio = electron_mass_c2/mass;
396     G4double tau = tkin/mass;                     397     G4double tau = tkin/mass;
397     G4double tmax = 2.0*electron_mass_c2*tau*(    398     G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
398       (1.0 + 2.0*ratio*(tau + 1.0) + ratio*rat    399       (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
399     cosTetMaxElec = 1.0 - std::min(cutEnergy,     400     cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
400   } else {                                        401   } else {
401                                                   402 
402     G4double tmax = (particle == theElectron)     403     G4double tmax = (particle == theElectron) ? 0.5*tkin : tkin;
403     G4double t = std::min(cutEnergy, tmax);       404     G4double t = std::min(cutEnergy, tmax);
404     G4double mom21 = t*(t + 2.0*electron_mass_    405     G4double mom21 = t*(t + 2.0*electron_mass_c2);
405     G4double t1 = tkin - t;                       406     G4double t1 = tkin - t;
406     //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax    407     //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= " 
407     //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4end    408     //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
408     if(t1 > 0.0) {                                409     if(t1 > 0.0) {
409       G4double mom22 = t1*(t1 + 2.0*mass);        410       G4double mom22 = t1*(t1 + 2.0*mass);
410       G4double ctm = (mom2 + mom22 - mom21)*0.    411       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
411       if(ctm <  1.0) { cosTetMaxElec = ctm; }     412       if(ctm <  1.0) { cosTetMaxElec = ctm; }
412       if(particle == theElectron && cosTetMaxE    413       if(particle == theElectron && cosTetMaxElec < 0.0) { 
413         cosTetMaxElec = 0.0;                      414         cosTetMaxElec = 0.0; 
414       }                                           415       }
415     }                                             416     }
416   }                                               417   }
417 }                                                 418 }
418                                                   419 
419 //....oooOO0OOooo........oooOO0OOooo........oo    420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
420                                                   421 
421 G4double                                          422 G4double 
422 G4WentzelOKandVIxSection::ComputeSecondTranspo    423 G4WentzelOKandVIxSection::ComputeSecondTransportMoment(G4double /*CosThetaMax*/)
423 {                                                 424 {
424   return 0.0;                                     425   return 0.0;
425 }                                                 426 }
426                                                   427 
427 //....oooOO0OOooo........oooOO0OOooo........oo    428 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
428                                                   429