235{
236
237
238
239
240
241
242
243
244
245
246
247 if (verboseLevel > 1)
248 G4cout <<
"Calling SampleSecondaries() of G4LivermoreGammaConversionModel" <<
G4endl;
249
252
254 G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
255
256
257 if (photonEnergy < smallEnergy )
258 {
259 epsilon = epsilon0Local + (0.5 - epsilon0Local) *
G4UniformRand();
260 }
261 else
262 {
263
264
267
268 if (element == 0)
269 {
270 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - element = 0"
272 return;
273 }
275 if (ionisation == 0)
276 {
277 G4cout <<
"G4LivermoreGammaConversionModel::SampleSecondaries - ionisation = 0"
279 return;
280 }
281
282
284 if (photonEnergy > 50. * MeV) fZ += 8. * (element->
GetfCoulomb());
285
286
288 G4double screenMax = std::exp ((42.24 - fZ)/8.368) - 0.952 ;
289 G4double screenMin = std::min(4.*screenFactor,screenMax) ;
290
291
292 G4double epsilon1 = 0.5 - 0.5 * std::sqrt(1. - screenMin / screenMax) ;
293 G4double epsilonMin = std::max(epsilon0Local,epsilon1);
294 G4double epsilonRange = 0.5 - epsilonMin ;
295
296
299
300 G4double f10 = ScreenFunction1(screenMin) - fZ;
301 G4double f20 = ScreenFunction2(screenMin) - fZ;
302 G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
303 G4double normF2 = std::max(1.5 * f20,0.);
304
305 do
306 {
308 {
309 epsilon = 0.5 - epsilonRange * std::pow(
G4UniformRand(), 0.333333) ;
310 screen = screenFactor / (epsilon * (1. - epsilon));
311 gReject = (ScreenFunction1(screen) - fZ) / f10 ;
312 }
313 else
314 {
316 screen = screenFactor / (epsilon * (1 - epsilon));
317 gReject = (ScreenFunction2(screen) - fZ) / f20 ;
318 }
320
321 }
322
323
324
327
329 {
330 electronTotEnergy = (1. - epsilon) * photonEnergy;
331 positronTotEnergy = epsilon * photonEnergy;
332 }
333 else
334 {
335 positronTotEnergy = (1. - epsilon) * photonEnergy;
336 electronTotEnergy = epsilon * photonEnergy;
337 }
338
339
340
341
342
346
347
348
350 {
352 }
353 else
354 {
356 }
357
358 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
359 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
361
362 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
363 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
364
365
366
367
368
369
370 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
371
373 electronDirection.rotateUz(photonDirection);
374
376 electronDirection,
377 electronKineEnergy);
378
379
380 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
381
383 positronDirection.rotateUz(photonDirection);
384
385
387 positronDirection,
388 positronKineEnergy);
389
390 fvect->push_back(particle1);
391 fvect->push_back(particle2);
392
393
396
397}
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
G4double GetfCoulomb() const
G4IonisParamElm * GetIonisation() const
G4double GetlogZ3() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4Positron * Positron()
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)