| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4WentzelOKandVIxSection.cc 74309 2013-10-03 06:42:30Z 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 //
>> 41 //
>> 42
40 // ------------------------------------------- 43 // -------------------------------------------------------------------
41 // 44 //
42 45
43 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
44 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
45 48
46 #include "G4WentzelOKandVIxSection.hh" 49 #include "G4WentzelOKandVIxSection.hh"
47 #include "G4ScreeningMottCrossSection.hh" <<
48 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh" 52 #include "Randomize.hh"
51 #include "G4Electron.hh" 53 #include "G4Electron.hh"
52 #include "G4Positron.hh" 54 #include "G4Positron.hh"
53 #include "G4Proton.hh" 55 #include "G4Proton.hh"
54 #include "G4EmParameters.hh" << 56 #include "G4LossTableManager.hh"
55 #include "G4Log.hh" 57 #include "G4Log.hh"
56 #include "G4Exp.hh" 58 #include "G4Exp.hh"
57 #include "G4AutoLock.hh" <<
58 59
59 //....oooOO0OOooo........oooOO0OOooo........oo 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 61
61 G4double G4WentzelOKandVIxSection::ScreenRSqua 62 G4double G4WentzelOKandVIxSection::ScreenRSquareElec[] = {0.0};
62 G4double G4WentzelOKandVIxSection::ScreenRSqua 63 G4double G4WentzelOKandVIxSection::ScreenRSquare[] = {0.0};
63 G4double G4WentzelOKandVIxSection::FormFactor[ 64 G4double G4WentzelOKandVIxSection::FormFactor[] = {0.0};
64 65
65 namespace <<
66 { <<
67 G4Mutex theWOKVIMutex = G4MUTEX_INITIALIZER; <<
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 <<
75 using namespace std; 66 using namespace std;
76 67
77 G4WentzelOKandVIxSection::G4WentzelOKandVIxSec << 68 G4WentzelOKandVIxSection::G4WentzelOKandVIxSection() :
78 temp(0.,0.,0.), << 69 numlimit(0.1),
79 isCombined(comb) << 70 nwarnings(0),
>> 71 nwarnlimit(50),
>> 72 alpha2(fine_structure_const*fine_structure_const)
80 { 73 {
81 fNistManager = G4NistManager::Instance(); 74 fNistManager = G4NistManager::Instance();
82 fG4pow = G4Pow::GetInstance(); 75 fG4pow = G4Pow::GetInstance();
83 <<
84 theElectron = G4Electron::Electron(); 76 theElectron = G4Electron::Electron();
85 thePositron = G4Positron::Positron(); 77 thePositron = G4Positron::Positron();
86 theProton = G4Proton::Proton(); 78 theProton = G4Proton::Proton();
>> 79 lowEnergyLimit = 1.0*eV;
>> 80 G4double p0 = electron_mass_c2*classic_electr_radius;
>> 81 coeff = twopi*p0*p0;
>> 82 particle = 0;
>> 83
>> 84 // Thomas-Fermi screening radii
>> 85 // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
>> 86
>> 87 if(0.0 == ScreenRSquare[0]) {
>> 88 G4double a0 = electron_mass_c2/0.88534;
>> 89 G4double constn = 6.937e-6/(MeV*MeV);
>> 90
>> 91 ScreenRSquare[0] = alpha2*a0*a0;
>> 92 ScreenRSquareElec[0] = ScreenRSquare[0];
>> 93 for(G4int j=1; j<100; ++j) {
>> 94 G4double x = a0*fG4pow->Z13(j);
>> 95 if(1 == j) { ScreenRSquare[j] = 0.5*alpha2*a0*a0; }
>> 96 else {
>> 97 ScreenRSquare[j] = 0.5*(1 + G4Exp(-j*j*0.001))*alpha2*x*x;
>> 98 ScreenRSquareElec[j] = 0.5*alpha2*x*x;
>> 99 }
>> 100 x = fNistManager->GetA27(j);
>> 101 FormFactor[j] = constn*x*x;
>> 102 }
>> 103 }
>> 104 currentMaterial = 0;
>> 105 elecXSRatio = factB = factD = formfactA = screenZ = 0.0;
>> 106 cosTetMaxElec = cosTetMaxNuc = invbeta2 = kinFactor = gam0pcmp = pcmp2 = 1.0;
87 107
88 G4double p0 = CLHEP::electron_mass_c2*CLHEP: << 108 factB1= 0.5*CLHEP::pi*fine_structure_const;
89 coeff = CLHEP::twopi*p0*p0; << 109
90 targetMass = CLHEP::proton_mass_c2; << 110 tkin = mom2 = momCM2 = factorA2 = mass = spin = chargeSquare = charge3 = 0.0;
>> 111 ecut = etag = DBL_MAX;
>> 112 targetZ = 0;
>> 113 cosThetaMax = 1.0;
>> 114 targetMass = proton_mass_c2;
>> 115 particle = 0;
91 } 116 }
92 117
93 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 119
95 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSe 120 G4WentzelOKandVIxSection::~G4WentzelOKandVIxSection()
96 { << 121 {}
97 delete fMottXSection; <<
98 } <<
99 122
100 //....oooOO0OOooo........oooOO0OOooo........oo 123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 124
102 void G4WentzelOKandVIxSection::Initialise(cons 125 void G4WentzelOKandVIxSection::Initialise(const G4ParticleDefinition* p,
103 G4do << 126 G4double CosThetaLim)
104 { 127 {
105 SetupParticle(p); 128 SetupParticle(p);
106 tkin = mom2 = momCM2 = 0.0; 129 tkin = mom2 = momCM2 = 0.0;
107 ecut = etag = DBL_MAX; 130 ecut = etag = DBL_MAX;
108 targetZ = 0; 131 targetZ = 0;
109 << 132 cosThetaMax = CosThetaLim;
110 // cosThetaMax is below 1.0 only when MSC is << 133 G4double a = G4LossTableManager::Instance()->FactorForAngleLimit()
111 if(isCombined) { cosThetaMax = cosThetaLim; << 134 *CLHEP::hbarc/CLHEP::fermi;
112 G4EmParameters* param = G4EmParameters::Inst <<
113 G4double a = param->FactorForAngleLimit()*CL <<
114 factorA2 = 0.5*a*a; 135 factorA2 = 0.5*a*a;
115 currentMaterial = nullptr; << 136 currentMaterial = 0;
116 << 137
117 fNucFormfactor = param->NuclearFormfactorTyp << 138 //G4cout << "G4WentzelOKandVIxSection::Initialise mass= " << mass
118 if(0.0 == ScreenRSquare[0]) { InitialiseA(); << 139 // << " " << p->GetParticleName() << G4endl;
119 << 140
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 } <<
131 <<
132 //....oooOO0OOooo........oooOO0OOooo........oo <<
133 <<
134 void G4WentzelOKandVIxSection::InitialiseA() <<
135 { <<
136 // Thomas-Fermi screening radii <<
137 // Formfactors from A.V. Butkevich et al., N <<
138 if(0.0 != ScreenRSquare[0]) { return; } <<
139 G4AutoLock l(&theWOKVIMutex); <<
140 if(0.0 == ScreenRSquare[0]) { <<
141 const G4double invmev2 = 1./(CLHEP::MeV*CL <<
142 G4double a0 = CLHEP::electron_mass_c2/0.88 <<
143 G4double constn = 6.937e-6*invmev2; <<
144 G4double fct = G4EmParameters::Instance()- <<
145 <<
146 G4double afact = 0.5*fct*alpha2*a0*a0; <<
147 ScreenRSquare[0] = afact; <<
148 ScreenRSquare[1] = afact; <<
149 ScreenRSquareElec[1] = afact; <<
150 FormFactor[1] = 3.097e-6*invmev2; <<
151 <<
152 for(G4int j=2; j<100; ++j) { <<
153 G4double x = fG4pow->Z13(j); <<
154 ScreenRSquare[j] = afact*(1 + G4Exp(-j*j <<
155 ScreenRSquareElec[j] = afact*x*x; <<
156 x = fNistManager->GetA27(j); <<
157 FormFactor[j] = constn*x*x; <<
158 } <<
159 } <<
160 l.unlock(); <<
161 } 141 }
162 142
163 //....oooOO0OOooo........oooOO0OOooo........oo 143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
164 144
165 void G4WentzelOKandVIxSection::SetupParticle(c 145 void G4WentzelOKandVIxSection::SetupParticle(const G4ParticleDefinition* p)
166 { 146 {
167 particle = p; 147 particle = p;
168 mass = particle->GetPDGMass(); 148 mass = particle->GetPDGMass();
169 spin = particle->GetPDGSpin(); 149 spin = particle->GetPDGSpin();
170 if(0.0 != spin) { spin = 0.5; } 150 if(0.0 != spin) { spin = 0.5; }
171 G4double q = std::abs(particle->GetPDGCharge << 151 G4double q = std::fabs(particle->GetPDGCharge()/eplus);
172 chargeSquare = q*q; 152 chargeSquare = q*q;
173 charge3 = chargeSquare*q; 153 charge3 = chargeSquare*q;
174 tkin = 0.0; 154 tkin = 0.0;
175 currentMaterial = nullptr; << 155 currentMaterial = 0;
176 targetZ = 0; 156 targetZ = 0;
177 } 157 }
178 158
179 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180 <<
181 G4double <<
182 G4WentzelOKandVIxSection::SetupKinematic(G4dou <<
183 { <<
184 if(ekin != tkin || mat != currentMaterial) { <<
185 currentMaterial = mat; <<
186 tkin = ekin; <<
187 mom2 = tkin*(tkin + 2.0*mass); <<
188 invbeta2 = 1.0 + mass*mass/mom2; <<
189 factB = spin/invbeta2; <<
190 cosTetMaxNuc = isCombined ? <<
191 std::max(cosThetaMax, 1.-factorA2*mat->G <<
192 : cosThetaMax; <<
193 } <<
194 return cosTetMaxNuc; <<
195 } <<
196 <<
197 //....oooOO0OOooo........oooOO0OOooo........oo <<
198 160
199 G4double 161 G4double
200 G4WentzelOKandVIxSection::SetupTarget(G4int Z, 162 G4WentzelOKandVIxSection::SetupTarget(G4int Z, G4double cut)
201 { 163 {
202 G4double cosTetMaxNuc2 = cosTetMaxNuc; 164 G4double cosTetMaxNuc2 = cosTetMaxNuc;
203 if(Z != targetZ || tkin != etag) { 165 if(Z != targetZ || tkin != etag) {
204 etag = tkin; 166 etag = tkin;
205 targetZ = std::min(Z, 99); << 167 targetZ = Z;
206 G4double massT = (1 == Z) ? CLHEP::proton_ << 168 if(targetZ > 99) { targetZ = 99; }
207 fNistManager->GetAtomicMassAmu(Z)*CLHEP: << 169 SetTargetMass(fNistManager->GetAtomicMassAmu(targetZ)*CLHEP::amu_c2);
208 SetTargetMass(massT); << 170 //G4double tmass2 = targetMass*targetMass;
>> 171 //G4double etot = tkin + mass;
>> 172 //G4double invmass2 = mass*mass + tmass2 + 2*etot*targetMass;
>> 173 //momCM2 = mom2*tmass2/invmass2;
>> 174 //gam0pcmp = (etot + targetMass)*targetMass/invmass2;
>> 175 //pcmp2 = tmass2/invmass2;
209 176
210 kinFactor = coeff*Z*chargeSquare*invbeta2/ 177 kinFactor = coeff*Z*chargeSquare*invbeta2/mom2;
211 if(particle == theElectron && fMottXSectio <<
212 fMottFactor = (1.0 + 2.0e-4*Z*Z); <<
213 } <<
214 178
215 if(1 == Z) { 179 if(1 == Z) {
216 screenZ = ScreenRSquare[targetZ]/mom2; 180 screenZ = ScreenRSquare[targetZ]/mom2;
217 } else if(mass > MeV) { 181 } else if(mass > MeV) {
218 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z 182 screenZ = std::min(Z*1.13,1.13 +3.76*Z*Z*invbeta2*alpha2*chargeSquare)*
219 ScreenRSquare[targetZ]/mom2; << 183 ScreenRSquare[targetZ]/mom2;
220 } else { 184 } else {
221 G4double tau = tkin/mass; 185 G4double tau = tkin/mass;
222 screenZ = std::min(Z*1.13,(1.13 +3.76*Z* 186 screenZ = std::min(Z*1.13,(1.13 +3.76*Z*Z
223 *invbeta2*alpha2*std::sqrt(tau/(tau << 187 *invbeta2*alpha2*std::sqrt(tau/(tau + fG4pow->Z23(targetZ)))))*
224 ScreenRSquareElec[targetZ]/mom2; << 188 ScreenRSquareElec[targetZ]/mom2;
225 } 189 }
226 if(targetZ == 1 && particle == theProton & << 190 if(targetZ == 1 && cosTetMaxNuc2 < 0.0 && particle == theProton) {
227 cosTetMaxNuc2 = 0.0; 191 cosTetMaxNuc2 = 0.0;
228 } 192 }
229 formfactA = FormFactor[targetZ]*mom2; 193 formfactA = FormFactor[targetZ]*mom2;
230 194
231 cosTetMaxElec = 1.0; 195 cosTetMaxElec = 1.0;
232 ComputeMaxElectronScattering(cut); 196 ComputeMaxElectronScattering(cut);
233 } 197 }
234 //G4cout << "SetupTarget: Z= " << targetZ < <<
235 // << " fMottFactor= " << fMottFactor << " <<
236 return cosTetMaxNuc2; 198 return cosTetMaxNuc2;
237 } 199 }
238 200
239 //....oooOO0OOooo........oooOO0OOooo........oo 201 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
240 202
241 G4double 203 G4double
242 G4WentzelOKandVIxSection::ComputeTransportCros 204 G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom(G4double cosTMax)
243 { 205 {
>> 206 G4double xsec = 0.0;
>> 207 if(cosTMax >= 1.0) { return xsec; }
>> 208
244 G4double xSection = 0.0; 209 G4double xSection = 0.0;
245 if(cosTMax >= 1.0) { return xSection; } << 210 G4double x = 0;
>> 211 G4double y = 0;
>> 212 G4double x1= 0;
>> 213 G4double x2= 0;
>> 214 G4double xlog = 0.0;
246 215
247 G4double costm = std::max(cosTMax,cosTetMaxE 216 G4double costm = std::max(cosTMax,cosTetMaxElec);
248 G4double fb = screenZ*factB; 217 G4double fb = screenZ*factB;
249 218
250 // scattering off electrons 219 // scattering off electrons
251 if(costm < 1.0) { 220 if(costm < 1.0) {
252 G4double x = (1.0 - costm)/screenZ; << 221 x = (1.0 - costm)/screenZ;
253 if(x < numlimit) { 222 if(x < numlimit) {
254 G4double x2 = 0.5*x*x; << 223 x2 = 0.5*x*x;
255 xSection = x2*((1.0 - 1.3333333*x + 3*x2 << 224 y = x2*(1.0 - 1.3333333*x + 3*x2);
>> 225 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
256 } else { 226 } else {
257 G4double x1= x/(1 + x); << 227 x1= x/(1 + x);
258 G4double xlog = G4Log(1.0 + x); << 228 xlog = G4Log(1.0 + x);
259 xSection = xlog - x1 - fb*(x + x1 - 2*xl << 229 y = xlog - x1;
>> 230 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
260 } 231 }
261 232
262 if(xSection < 0.0) { << 233 if(y < 0.0) {
263 ++nwarnings; 234 ++nwarnings;
264 if(nwarnings < nwarnlimit) { 235 if(nwarnings < nwarnlimit) {
265 G4cout << "G4WentzelOKandVIxSection::C << 236 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
266 << " scattering on e- <0" << 237 << G4endl;
267 << G4endl; << 238 G4cout << "y= " << y
268 G4cout << "cross= " << xSection << 239 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
269 << " e(MeV)= " << tkin << " p(M << 240 << " Z= " << targetZ << " "
270 << " Z= " << targetZ << " " << 241 << particle->GetParticleName() << G4endl;
271 << particle->GetParticleName() << 242 G4cout << " 1-costm= " << 1.0-costm << " screenZ= " << screenZ
272 G4cout << " 1-costm= " << 1.0-costm << << 243 << " x= " << x << G4endl;
273 << " x= " << x << G4endl; <<
274 } 244 }
275 xSection = 0.0; << 245 y = 0.0;
276 } 246 }
>> 247 xSection = y;
277 } 248 }
278 /* << 249 /*
279 G4cout << "G4WentzelOKandVIxSection::Com << 250 G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
280 << " Z= " << targetZ << 251 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
281 << " e(MeV)= " << tkin/MeV << " XSel= " << 252 << " cut(MeV)= " << ecut/MeV
282 << " zmaxE= " << (1.0 - cosTetMaxElec)/s << 253 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
283 << " zmaxN= " << (1.0 - cosThetaMax)/scr << 254 << " zmaxN= " << (1.0 - cosThetaMax)/screenZ
284 << " 1-costm= " << 1.0 - cosThetaMax << << 255 << " 1-costm= " << 1.0 - cosThetaMax << G4endl;
285 */ 256 */
286 // scattering off nucleus 257 // scattering off nucleus
287 if(cosTMax < 1.0) { 258 if(cosTMax < 1.0) {
288 G4double x = (1.0 - cosTMax)/screenZ; << 259 x = (1.0 - cosTMax)/screenZ;
289 G4double y; <<
290 if(x < numlimit) { 260 if(x < numlimit) {
291 G4double x2 = 0.5*x*x; << 261 x2 = 0.5*x*x;
292 y = x2*((1.0 - 1.3333333*x + 3*x2) - fb* << 262 y = x2*(1.0 - 1.3333333*x + 3*x2);
>> 263 if(0.0 < factB) { y -= fb*x2*x*(0.6666667 - x); }
293 } else { 264 } else {
294 G4double x1= x/(1 + x); << 265 x1= x/(1 + x);
295 G4double xlog = G4Log(1.0 + x); << 266 xlog = G4Log(1.0 + x);
296 y = xlog - x1 - fb*(x + x1 - 2*xlog); << 267 y = xlog - x1;
>> 268 if(0.0 < factB) { y -= fb*(x + x1 - 2*xlog); }
297 } 269 }
298 270
299 if(y < 0.0) { 271 if(y < 0.0) {
300 ++nwarnings; 272 ++nwarnings;
301 if(nwarnings < nwarnlimit) { 273 if(nwarnings < nwarnlimit) {
302 G4cout << "G4WentzelOKandVIxSection::C << 274 G4cout << "G4WentzelOKandVIxSection::ComputeTransportCrossSectionPerAtom scattering on e- <0"
303 << " scattering on nucleus <0" << 275 << G4endl;
304 << G4endl; << 276 G4cout << "y= " << y
305 G4cout << "y= " << y << 277 << " e(MeV)= " << tkin << " Z= " << targetZ << " "
306 << " e(MeV)= " << tkin << " Z= << 278 << particle->GetParticleName() << G4endl;
307 << particle->GetParticleName() << 279 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
308 G4cout << " formfactA= " << formfactA << 280 << " x= " << " x1= " << x1 <<G4endl;
309 << " x= " << x <<G4endl; <<
310 } 281 }
311 y = 0.0; 282 y = 0.0;
312 } 283 }
313 xSection += y*targetZ; 284 xSection += y*targetZ;
314 } 285 }
315 xSection *= kinFactor; 286 xSection *= kinFactor;
316 << 287 /*
317 /* <<
318 G4cout << "Z= " << targetZ << " XStot= " << 288 G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
319 << " screenZ= " << screenZ << " formF << 289 << " screenZ= " << screenZ << " formF= " << formfactA
320 << " for " << particle->GetParticleNa << 290 << " for " << particle->GetParticleName()
321 << " m= " << mass << " 1/v= " << sqrt(invbe << 291 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
322 << " p= " << sqrt(mom2) << 292 << " x= " << x
323 << " x= " << x << G4endl; << 293 << G4endl;
324 */ 294 */
325 return xSection; 295 return xSection;
326 } 296 }
327 297
328 //....oooOO0OOooo........oooOO0OOooo........oo 298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 299
330 G4ThreeVector& << 300 G4ThreeVector
331 G4WentzelOKandVIxSection::SampleSingleScatteri 301 G4WentzelOKandVIxSection::SampleSingleScattering(G4double cosTMin,
332 << 302 G4double cosTMax,
333 << 303 G4double elecRatio)
334 { 304 {
335 temp.set(0.0,0.0,1.0); << 305 G4ThreeVector v(0.0,0.0,1.0);
336 CLHEP::HepRandomEngine* rndmEngineMod = G4Ra <<
337 306
338 G4double formf = formfactA; 307 G4double formf = formfactA;
339 G4double cost1 = cosTMin; 308 G4double cost1 = cosTMin;
340 G4double cost2 = cosTMax; 309 G4double cost2 = cosTMax;
341 if(elecRatio > 0.0) { 310 if(elecRatio > 0.0) {
342 if(rndmEngineMod->flat() <= elecRatio) { << 311 if(G4UniformRand() <= elecRatio) {
343 formf = 0.0; 312 formf = 0.0;
344 cost1 = std::max(cost1,cosTetMaxElec); 313 cost1 = std::max(cost1,cosTetMaxElec);
345 cost2 = std::max(cost2,cosTetMaxElec); 314 cost2 = std::max(cost2,cosTetMaxElec);
346 } 315 }
347 } 316 }
348 if(cost1 > cost2) { << 317 if(cost1 < cost2) { return v; }
349 G4double w1 = 1. - cost1; << 318
350 G4double w2 = 1. - cost2; << 319 G4double w1 = 1. - cost1 + screenZ;
351 G4double w3 = rndmEngineMod->flat()*(w2 - << 320 G4double w2 = 1. - cost2 + screenZ;
352 G4double z1 = ((w2 - w3)*screenZ + w1*w2)/ << 321 G4double z1 = w1*w2/(w1 + G4UniformRand()*(w2 - w1)) - screenZ;
353 G4double fm = 1.0; << 322
354 << 323 //G4double fm = 1.0 + formf*z1/(1.0 + (mass + tkin)*z1/targetMass);
355 if(fNucFormfactor == fExponentialNF) { << 324 G4double fm = 1.0 + formf*z1;
356 fm += formf*z1; << 325 //G4double grej = (1. - z1*factB)/( (1.0 + z1*factD)*fm*fm );
357 fm = 1.0/(fm*fm); << 326 G4double grej = (1. - z1*factB + factB1*targetZ*sqrt(z1*factB)*(2 - z1))/( (1.0 + z1*factD)*fm*fm );
358 } else if(fNucFormfactor == fGaussianNF) { << 327 // "false" scattering
359 fm = G4Exp(-2*formf*z1); << 328 if( G4UniformRand() > grej ) { return v; }
360 } else if(fNucFormfactor == fFlatNF) { << 329 // }
361 static const G4double ccoef = 0.00508/CL << 330 G4double cost = 1.0 - z1;
362 G4double x = std::sqrt(2.*mom2*z1)*ccoef << 331
363 fm = FlatFormfactor(x); << 332 if(cost > 1.0) { cost = 1.0; }
364 fm *= FlatFormfactor(x*0.6*fG4pow->A13(f << 333 else if(cost < -1.0) { cost =-1.0; }
365 } << 334 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
366 // G4cout << " fm=" << fm << " " << fMott << 335 //G4cout << "sint= " << sint << G4endl;
367 G4double grej; << 336 G4double phi = twopi*G4UniformRand();
368 if(nullptr != fMottXSection) { << 337 G4double vx1 = sint*cos(phi);
369 fMottXSection->SetupKinematic(tkin, targ << 338 G4double vy1 = sint*sin(phi);
370 grej = fMottXSection->RatioMottRutherfor << 339
371 } else { << 340 // only direction is changed
372 grej = (1. - z1*factB + factB1*targetZ*s << 341 v.set(vx1,vy1,cost);
373 *fm/(1.0 + z1*factD); << 342 return v;
374 } <<
375 if(fMottFactor*rndmEngineMod->flat() <= gr <<
376 // exclude "false" scattering due to for <<
377 G4double cost = 1.0 - z1; <<
378 if(cost > 1.0) { cost = 1.0; } <<
379 else if(cost < -1.0) { cost =-1.0; } <<
380 G4double sint = sqrt((1.0 - cost)*(1.0 + <<
381 //G4cout << "sint= " << sint << G4endl; <<
382 G4double phi = twopi*rndmEngineMod->fla <<
383 temp.set(sint*cos(phi),sint*sin(phi),cos <<
384 } <<
385 } <<
386 return temp; <<
387 } 343 }
388 344
389 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
390 346
391 void 347 void
392 G4WentzelOKandVIxSection::ComputeMaxElectronSc 348 G4WentzelOKandVIxSection::ComputeMaxElectronScattering(G4double cutEnergy)
393 { 349 {
394 if(mass > MeV) { 350 if(mass > MeV) {
395 G4double ratio = electron_mass_c2/mass; 351 G4double ratio = electron_mass_c2/mass;
396 G4double tau = tkin/mass; 352 G4double tau = tkin/mass;
397 G4double tmax = 2.0*electron_mass_c2*tau*( 353 G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
398 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*rat 354 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
>> 355 //tmax = std::min(tmax, targetZ*targetZ*10*eV);
399 cosTetMaxElec = 1.0 - std::min(cutEnergy, 356 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
400 } else { 357 } else {
401 358
402 G4double tmax = (particle == theElectron) << 359 G4double tmax = tkin;
>> 360 if(particle == theElectron) { tmax *= 0.5; }
>> 361 //tmax = std::min(tmax, targetZ*targetZ*10*eV);
403 G4double t = std::min(cutEnergy, tmax); 362 G4double t = std::min(cutEnergy, tmax);
404 G4double mom21 = t*(t + 2.0*electron_mass_ 363 G4double mom21 = t*(t + 2.0*electron_mass_c2);
405 G4double t1 = tkin - t; 364 G4double t1 = tkin - t;
406 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax 365 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
407 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4end 366 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
408 if(t1 > 0.0) { 367 if(t1 > 0.0) {
409 G4double mom22 = t1*(t1 + 2.0*mass); 368 G4double mom22 = t1*(t1 + 2.0*mass);
410 G4double ctm = (mom2 + mom22 - mom21)*0. 369 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
411 if(ctm < 1.0) { cosTetMaxElec = ctm; } 370 if(ctm < 1.0) { cosTetMaxElec = ctm; }
412 if(particle == theElectron && cosTetMaxE << 371 if(particle == theElectron && cosTetMaxElec < 0.0) { cosTetMaxElec = 0.0; }
413 cosTetMaxElec = 0.0; <<
414 } <<
415 } 372 }
416 } 373 }
417 } <<
418 <<
419 //....oooOO0OOooo........oooOO0OOooo........oo <<
420 <<
421 G4double <<
422 G4WentzelOKandVIxSection::ComputeSecondTranspo <<
423 { <<
424 return 0.0; <<
425 } 374 }
426 375
427 //....oooOO0OOooo........oooOO0OOooo........oo 376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
428 377