| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
>> 23 // $Id: G4MollerBhabhaModel.cc,v 1.15 2004/12/01 19:37:14 vnivanch Exp $
>> 24 // GEANT4 tag $Name: geant4-07-00-cand-03 $
26 // 25 //
27 // ------------------------------------------- 26 // -------------------------------------------------------------------
28 // 27 //
29 // GEANT4 Class file 28 // GEANT4 Class file
30 // 29 //
31 // 30 //
32 // File name: G4MollerBhabhaModel 31 // File name: G4MollerBhabhaModel
33 // 32 //
34 // Author: Vladimir Ivanchenko on base 33 // Author: Vladimir Ivanchenko on base of Laszlo Urban code
35 // 34 //
36 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002
37 // 36 //
38 // Modifications: 37 // Modifications:
39 // 38 //
40 // 13-11-02 Minor fix - use normalised directi 39 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko)
41 // 04-12-02 Change G4DynamicParticle construct 40 // 04-12-02 Change G4DynamicParticle constructor in PostStepDoIt (V.Ivanchenko)
42 // 23-12-02 Change interface in order to move 41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko)
43 // 27-01-03 Make models region aware (V.Ivanch 42 // 27-01-03 Make models region aware (V.Ivanchenko)
44 // 13-02-03 Add name (V.Ivanchenko) 43 // 13-02-03 Add name (V.Ivanchenko)
45 // 08-04-05 Major optimisation of internal int <<
46 // 25-07-05 Add protection in calculation of r <<
47 // of complete energy transfer from e <<
48 // 06-02-06 ComputeCrossSectionPerElectron, Co <<
49 // 15-05-06 Fix MinEnergyCut (V.Ivanchenko) <<
50 // 44 //
51 // 45 //
52 // Class Description: 46 // Class Description:
53 // 47 //
54 // Implementation of energy loss and delta-ele 48 // Implementation of energy loss and delta-electron production by e+/e-
55 // 49 //
56 // ------------------------------------------- 50 // -------------------------------------------------------------------
57 // 51 //
58 //....oooOO0OOooo........oooOO0OOooo........oo << 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
59 //....oooOO0OOooo........oooOO0OOooo........oo << 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60 54
61 #include "G4MollerBhabhaModel.hh" 55 #include "G4MollerBhabhaModel.hh"
62 #include "G4PhysicalConstants.hh" <<
63 #include "G4SystemOfUnits.hh" <<
64 #include "G4Electron.hh" 56 #include "G4Electron.hh"
65 #include "G4Positron.hh" 57 #include "G4Positron.hh"
66 #include "Randomize.hh" 58 #include "Randomize.hh"
67 #include "G4ParticleChangeForLoss.hh" <<
68 #include "G4Log.hh" <<
69 #include "G4DeltaAngle.hh" <<
70 59
71 //....oooOO0OOooo........oooOO0OOooo........oo << 60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 61
73 using namespace std; 62 using namespace std;
74 63
75 G4MollerBhabhaModel::G4MollerBhabhaModel(const 64 G4MollerBhabhaModel::G4MollerBhabhaModel(const G4ParticleDefinition* p,
76 const 65 const G4String& nam)
77 : G4VEmModel(nam), 66 : G4VEmModel(nam),
78 particle(nullptr), << 67 particle(0),
79 isElectron(true), << 68 highKinEnergy(100.*TeV),
80 twoln10(2.0*G4Log(10.0)), << 69 lowKinEnergy(0.1*keV),
81 lowLimit(0.02*keV), << 70 twoln10(2.0*log(10.0)),
82 isInitialised(false) << 71 lowLimit(0.2*keV),
>> 72 isElectron(true)
83 { 73 {
>> 74 if(p) SetParticle(p);
84 theElectron = G4Electron::Electron(); 75 theElectron = G4Electron::Electron();
85 if(nullptr != p) { SetParticle(p); } <<
86 fParticleChange = nullptr; <<
87 } 76 }
88 77
89 //....oooOO0OOooo........oooOO0OOooo........oo << 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 79
91 G4MollerBhabhaModel::~G4MollerBhabhaModel() = << 80 G4MollerBhabhaModel::~G4MollerBhabhaModel()
>> 81 {}
92 82
93 //....oooOO0OOooo........oooOO0OOooo........oo << 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94 84
95 G4double G4MollerBhabhaModel::MaxSecondaryEner << 85 void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p)
96 <<
97 { 86 {
98 G4double tmax = kinEnergy; << 87 particle = p;
99 if(isElectron) { tmax *= 0.5; } << 88 if(p != theElectron) isElectron = false;
100 return tmax; <<
101 } 89 }
102 90
103 //....oooOO0OOooo........oooOO0OOooo........oo << 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104 92
105 void G4MollerBhabhaModel::Initialise(const G4P << 93 G4double G4MollerBhabhaModel::HighEnergyLimit(const G4ParticleDefinition*)
106 const G4D <<
107 { 94 {
108 if(p != particle) { SetParticle(p); } << 95 return highKinEnergy;
109 <<
110 if(isInitialised) { return; } <<
111 <<
112 isInitialised = true; <<
113 fParticleChange = GetParticleChangeForLoss() <<
114 if(UseAngularGeneratorFlag() && !GetAngularD <<
115 SetAngularDistribution(new G4DeltaAngle()) <<
116 } <<
117 } 96 }
118 97
119 //....oooOO0OOooo........oooOO0OOooo........oo << 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120 99
121 G4double G4MollerBhabhaModel::ComputeCrossSect << 100 G4double G4MollerBhabhaModel::LowEnergyLimit(const G4ParticleDefinition*)
122 const G4ParticleDefinition* p, G4doub <<
123 G4double cutEnergy, G4double maxEnergy) <<
124 { 101 {
125 if(p != particle) { SetParticle(p); } << 102 return lowKinEnergy;
126 << 103 }
127 G4double cross = 0.0; <<
128 G4double tmax = MaxSecondaryEnergy(p, kineti <<
129 tmax = std::min(maxEnergy, tmax); <<
130 //G4cout << "E= " << kineticEnergy << " cut= <<
131 // << " Emax= " << tmax << G4endl; <<
132 if(cutEnergy < tmax) { <<
133 <<
134 G4double xmin = cutEnergy/kineticEnergy; <<
135 G4double xmax = tmax/kineticEnergy; <<
136 G4double tau = kineticEnergy/electron_ma <<
137 G4double gam = tau + 1.0; <<
138 G4double gamma2= gam*gam; <<
139 G4double beta2 = tau*(tau + 2)/gamma2; <<
140 <<
141 //Moller (e-e-) scattering <<
142 if (isElectron) { <<
143 <<
144 G4double gg = (2.0*gam - 1.0)/gamma2; <<
145 cross = ((xmax - xmin)*(1.0 - gg + 1.0/( <<
146 + 1.0/((1.0-xmin <<
147 - gg*G4Log( xmax*(1.0 - xmin)/(xmi <<
148 <<
149 //Bhabha (e+e-) scattering <<
150 } else { <<
151 <<
152 G4double y = 1.0/(1.0 + gam); <<
153 G4double y2 = y*y; <<
154 G4double y12 = 1.0 - 2.0*y; <<
155 G4double b1 = 2.0 - y2; <<
156 G4double b2 = y12*(3.0 + y2); <<
157 G4double y122= y12*y12; <<
158 G4double b4 = y122*y12; <<
159 G4double b3 = b4 + y122; <<
160 104
161 cross = (xmax - xmin)*(1.0/(beta2*xmin*x << 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 - 0.5*b3*(xmin + xmax) <<
163 + b4*(xmin*xmin + xmin*xmax + xmax <<
164 - b1*G4Log(xmax/xmin); <<
165 } <<
166 106
167 cross *= twopi_mc2_rcl2/kineticEnergy; << 107 G4double G4MollerBhabhaModel::MinEnergyCut(const G4ParticleDefinition*,
168 } << 108 const G4MaterialCutsCouple* couple)
169 return cross; << 109 {
>> 110 return couple->GetMaterial()->GetIonisation()->GetMeanExcitationEnergy();
170 } 111 }
171 112
172 //....oooOO0OOooo........oooOO0OOooo........oo << 113 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 114
174 G4double G4MollerBhabhaModel::ComputeCrossSect << 115 G4bool G4MollerBhabhaModel::IsInCharge(const G4ParticleDefinition* p)
175 con <<
176 <<
177 <<
178 <<
179 <<
180 { 116 {
181 return Z*ComputeCrossSectionPerElectron(p,ki << 117 return (p == theElectron || p == G4Positron::Positron());
182 } 118 }
183 119
184 //....oooOO0OOooo........oooOO0OOooo........oo << 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185 121
186 G4double G4MollerBhabhaModel::CrossSectionPerV << 122 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p,
187 con << 123 const G4DataVector&)
188 con <<
189 <<
190 <<
191 <<
192 { 124 {
193 G4double eDensity = material->GetElectronDen << 125 if(!particle) SetParticle(p);
194 return eDensity*ComputeCrossSectionPerElectr <<
195 } 126 }
196 127
197 //....oooOO0OOooo........oooOO0OOooo........oo << 128 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 129
199 G4double G4MollerBhabhaModel::ComputeDEDXPerVo << 130 G4double G4MollerBhabhaModel::ComputeDEDX(const G4MaterialCutsCouple* couple,
200 cons <<
201 cons 131 const G4ParticleDefinition* p,
202 132 G4double kineticEnergy,
203 << 133 G4double cutEnergy)
204 { 134 {
205 if(p != particle) { SetParticle(p); } << 135 if(!particle) SetParticle(p);
206 // calculate the dE/dx due to the ionization 136 // calculate the dE/dx due to the ionization by Seltzer-Berger formula
207 // checl low-energy limit << 137 const G4Material* material = couple->GetMaterial();
208 G4double electronDensity = material->GetElec 138 G4double electronDensity = material->GetElectronDensity();
209 << 139 G4double Zeff = electronDensity/material->GetTotNbOfAtomsPerVolume();
210 G4double Zeff = material->GetIonisation()-> <<
211 G4double th = 0.25*sqrt(Zeff)*keV; 140 G4double th = 0.25*sqrt(Zeff)*keV;
212 G4double tkin = std::max(kineticEnergy, th); << 141 G4double tkin = kineticEnergy;
213 << 142 G4bool lowEnergy = false;
>> 143 if (kineticEnergy < th) {
>> 144 tkin = th;
>> 145 lowEnergy = true;
>> 146 }
214 G4double tau = tkin/electron_mass_c2; 147 G4double tau = tkin/electron_mass_c2;
215 G4double gam = tau + 1.0; 148 G4double gam = tau + 1.0;
216 G4double gamma2= gam*gam; 149 G4double gamma2= gam*gam;
217 G4double bg2 = tau*(tau + 2); << 150 G4double beta2 = 1. - 1./gamma2;
218 G4double beta2 = bg2/gamma2; << 151 G4double bg2 = beta2*gamma2;
219 152
220 G4double eexc = material->GetIonisation()-> 153 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
221 eexc /= electron_mass_c2; 154 eexc /= electron_mass_c2;
222 G4double eexc2 = eexc*eexc; 155 G4double eexc2 = eexc*eexc;
223 156
224 G4double d = std::min(cut, MaxSecondaryEnerg << 157 G4double d = min(cutEnergy, MaxSecondaryEnergy(p, tkin))/electron_mass_c2;
225 G4double dedx; 158 G4double dedx;
226 159
227 // electron 160 // electron
228 if (isElectron) { 161 if (isElectron) {
229 162
230 dedx = G4Log(2.0*(tau + 2.0)/eexc2) - 1.0 << 163 dedx = log(2.0*(tau + 2.0)/eexc2) - 1.0 - beta2
231 + G4Log((tau-d)*d) + tau/(tau-d) << 164 + log((tau-d)*d) + tau/(tau-d)
232 + (0.5*d*d + (2.0*tau + 1.)*G4Log(1. << 165 + (0.5*d*d + (2.0*tau + 1.)*log(1. - d/tau))/gamma2;
233 166
234 //positron 167 //positron
235 } else { 168 } else {
236 169
237 G4double d2 = d*d*0.5; 170 G4double d2 = d*d*0.5;
238 G4double d3 = d2*d/1.5; 171 G4double d3 = d2*d/1.5;
239 G4double d4 = d3*d*0.75; << 172 G4double d4 = d3*d*3.75;
240 G4double y = 1.0/(1.0 + gam); 173 G4double y = 1.0/(1.0 + gam);
241 dedx = G4Log(2.0*(tau + 2.0)/eexc2) + G4Lo << 174 dedx = log(2.0*(tau + 2.0)/eexc2) + log(tau*d)
242 - beta2*(tau + 2.0*d - y*(3.0*d2 175 - beta2*(tau + 2.0*d - y*(3.0*d2
243 + y*(d - d3 + y*(d2 - tau*d3 + d4)))) 176 + y*(d - d3 + y*(d2 - tau*d3 + d4))))/tau;
244 } 177 }
245 178
246 //density correction 179 //density correction
247 G4double x = G4Log(bg2)/twoln10; << 180 G4double cden = material->GetIonisation()->GetCdensity();
248 dedx -= material->GetIonisation()->DensityCo << 181 G4double mden = material->GetIonisation()->GetMdensity();
>> 182 G4double aden = material->GetIonisation()->GetAdensity();
>> 183 G4double x0den = material->GetIonisation()->GetX0density();
>> 184 G4double x1den = material->GetIonisation()->GetX1density();
>> 185 G4double x = log(bg2)/twoln10;
>> 186
>> 187 if (x >= x0den) {
>> 188 dedx -= twoln10*x - cden;
>> 189 if (x < x1den) dedx -= aden*pow(x1den-x, mden);
>> 190 }
249 191
250 // now you can compute the total ionization 192 // now you can compute the total ionization loss
251 dedx *= twopi_mc2_rcl2*electronDensity/beta2 193 dedx *= twopi_mc2_rcl2*electronDensity/beta2;
252 if (dedx < 0.0) { dedx = 0.0; } << 194 if (dedx < 0.0) dedx = 0.0;
253 << 195
254 // lowenergy extrapolation 196 // lowenergy extrapolation
255 << 197
256 if (kineticEnergy < th) { << 198 if (lowEnergy) {
257 x = kineticEnergy/th; << 199
258 if(x > 0.25) { dedx /= sqrt(x); } << 200 if (kineticEnergy >= lowLimit) dedx *= sqrt(tkin/kineticEnergy);
259 else { dedx *= 1.4*sqrt(x)/(0.1 + x); } << 201 else dedx *= sqrt(tkin*kineticEnergy)/lowLimit;
>> 202
260 } 203 }
261 return dedx; 204 return dedx;
262 } 205 }
263 206
264 //....oooOO0OOooo........oooOO0OOooo........oo << 207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
265 208
266 void << 209 G4double G4MollerBhabhaModel::CrossSection(const G4MaterialCutsCouple* couple,
267 G4MollerBhabhaModel::SampleSecondaries(std::ve << 210 const G4ParticleDefinition* p,
268 const G << 211 G4double kineticEnergy,
269 const G << 212 G4double cutEnergy,
270 G4doubl << 213 G4double maxEnergy)
271 G4doubl <<
272 { 214 {
273 G4double kineticEnergy = dp->GetKineticEnerg << 215 if(!particle) SetParticle(p);
274 //G4cout << "G4MollerBhabhaModel::SampleSeco << 216
275 // << " in " << couple->GetMaterial()- << 217 G4double cross = 0.0;
276 G4double tmax; << 218 G4double tmax = MaxSecondaryEnergy(p, kineticEnergy);
277 G4double tmin = cutEnergy; << 219 tmax = min(maxEnergy, tmax);
278 if(isElectron) { << 220
279 tmax = 0.5*kineticEnergy; << 221 if(cutEnergy < tmax) {
280 } else { << 222
281 tmax = kineticEnergy; << 223 G4double xmin = cutEnergy/kineticEnergy;
>> 224 G4double xmax = tmax/kineticEnergy;
>> 225 G4double gam = kineticEnergy/electron_mass_c2 + 1.0;
>> 226 G4double gamma2= gam*gam;
>> 227 G4double beta2 = 1.0 - 1.0/gamma2;
>> 228
>> 229 //Moller (e-e-) scattering
>> 230 if (isElectron) {
>> 231
>> 232 G4double g = (2.0*gam - 1.0)/gamma2;
>> 233 cross = ((xmax - xmin)*(1.0 - g + 1.0/(xmin*xmax)
>> 234 + 1.0/((1.0-xmin)*(1.0 - xmax)))
>> 235 - g*log( xmax*(1.0 - xmin)/(xmin*(1.0 - xmax)) ) ) / beta2;
>> 236
>> 237 //Bhabha (e+e-) scattering
>> 238 } else {
>> 239
>> 240 G4double y = 1.0/(1.0 + gam);
>> 241 G4double y2 = y*y;
>> 242 G4double y12 = 1.0 - 2.0*y;
>> 243 G4double b1 = 2.0 - y2;
>> 244 G4double b2 = y12*(3.0 + y2);
>> 245 G4double y122= y12*y12;
>> 246 G4double b4 = y122*y12;
>> 247 G4double b3 = b4 + y122;
>> 248
>> 249 cross = (xmax - xmin)*(1.0/(beta2*xmin*xmax) + b2
>> 250 - 0.5*b3*(xmin + xmax)
>> 251 + b4*(xmin*xmin + xmin*xmax + xmax*xmax)/3.0)
>> 252 - b1*log(xmax/xmin);
>> 253 }
>> 254
>> 255 cross *= twopi_mc2_rcl2*(couple->GetMaterial()->GetElectronDensity())/kineticEnergy;
282 } 256 }
283 if(maxEnergy < tmax) { tmax = maxEnergy; } << 257 return cross;
284 if(tmin >= tmax) { return; } << 258 }
285 259
>> 260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 261
>> 262 G4DynamicParticle* G4MollerBhabhaModel::SampleSecondary(
>> 263 const G4MaterialCutsCouple*,
>> 264 const G4DynamicParticle* dp,
>> 265 G4double tmin,
>> 266 G4double maxEnergy)
>> 267 {
>> 268 G4double tmax = min(maxEnergy, MaxSecondaryEnergy(dp));
>> 269 if(tmin > tmax) tmin = tmax;
>> 270
>> 271 G4double kineticEnergy = dp->GetKineticEnergy();
286 G4double energy = kineticEnergy + electron_m 272 G4double energy = kineticEnergy + electron_mass_c2;
>> 273 G4double totalMomentum = sqrt(kineticEnergy*(energy + electron_mass_c2));
287 G4double xmin = tmin/kineticEnergy; 274 G4double xmin = tmin/kineticEnergy;
288 G4double xmax = tmax/kineticEnergy; 275 G4double xmax = tmax/kineticEnergy;
289 G4double gam = energy/electron_mass_c2; 276 G4double gam = energy/electron_mass_c2;
290 G4double gamma2 = gam*gam; 277 G4double gamma2 = gam*gam;
291 G4double beta2 = 1.0 - 1.0/gamma2; 278 G4double beta2 = 1.0 - 1.0/gamma2;
292 G4double x, z, grej; << 279 G4double x, z, q, grej;
293 CLHEP::HepRandomEngine* rndmEngine = G4Rando << 280
294 G4double rndm[2]; << 281 G4ThreeVector momentum = dp->GetMomentumDirection();
295 282
296 //Moller (e-e-) scattering 283 //Moller (e-e-) scattering
297 if (isElectron) { 284 if (isElectron) {
298 285
299 G4double gg = (2.0*gam - 1.0)/gamma2; << 286 G4double g = (2.0*gam - 1.0)/gamma2;
300 G4double y = 1.0 - xmax; 287 G4double y = 1.0 - xmax;
301 grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg << 288 grej = 1.0 - g*xmax + xmax*xmax*(1.0 - g + (1.0 - g*y)/(y*y));
302 289
303 do { 290 do {
304 rndmEngine->flatArray(2, rndm); << 291 q = G4UniformRand();
305 x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm << 292 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
306 y = 1.0 - x; 293 y = 1.0 - x;
307 z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 - << 294 z = 1.0 - g*x + x*x*(1.0 - g + (1.0 - g*y)/(y*y));
308 /* 295 /*
309 if(z > grej) { 296 if(z > grej) {
310 G4cout << "G4MollerBhabhaModel::Sample 297 G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
311 << "Majorant " << grej << " < " 298 << "Majorant " << grej << " < "
312 << z << " for x= " << x 299 << z << " for x= " << x
313 << " e-e- scattering" 300 << " e-e- scattering"
314 << G4endl; 301 << G4endl;
315 } 302 }
316 */ 303 */
317 // Loop checking, 03-Aug-2015, Vladimir << 304 } while(grej * G4UniformRand() > z);
318 } while(grej * rndm[1] > z); <<
319 305
320 //Bhabha (e+e-) scattering 306 //Bhabha (e+e-) scattering
321 } else { 307 } else {
322 308
323 G4double y = 1.0/(1.0 + gam); 309 G4double y = 1.0/(1.0 + gam);
324 G4double y2 = y*y; 310 G4double y2 = y*y;
325 G4double y12 = 1.0 - 2.0*y; 311 G4double y12 = 1.0 - 2.0*y;
326 G4double b1 = 2.0 - y2; 312 G4double b1 = 2.0 - y2;
327 G4double b2 = y12*(3.0 + y2); 313 G4double b2 = y12*(3.0 + y2);
328 G4double y122= y12*y12; 314 G4double y122= y12*y12;
329 G4double b4 = y122*y12; 315 G4double b4 = y122*y12;
330 G4double b3 = b4 + y122; 316 G4double b3 = b4 + y122;
331 317
332 y = xmax*xmax; << 318 y = xmax*xmax;
333 grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + << 319 grej = -xmin*b1;
>> 320 grej += y*b2;
>> 321 grej -= xmin*xmin*xmin*b3;
>> 322 grej += y*y*b4;
>> 323 grej *= beta2;
>> 324 grej += 1.0;
334 do { 325 do {
335 rndmEngine->flatArray(2, rndm); << 326 q = G4UniformRand();
336 x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xm << 327 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
337 y = x*x; << 328 z = -x*b1;
338 z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1 << 329 y = x*x;
>> 330 z += y*b2;
>> 331 y *= x;
>> 332 z -= y*b3;
>> 333 y *= x;
>> 334 z += y*b4;
>> 335 z *= beta2;
>> 336 z += 1.0;
339 /* 337 /*
340 if(z > grej) { 338 if(z > grej) {
341 G4cout << "G4MollerBhabhaModel::Sample 339 G4cout << "G4MollerBhabhaModel::SampleSecondary Warning! "
342 << "Majorant " << grej << " < " 340 << "Majorant " << grej << " < "
343 << z << " for x= " << x 341 << z << " for x= " << x
344 << " e+e- scattering" 342 << " e+e- scattering"
345 << G4endl; 343 << G4endl;
346 } 344 }
347 */ 345 */
348 // Loop checking, 03-Aug-2015, Vladimir << 346 } while(grej * G4UniformRand() > z);
349 } while(grej * rndm[1] > z); <<
350 } 347 }
351 348
352 G4double deltaKinEnergy = x * kineticEnergy; 349 G4double deltaKinEnergy = x * kineticEnergy;
353 350
354 G4ThreeVector deltaDirection; << 351 G4double deltaMomentum =
355 << 352 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
356 if(UseAngularGeneratorFlag()) { << 353 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
357 const G4Material* mat = couple->GetMateri << 354 (deltaMomentum * totalMomentum);
358 G4int Z = SelectRandomAtomNumber(mat); << 355 G4double sint = sqrt(1.0 - cost*cost);
359 356
360 deltaDirection = << 357 G4double phi = twopi * G4UniformRand() ;
361 GetAngularDistribution()->SampleDirectio <<
362 358
363 } else { << 359 G4ThreeVector deltaDirection(sint*cos(phi),sint*sin(phi), cost) ;
364 << 360 deltaDirection.rotateUz(momentum);
365 G4double deltaMomentum = <<
366 sqrt(deltaKinEnergy * (deltaKinEnergy + <<
367 G4double cost = deltaKinEnergy * (energy + <<
368 (deltaMomentum * dp->GetTotalMomentum()) <<
369 if(cost > 1.0) { cost = 1.0; } <<
370 G4double sint = sqrt((1.0 - cost)*(1.0 + c <<
371 <<
372 G4double phi = twopi * rndmEngine->flat() <<
373 <<
374 deltaDirection.set(sint*cos(phi),sint*sin( <<
375 deltaDirection.rotateUz(dp->GetMomentumDir <<
376 } <<
377 361
378 // create G4DynamicParticle object for delta 362 // create G4DynamicParticle object for delta ray
379 auto delta = new G4DynamicParticle(theElectr << 363 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
380 vdp->push_back(delta); << 364 return delta;
>> 365 }
381 366
382 // primary change << 367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
383 kineticEnergy -= deltaKinEnergy; << 368
384 G4ThreeVector finalP = dp->GetMomentum() - d << 369 vector<G4DynamicParticle*>* G4MollerBhabhaModel::SampleSecondaries(
385 finalP = finalP.unit(); << 370 const G4MaterialCutsCouple* couple,
>> 371 const G4DynamicParticle* dp,
>> 372 G4double tmin,
>> 373 G4double maxEnergy)
>> 374 {
>> 375 vector<G4DynamicParticle*>* vdp = new vector<G4DynamicParticle*>;
>> 376 G4DynamicParticle* delta = SampleSecondary(couple, dp, tmin, maxEnergy);
>> 377 vdp->push_back(delta);
386 378
387 fParticleChange->SetProposedKineticEnergy(ki << 379 return vdp;
388 fParticleChange->SetProposedMomentumDirectio <<
389 } 380 }
390 381
391 //....oooOO0OOooo........oooOO0OOooo........oo << 382 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
392 383