172{
173
174
175
176
177
178
179 if (verboseLevel > 3)
180 G4cout <<
"Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" <<
G4endl;
181
183
184
185 if(photonEnergy <= lowEnergyLimit)
186 {
189 return;
190 }
191
192
195
196
197
198
199 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.
mag()==0))
200 {
201 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
202 }
203 else
204 {
206 {
207 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
208 }
209 }
210
211
212
213
215 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
216
217
218
219 if (photonEnergy < smallEnergy )
220 {
221 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
222 }
223 else
224 {
225
226
227
228
229
230
233
234
235
236
237
238
239
240
242
243
244
245
246
247
248
249
250
251
252
254 if (photonEnergy > 50. * MeV) fZ += 8. * (element->
GetfCoulomb());
255
256
258 G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
259 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
260
261
262 G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
263 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
264 G4double epsilonRange = 0.5 - epsilonMin ;
265
266
269
270 G4double f10 = ScreenFunction1(screenMin) - fZ;
271 G4double f20 = ScreenFunction2(screenMin) - fZ;
272 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
273 G4double normF2 = std::max(1.5 * f20,0.);
274
275 do {
277 {
278 epsilon = 0.5 - epsilonRange * pow(
G4UniformRand(), 0.3333) ;
279 screen = screenFactor / (epsilon * (1. - epsilon));
280 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
281 }
282 else
283 {
285 screen = screenFactor / (epsilon * (1 - epsilon));
286 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
287
288
289 }
291
292 }
293
294
295
298
299
301 {
302 electronTotEnergy = (1. - epsilon) * photonEnergy;
303 positronTotEnergy = epsilon * photonEnergy;
304 }
305 else
306 {
307 positronTotEnergy = (1. - epsilon) * photonEnergy;
308 electronTotEnergy = epsilon * photonEnergy;
309 }
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330 G4double Ene = electronTotEnergy/electron_mass_c2;
331
334
335 SetTheta(&cosTheta,&sinTheta,Ene);
336
337
338
339
341
342
343
344
345
346 phi = SetPhi(photonEnergy);
347 psi = SetPsi(photonEnergy,phi);
348
349
350
351
352
353
354
355 Psi = psi;
356 Phi = phi;
357
358
359
361
366
367
368
369
370
371
372 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
373
374 SystemOfRefChange(gammaDirection0,electronDirection,
375 gammaPolarization0);
376
378 electronDirection,
379 electronKineEnergy);
380
381
382
383 Ene = positronTotEnergy/electron_mass_c2;
384
385 cosTheta = 0.;
386 sinTheta = 0.;
387
388 SetTheta(&cosTheta,&sinTheta,Ene);
390
391 dirX = sinTheta*cos(phip);
392 dirY = sinTheta*sin(phip);
393 dirZ = cosTheta;
395
396 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
397 SystemOfRefChange(gammaDirection0,positronDirection,
398 gammaPolarization0);
399
400
402 positronDirection, positronKineEnergy);
403
404
405 fvect->push_back(particle1);
406 fvect->push_back(particle2);
407
408
409
410
411
412
413
414
415
416
417
418
422
423}
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
double howOrthogonal(const Hep3Vector &v) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4ThreeVector & GetPolarization() const
static G4Electron * Electron()
G4double GetfCoulomb() const
G4IonisParamElm * GetIonisation() const
G4double GetlogZ3() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
static G4Positron * Positron()
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)