203{
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
220
221 if (verboseLevel > 3) {
222 G4cout <<
"G4LowEPComptonModel::SampleSecondaries() E(MeV)= "
225 }
226
227
228 if (photonEnergy0 <= lowEnergyLimit)
229 {
233 return ;
234 }
235
236 G4double e0m = photonEnergy0 / electron_mass_c2 ;
238
239
243
244 G4double LowEPCepsilon0 = 1. / (1. + 2. * e0m);
245 G4double LowEPCepsilon0Sq = LowEPCepsilon0 * LowEPCepsilon0;
246 G4double alpha1 = -std::log(LowEPCepsilon0);
247 G4double alpha2 = 0.5 * (1. - LowEPCepsilon0Sq);
248
249 G4double wlPhoton = h_Planck*c_light/photonEnergy0;
250
251
257
258 do
259 {
261 {
263 LowEPCepsilonSq = LowEPCepsilon * LowEPCepsilon;
264 }
265 else
266 {
267 LowEPCepsilonSq = LowEPCepsilon0Sq + (1. - LowEPCepsilon0Sq) *
G4UniformRand();
268 LowEPCepsilon = std::sqrt(LowEPCepsilonSq);
269 }
270
271 oneCosT = (1. - LowEPCepsilon) / ( LowEPCepsilon * e0m);
272 sinT2 = oneCosT * (2. - oneCosT);
273 G4double x = std::sqrt(oneCosT/2.) / (wlPhoton/cm);
275 gReject = (1. - LowEPCepsilon * sinT2 / (1. + LowEPCepsilonSq)) * scatteringFunction;
276
278
280 G4double sinTheta = std::sqrt(sinT2);
282 G4double dirx = sinTheta * std::cos(phi);
283 G4double diry = sinTheta * std::sin(phi);
285
286
287
288
289
290
291
292
293
295 G4double momentum_au_to_nat = (
pi/2.0)*1.992851740*std::pow(10.,-24.);
296 G4double e_mass_kg = 9.10938188 * std::pow(10.,-31.);
297
298 G4int maxDopplerIterations = 1000;
300 G4double pEIncident = photonEnergy0 ;
305
308
315 do{
316
317
318
319
320
321
322 do
323 {
324 iteration++;
325
326
327
328
329
330
331
332
335
336
337
338
340
341
342
343 G4double ePSI = ePAU * momentum_au_to_nat;
344
345
346
347 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / ((e_mass_kg*e_mass_kg)*(vel_c*vel_c)+(ePSI*ePSI)))/vel_c;
348
349
350
353
354
355
356 G4double eEIncident = electron_mass_c2 / sqrt( 1 - (u_temp*u_temp));
357 G4double systemE = eEIncident + pEIncident;
358
359
360 G4double gamma_temp = 1.0 / sqrt( 1 - (u_temp*u_temp));
361 G4double numerator = gamma_temp*electron_mass_c2*(1 - u_temp * std::cos(e_alpha));
362 G4double subdenom1 = u_temp*cosTheta*std::cos(e_alpha);
363 G4double subdenom2 = u_temp*sinTheta*std::sin(e_alpha)*std::cos(e_beta);
364 G4double denominator = (1.0 - cosTheta) + (gamma_temp*electron_mass_c2*(1 - subdenom1 - subdenom2) / pEIncident);
365 pERecoil = (numerator/denominator);
366 eERecoil = systemE - pERecoil;
367 CE_emission_flag = pEIncident - pERecoil;
368 } while ( (iteration <= maxDopplerIterations) && (CE_emission_flag < bindingE));
369
370
371
372
373
374
375
376
377
378
379
380
381
382 G4double a_temp = eERecoil / electron_mass_c2;
383 G4double u_p_temp = sqrt(1 - (1 / (a_temp*a_temp)));
384
385
386
387 G4double sinAlpha = std::sin(e_alpha);
388 G4double cosAlpha = std::cos(e_alpha);
389 G4double sinBeta = std::sin(e_beta);
390 G4double cosBeta = std::cos(e_beta);
391
392 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_temp));
393 G4double gamma_p = 1.0 / sqrt(1 - (u_p_temp*u_p_temp));
394
395 G4double var_A = pERecoil*u_p_temp*sinTheta;
396 G4double var_B = u_p_temp* (pERecoil*cosTheta-pEIncident);
397 G4double var_C = (pERecoil-pEIncident) - ( (pERecoil*pEIncident) / (gamma_p*electron_mass_c2))*(1 - cosTheta);
398
399 G4double var_D1 = gamma*electron_mass_c2*pERecoil;
400 G4double var_D2 = (1 - (u_temp*cosTheta*cosAlpha) - (u_temp*sinTheta*cosBeta*sinAlpha));
401 G4double var_D3 = ((electron_mass_c2*electron_mass_c2)*(gamma*gamma_p - 1)) - (gamma_p*electron_mass_c2*pERecoil);
402 G4double var_D = var_D1*var_D2 + var_D3;
403
404 G4double var_E1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosAlpha);
405 G4double var_E2 = gamma_p*electron_mass_c2*pERecoil*u_p_temp*cosTheta;
407
408 G4double var_F1 = ((gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*cosBeta*sinAlpha);
409 G4double var_F2 = (gamma_p*electron_mass_c2*pERecoil*u_p_temp*sinTheta);
411
412 G4double var_G = (gamma*gamma_p)*(electron_mass_c2*electron_mass_c2)*(u_temp*u_p_temp)*sinBeta*sinAlpha;
413
414
415
416
417 G4double var_W1 = (var_F*var_B - var_E*var_A)*(var_F*var_B - var_E*var_A);
418 G4double var_W2 = (var_G*var_G)*(var_A*var_A) + (var_G*var_G)*(var_B*var_B);
420
421 G4double var_Y = 2.0*(((var_A*var_D-var_F*var_C)*(var_F*var_B-var_E*var_A)) - ((var_G*var_G)*var_B*var_C));
422
423 G4double var_Z1 = (var_A*var_D - var_F*var_C)*(var_A*var_D - var_F*var_C);
424 G4double var_Z2 = (var_G*var_G)*(var_C*var_C) - (var_G*var_G)*(var_A*var_A);
426
427 G4double diff = ((var_Y*var_Y)-4*var_W*var_Z);
428
429
430
431 if (diff < 0.0)
432 {
433
436
438 funorder = abs(
G4lrint(std::log10(
static_cast<double>(abs(diff1)))+std::numeric_limits<double>::epsilon()))+sf;
439 diff1 =
G4lrint(diff1*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
440
442 funorder = abs(
G4lrint(std::log10(
static_cast<double>(abs(diff2)))+std::numeric_limits<double>::epsilon()))+sf;
443 diff2 =
G4lrint(diff2*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
444
445 diff = diff1 -diff2;
446
447 }
448
449
450
451
452 G4double X_p = (-var_Y + sqrt (diff))/(2*var_W);
453 G4double X_m = (-var_Y - sqrt (diff))/(2*var_W);
454
455
456
459
460 if (sol_select < 0.5)
461 {
462 ThetaE = std::acos(X_p);
463 }
464 if (sol_select > 0.5)
465 {
466 ThetaE = std::acos(X_m);
467 }
468
469 cosThetaE = std::cos(ThetaE);
470 sinThetaE = std::sin(ThetaE);
471 G4double Theta = std::acos(cosTheta);
472
473
474 G4double iSinThetaE = std::sqrt(1+std::tan((pi/2.0)-ThetaE)*std::tan((pi/2.0)-ThetaE));
475 G4double iSinTheta = std::sqrt(1+std::tan((pi/2.0)-Theta)*std::tan((pi/2.0)-Theta));
476 G4double ivar_A = iSinTheta/ (pERecoil*u_p_temp);
477
478 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_A*iSinThetaE);
479
480
481
482 if (abs(cosPhiE) > 1.0)
483 {
486 funorder = abs(
G4lrint(std::log10(
static_cast<double>(abs(cosPhiE)))+std::numeric_limits<double>::epsilon()))+sf;
487 cosPhiE =
G4lrint(cosPhiE*std::pow(10.0,funorder))/(1.0*std::pow(10.0,funorder));
488
489 }
490
491
492
493
494 } while ( (iteration <= maxDopplerIterations) && (abs(cosPhiE) > 1));
495
496
497
498
499 if (iteration >= maxDopplerIterations)
500 {
501 pERecoil = photonEnergy0 ;
502 bindingE = 0.;
503 dirx=0.0;
504 diry=0.0;
505 dirz=1.0;
506 }
507
508
509
511 photonDirection1.rotateUz(photonDirection0);
513
514 if (pERecoil > 0.)
515 {
517
518
520 G4double eDirX = sinThetaE * std::cos(phi+PhiE);
521 G4double eDirY = sinThetaE * std::sin(phi+PhiE);
523
524 G4double eKineticEnergy = pEIncident - pERecoil - bindingE;
525
527 eDirection.rotateUz(photonDirection0);
529 eDirection,eKineticEnergy) ;
530 fvect->push_back(dp);
531
532 }
533 else
534 {
537 }
538
539
540
541
542
543
544 if(fAtomDeexcitation && iteration < maxDopplerIterations) {
547 size_t nbefore = fvect->size();
551 size_t nafter = fvect->size();
552 if(nafter > nbefore) {
553 for (size_t i=nbefore; i<nafter; ++i) {
554 bindingE -= ((*fvect)[i])->GetKineticEnergy();
555 }
556 }
557 }
558 }
559 if(bindingE < 0.0) { bindingE = 0.0; }
561
562}
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4String & GetName() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double BindingEnergy(G4int Z, G4int shellIndex) const
G4int SelectRandomShell(G4int Z) const
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)