Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NuclNuclDiffuseElastic.hh
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id$
28//
29//
30// G4 Model: optical elastic scattering with 4-momentum balance
31//
32// Class Description
33// Final state production model for nucleus-nucleus elastic scattering;
34// Coulomb amplitude is not considered as correction
35// (as in G4DiffuseElastic)
36// Class Description - End
37//
38//
39// 17.03.09 V. Grichine implementation for Coulomb elastic scattering
40
41
42#ifndef G4NuclNuclDiffuseElastic_h
43#define G4NuclNuclDiffuseElastic_h 1
44
45#include <complex>
46#include <CLHEP/Units/PhysicalConstants.h>
47#include "globals.hh"
48#include "G4Integrator.hh"
49#include "G4HadronElastic.hh"
50#include "G4HadProjectile.hh"
51#include "G4Nucleus.hh"
52
53using namespace std;
54
55class G4ParticleDefinition;
56class G4PhysicsTable;
57class G4PhysicsLogVector;
58
59class G4NuclNuclDiffuseElastic : public G4HadronElastic // G4HadronicInteraction
60{
61public:
62
63 G4NuclNuclDiffuseElastic();
64
65 // G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
66
67
68
69
70
71 virtual ~G4NuclNuclDiffuseElastic();
72
73 void Initialise();
74
75 void InitialiseOnFly(G4double Z, G4double A);
76
77 void BuildAngleTable();
78
79
80 // G4HadFinalState * ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
81
82 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
83 G4double plab,
84 G4int Z, G4int A);
85
86 void SetPlabLowLimit(G4double value);
87
88 void SetHEModelLowLimit(G4double value);
89
90 void SetQModelLowLimit(G4double value);
91
92 void SetLowestEnergyLimit(G4double value);
93
94 void SetRecoilKinEnergyLimit(G4double value);
95
96 G4double SampleT(const G4ParticleDefinition* aParticle,
97 G4double p, G4double A);
98
99 G4double SampleTableT(const G4ParticleDefinition* aParticle,
100 G4double p, G4double Z, G4double A);
101
102 G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
103
104 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
105 G4double Z, G4double A);
106
107 G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
108
109 G4double SampleThetaLab(const G4HadProjectile* aParticle,
110 G4double tmass, G4double A);
111
112 G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
113 G4double theta,
114 G4double momentum,
115 G4double A );
116
117 G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
118 G4double theta,
119 G4double momentum,
120 G4double A, G4double Z );
121
122 G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
123 G4double theta,
124 G4double momentum,
125 G4double A, G4double Z );
126
127 G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
128 G4double tMand,
129 G4double momentum,
130 G4double A, G4double Z );
131
132 G4double IntegralElasticProb( const G4ParticleDefinition* particle,
133 G4double theta,
134 G4double momentum,
135 G4double A );
136
137
138 G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
139 G4double theta,
140 G4double momentum,
141 G4double Z );
142
143 G4double GetRutherfordXsc( G4double theta );
144
145 G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
146 G4double tMand,
147 G4double momentum,
148 G4double A, G4double Z );
149
150 G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
151 G4double momentum, G4double Z );
152
153 G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
154 G4double momentum, G4double Z,
155 G4double theta1, G4double theta2 );
156
157
158 G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
159 G4double momentum );
160
161 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
162
163 G4double CalculateAm( G4double momentum, G4double n, G4double Z);
164
165 G4double CalculateNuclearRad( G4double A);
166
167 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
168 G4double tmass, G4double thetaCMS);
169
170 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
171 G4double tmass, G4double thetaLab);
172
173 void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
174 G4double Z, G4double A);
175
176
177
178 G4double BesselJzero(G4double z);
179 G4double BesselJone(G4double z);
180 G4double DampFactor(G4double z);
181 G4double BesselOneByArg(G4double z);
182
183 G4double GetDiffElasticProb(G4double theta);
184 G4double GetDiffElasticSumProb(G4double theta);
185 G4double GetDiffElasticSumProbA(G4double alpha);
186 G4double GetIntegrandFunction(G4double theta);
187
188 G4double GetNuclearRadius(){return fNuclearRadius;};
189
190
191 // Technical math functions for strong Coulomb contribution
192
193 G4complex GammaLogarithm(G4complex xx);
194 G4complex GammaLogB2n(G4complex xx);
195
196 G4double GetErf(G4double x);
197
198 G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
199 G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
200
201 G4double GetCint(G4double x);
202 G4double GetSint(G4double x);
203
204
205 G4complex GetErfcComp(G4complex z, G4int nMax);
206 G4complex GetErfcSer(G4complex z, G4int nMax);
207 G4complex GetErfcInt(G4complex z); // , G4int nMax);
208
209 G4complex GetErfComp(G4complex z, G4int nMax); // AandS algorithm != Ser, Int
210 G4complex GetErfSer(G4complex z, G4int nMax);
211
212 G4double GetExpCos(G4double x);
213 G4double GetExpSin(G4double x);
214 G4complex GetErfInt(G4complex z); // , G4int nMax);
215
216
217
218
219 G4double GetLegendrePol(G4int n, G4double x);
220
221 G4complex TestErfcComp(G4complex z, G4int nMax);
222 G4complex TestErfcSer(G4complex z, G4int nMax);
223 G4complex TestErfcInt(G4complex z); // , G4int nMax);
224
225 G4complex CoulombAmplitude(G4double theta);
226 G4double CoulombAmplitudeMod2(G4double theta);
227
228
229 void CalculateCoulombPhaseZero();
230 G4double CalculateCoulombPhase(G4int n);
231 void CalculateRutherfordAnglePar();
232
233 G4double ProfileNear(G4double theta);
234 G4double ProfileFar(G4double theta);
235 G4double Profile(G4double theta);
236
237 G4complex PhaseNear(G4double theta);
238 G4complex PhaseFar(G4double theta);
239
240 G4complex GammaLess(G4double theta);
241 G4complex GammaMore(G4double theta);
242
243 G4complex AmplitudeNear(G4double theta);
244 G4complex AmplitudeFar(G4double theta);
245
246 G4complex Amplitude(G4double theta);
247 G4double AmplitudeMod2(G4double theta);
248
249 G4complex AmplitudeSim(G4double theta);
250 G4double AmplitudeSimMod2(G4double theta);
251
252 G4double GetRatioSim(G4double theta);
253 G4double GetRatioGen(G4double theta);
254
255 G4double GetFresnelDiffuseXsc(G4double theta);
256 G4double GetFresnelIntegrandXsc(G4double alpha);
257
258
259 G4complex AmplitudeGla(G4double theta);
260 G4double AmplitudeGlaMod2(G4double theta);
261
262 G4complex AmplitudeGG(G4double theta);
263 G4double AmplitudeGGMod2(G4double theta);
264
265 void InitParameters(const G4ParticleDefinition* theParticle,
266 G4double partMom, G4double Z, G4double A);
267
268 void InitDynParameters(const G4ParticleDefinition* theParticle,
269 G4double partMom);
270
271 void InitParametersGla(const G4DynamicParticle* aParticle,
272 G4double partMom, G4double Z, G4double A);
273
274 G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
275 G4double pTkin,
276 G4ParticleDefinition* tParticle);
277
278 G4double CalcMandelstamS( const G4double mp ,
279 const G4double mt ,
280 const G4double Plab );
281
282 G4double GetProfileLambda(){return fProfileLambda;};
283
284 void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
285 void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
286 void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
287 void SetCofLambda(G4double pa){fCofLambda = pa;};
288
289 void SetCofAlpha(G4double pa){fCofAlpha = pa;};
290 void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
291 void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
292
293 void SetCofDelta(G4double pa){fCofDelta = pa;};
294 void SetCofPhase(G4double pa){fCofPhase = pa;};
295 void SetCofFar(G4double pa){fCofFar = pa;};
296 void SetEtaRatio(G4double pa){fEtaRatio = pa;};
297 void SetMaxL(G4int l){fMaxL = l;};
298 void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
299
300 G4double GetCofAlphaMax(){return fCofAlphaMax;};
301 G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
302
303private:
304
305
306 G4ParticleDefinition* theProton;
307 G4ParticleDefinition* theNeutron;
308 G4ParticleDefinition* theDeuteron;
309 G4ParticleDefinition* theAlpha;
310
311 const G4ParticleDefinition* thePionPlus;
312 const G4ParticleDefinition* thePionMinus;
313
314 G4double lowEnergyRecoilLimit;
315 G4double lowEnergyLimitHE;
316 G4double lowEnergyLimitQ;
317 G4double lowestEnergyLimit;
318 G4double plabLowLimit;
319
320 G4int fEnergyBin;
321 G4int fAngleBin;
322
323 G4PhysicsLogVector* fEnergyVector;
324 G4PhysicsTable* fAngleTable;
325 std::vector<G4PhysicsTable*> fAngleBank;
326
327 std::vector<G4double> fElementNumberVector;
328 std::vector<G4String> fElementNameVector;
329
330 const G4ParticleDefinition* fParticle;
331
332 G4double fWaveVector;
333 G4double fAtomicWeight;
334 G4double fAtomicNumber;
335
336 G4double fNuclearRadius1;
337 G4double fNuclearRadius2;
338 G4double fNuclearRadius;
339 G4double fNuclearRadiusSquare;
340 G4double fNuclearRadiusCof;
341
342 G4double fBeta;
343 G4double fZommerfeld;
344 G4double fRutherfordRatio;
345 G4double fAm;
346 G4bool fAddCoulomb;
347
348 G4double fCoulombPhase0;
349 G4double fHalfRutThetaTg;
350 G4double fHalfRutThetaTg2;
351 G4double fRutherfordTheta;
352
353 G4double fProfileLambda;
354 G4double fProfileDelta;
355 G4double fProfileAlpha;
356
357 G4double fCofLambda;
358 G4double fCofAlpha;
359 G4double fCofDelta;
360 G4double fCofPhase;
361 G4double fCofFar;
362
363 G4double fCofAlphaMax;
364 G4double fCofAlphaCoulomb;
365
366 G4int fMaxL;
367 G4double fSumSigma;
368 G4double fEtaRatio;
369
370 G4double fReZ;
371
372};
373
374
375inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
376{
377 lowEnergyRecoilLimit = value;
378}
379
380inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
381{
382 plabLowLimit = value;
383}
384
385inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
386{
387 lowEnergyLimitHE = value;
388}
389
390inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
391{
392 lowEnergyLimitQ = value;
393}
394
395inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
396{
397 lowestEnergyLimit = value;
398}
399
400
401/////////////////////////////////////////////////////////////
402//
403// Bessel J0 function based on rational approximation from
404// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
405
406inline G4double G4NuclNuclDiffuseElastic::BesselJzero(G4double value)
407{
408 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
409
410 modvalue = fabs(value);
411
412 if ( value < 8.0 && value > -8.0 )
413 {
414 value2 = value*value;
415
416 fact1 = 57568490574.0 + value2*(-13362590354.0
417 + value2*( 651619640.7
418 + value2*(-11214424.18
419 + value2*( 77392.33017
420 + value2*(-184.9052456 ) ) ) ) );
421
422 fact2 = 57568490411.0 + value2*( 1029532985.0
423 + value2*( 9494680.718
424 + value2*(59272.64853
425 + value2*(267.8532712
426 + value2*1.0 ) ) ) );
427
428 bessel = fact1/fact2;
429 }
430 else
431 {
432 arg = 8.0/modvalue;
433
434 value2 = arg*arg;
435
436 shift = modvalue-0.785398164;
437
438 fact1 = 1.0 + value2*(-0.1098628627e-2
439 + value2*(0.2734510407e-4
440 + value2*(-0.2073370639e-5
441 + value2*0.2093887211e-6 ) ) );
442
443 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
444 + value2*(-0.6911147651e-5
445 + value2*(0.7621095161e-6
446 - value2*0.934945152e-7 ) ) );
447
448 bessel = sqrt(0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2 );
449 }
450 return bessel;
451}
452
453/////////////////////////////////////////////////////////////
454//
455// Bessel J1 function based on rational approximation from
456// J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
457
458inline G4double G4NuclNuclDiffuseElastic::BesselJone(G4double value)
459{
460 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
461
462 modvalue = fabs(value);
463
464 if ( modvalue < 8.0 )
465 {
466 value2 = value*value;
467
468 fact1 = value*(72362614232.0 + value2*(-7895059235.0
469 + value2*( 242396853.1
470 + value2*(-2972611.439
471 + value2*( 15704.48260
472 + value2*(-30.16036606 ) ) ) ) ) );
473
474 fact2 = 144725228442.0 + value2*(2300535178.0
475 + value2*(18583304.74
476 + value2*(99447.43394
477 + value2*(376.9991397
478 + value2*1.0 ) ) ) );
479 bessel = fact1/fact2;
480 }
481 else
482 {
483 arg = 8.0/modvalue;
484
485 value2 = arg*arg;
486
487 shift = modvalue - 2.356194491;
488
489 fact1 = 1.0 + value2*( 0.183105e-2
490 + value2*(-0.3516396496e-4
491 + value2*(0.2457520174e-5
492 + value2*(-0.240337019e-6 ) ) ) );
493
494 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
495 + value2*( 0.8449199096e-5
496 + value2*(-0.88228987e-6
497 + value2*0.105787412e-6 ) ) );
498
499 bessel = sqrt( 0.636619772/modvalue)*(cos(shift)*fact1 - arg*sin(shift)*fact2);
500
501 if (value < 0.0) bessel = -bessel;
502 }
503 return bessel;
504}
505
506////////////////////////////////////////////////////////////////////
507//
508// damp factor in diffraction x/sh(x), x was already *pi
509
510inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
511{
512 G4double df;
513 G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
514
515 // x *= pi;
516
517 if( std::fabs(x) < 0.01 )
518 {
519 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
520 }
521 else
522 {
523 df = x/std::sinh(x);
524 }
525 return df;
526}
527
528
529////////////////////////////////////////////////////////////////////
530//
531// return J1(x)/x with special case for small x
532
533inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
534{
535 G4double x2, result;
536
537 if( std::fabs(x) < 0.01 )
538 {
539 x *= 0.5;
540 x2 = x*x;
541 result = 2. - x2 + x2*x2/6.;
542 }
543 else
544 {
545 result = BesselJone(x)/x;
546 }
547 return result;
548}
549
550////////////////////////////////////////////////////////////////////
551//
552// return particle beta
553
554inline G4double G4NuclNuclDiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
555 G4double momentum )
556{
557 G4double mass = particle->GetPDGMass();
558 G4double a = momentum/mass;
559 fBeta = a/std::sqrt(1+a*a);
560
561 return fBeta;
562}
563
564////////////////////////////////////////////////////////////////////
565//
566// return Zommerfeld parameter for Coulomb scattering
567
568inline G4double G4NuclNuclDiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
569{
570 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
571
572 return fZommerfeld;
573}
574
575////////////////////////////////////////////////////////////////////
576//
577// return Wentzel correction for Coulomb scattering
578
579inline G4double G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
580{
581 G4double k = momentum/CLHEP::hbarc;
582 G4double ch = 1.13 + 3.76*n*n;
583 G4double zn = 1.77*k*std::pow(Z,-1./3.)*CLHEP::Bohr_radius;
584 G4double zn2 = zn*zn;
585 fAm = ch/zn2;
586
587 return fAm;
588}
589
590////////////////////////////////////////////////////////////////////
591//
592// calculate nuclear radius for different atomic weights using different approximations
593
594inline G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
595{
596 G4double r0 = 1.*CLHEP::fermi, radius;
597 // r0 *= 1.12;
598 // r0 *= 1.44;
599 r0 *= fNuclearRadiusCof;
600
601 /*
602 if( A < 50. )
603 {
604 if( A > 10. ) r0 = 1.16*( 1 - std::pow(A, -2./3.) )*CLHEP::fermi; // 1.08*fermi;
605 else r0 = 1.1*CLHEP::fermi;
606
607 radius = r0*std::pow(A, 1./3.);
608 }
609 else
610 {
611 r0 = 1.7*CLHEP::fermi; // 1.7*fermi;
612
613 radius = r0*std::pow(A, 0.27); // 0.27);
614 }
615 */
616 radius = r0*std::pow(A, 1./3.);
617
618 return radius;
619}
620
621////////////////////////////////////////////////////////////////////
622//
623// return Coulomb scattering differential xsc with Wentzel correction. Test function
624
625inline G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
626 G4double theta,
627 G4double momentum,
628 G4double Z )
629{
630 G4double sinHalfTheta = std::sin(0.5*theta);
631 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
632 G4double beta = CalculateParticleBeta( particle, momentum);
633 G4double z = particle->GetPDGCharge();
634 G4double n = CalculateZommerfeld( beta, z, Z );
635 G4double am = CalculateAm( momentum, n, Z);
636 G4double k = momentum/CLHEP::hbarc;
637 G4double ch = 0.5*n/k;
638 G4double ch2 = ch*ch;
639 G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
640
641 return xsc;
642}
643
644////////////////////////////////////////////////////////////////////
645//
646// return Rutherford scattering differential xsc with Wentzel correction. For Sampling.
647
648inline G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc( G4double theta )
649{
650 G4double sinHalfTheta = std::sin(0.5*theta);
651 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
652
653 G4double ch2 = fRutherfordRatio*fRutherfordRatio;
654
655 G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
656
657 return xsc;
658}
659
660
661////////////////////////////////////////////////////////////////////
662//
663// return Coulomb scattering total xsc with Wentzel correction
664
665inline G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
666 G4double momentum, G4double Z )
667{
668 G4double beta = CalculateParticleBeta( particle, momentum);
669 G4cout<<"beta = "<<beta<<G4endl;
670 G4double z = particle->GetPDGCharge();
671 G4double n = CalculateZommerfeld( beta, z, Z );
672 G4cout<<"fZomerfeld = "<<n<<G4endl;
673 G4double am = CalculateAm( momentum, n, Z);
674 G4cout<<"cof Am = "<<am<<G4endl;
675 G4double k = momentum/CLHEP::hbarc;
676 G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
677 G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
678 G4double ch = n/k;
679 G4double ch2 = ch*ch;
680 G4double xsc = ch2*CLHEP::pi/(am +am*am);
681
682 return xsc;
683}
684
685////////////////////////////////////////////////////////////////////
686//
687// return Coulomb scattering xsc with Wentzel correction integrated between
688// theta1 and < theta2
689
690inline G4double G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
691 G4double momentum, G4double Z,
692 G4double theta1, G4double theta2 )
693{
694 G4double c1 = std::cos(theta1);
695 G4cout<<"c1 = "<<c1<<G4endl;
696 G4double c2 = std::cos(theta2);
697 G4cout<<"c2 = "<<c2<<G4endl;
698 G4double beta = CalculateParticleBeta( particle, momentum);
699 // G4cout<<"beta = "<<beta<<G4endl;
700 G4double z = particle->GetPDGCharge();
701 G4double n = CalculateZommerfeld( beta, z, Z );
702 // G4cout<<"fZomerfeld = "<<n<<G4endl;
703 G4double am = CalculateAm( momentum, n, Z);
704 // G4cout<<"cof Am = "<<am<<G4endl;
705 G4double k = momentum/CLHEP::hbarc;
706 // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
707 // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
708 G4double ch = n/k;
709 G4double ch2 = ch*ch;
710 am *= 2.;
711 G4double xsc = ch2*CLHEP::twopi*(c1-c2);
712 xsc /= (1 - c1 + am)*(1 - c2 + am);
713
714 return xsc;
715}
716
717///////////////////////////////////////////////////////////////////
718//
719// For the calculation of arg Gamma(z) one needs complex extension
720// of ln(Gamma(z))
721
722inline G4complex G4NuclNuclDiffuseElastic::GammaLogarithm(G4complex zz)
723{
724 static G4double cof[6] = { 76.18009172947146, -86.50532032941677,
725 24.01409824083091, -1.231739572450155,
726 0.1208650973866179e-2, -0.5395239384953e-5 } ;
727 register G4int j;
728 G4complex z = zz - 1.0;
729 G4complex tmp = z + 5.5;
730 tmp -= (z + 0.5) * std::log(tmp);
731 G4complex ser = G4complex(1.000000000190015,0.);
732
733 for ( j = 0; j <= 5; j++ )
734 {
735 z += 1.0;
736 ser += cof[j]/z;
737 }
738 return -tmp + std::log(2.5066282746310005*ser);
739}
740
741///////////////////////////////////////////////////////////////////
742//
743// For the calculation of arg Gamma(z) one needs complex extension
744// of ln(Gamma(z)) here is approximate algorithm
745
746inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
747{
748 G4complex z1 = 12.*z;
749 G4complex z2 = z*z;
750 G4complex z3 = z2*z;
751 G4complex z5 = z2*z3;
752 G4complex z7 = z2*z5;
753
754 z3 *= 360.;
755 z5 *= 1260.;
756 z7 *= 1680.;
757
758 G4complex result = (z-0.5)*std::log(z)-z+0.5*std::log(CLHEP::twopi);
759 result += 1./z1 - 1./z3 +1./z5 -1./z7;
760 return result;
761}
762
763/////////////////////////////////////////////////////////////////
764//
765//
766
767inline G4double G4NuclNuclDiffuseElastic::GetErf(G4double x)
768{
769 G4double t, z, tmp, result;
770
771 z = std::fabs(x);
772 t = 1.0/(1.0+0.5*z);
773
774 tmp = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
775 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
776 t*(-0.82215223+t*0.17087277)))))))));
777
778 if( x >= 0.) result = 1. - tmp;
779 else result = 1. + tmp;
780
781 return result;
782}
783
784/////////////////////////////////////////////////////////////////
785//
786//
787
788inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
789{
790 G4complex erfcz = 1. - GetErfComp( z, nMax);
791 return erfcz;
792}
793
794/////////////////////////////////////////////////////////////////
795//
796//
797
798inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
799{
800 G4complex erfcz = 1. - GetErfSer( z, nMax);
801 return erfcz;
802}
803
804/////////////////////////////////////////////////////////////////
805//
806//
807
808inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
809{
810 G4complex erfcz = 1. - GetErfInt( z); // , nMax);
811 return erfcz;
812}
813
814inline G4double G4NuclNuclDiffuseElastic::GetLegendrePol(G4int n, G4double theta)
815{
816 G4double legPol, epsilon = 1.e-6;
817 G4double x = std::cos(theta);
818
819 if ( n < 0 ) legPol = 0.;
820 else if( n == 0 ) legPol = 1.;
821 else if( n == 1 ) legPol = x;
822 else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
823 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
824 else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
825 else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
826 else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
827 else
828 {
829 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
830
831 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
832 }
833 return legPol;
834}
835
836
837
838/////////////////////////////////////////////////////////////////
839//
840//
841
842inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
843{
844 G4complex miz = G4complex( z.imag(), -z.real() );
845 G4complex erfcz = 1. - GetErfComp( miz, nMax);
846 G4complex w = std::exp(-z*z)*erfcz;
847 return w;
848}
849
850/////////////////////////////////////////////////////////////////
851//
852//
853
854inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
855{
856 G4complex miz = G4complex( z.imag(), -z.real() );
857 G4complex erfcz = 1. - GetErfSer( miz, nMax);
858 G4complex w = std::exp(-z*z)*erfcz;
859 return w;
860}
861
862/////////////////////////////////////////////////////////////////
863//
864//
865
866inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
867{
868 G4complex miz = G4complex( z.imag(), -z.real() );
869 G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
870 G4complex w = std::exp(-z*z)*erfcz;
871 return w;
872}
873
874/////////////////////////////////////////////////////////////////
875//
876//
877
878inline G4complex G4NuclNuclDiffuseElastic::GetErfComp(G4complex z, G4int nMax)
879{
880 G4int n;
881 G4double n2, cofn, shny, chny, fn, gn;
882
883 G4double x = z.real();
884 G4double y = z.imag();
885
886 G4double outRe = 0., outIm = 0.;
887
888 G4double twox = 2.*x;
889 G4double twoxy = twox*y;
890 G4double twox2 = twox*twox;
891
892 G4double cof1 = std::exp(-x*x)/CLHEP::pi;
893
894 G4double cos2xy = std::cos(twoxy);
895 G4double sin2xy = std::sin(twoxy);
896
897 G4double twoxcos2xy = twox*cos2xy;
898 G4double twoxsin2xy = twox*sin2xy;
899
900 for( n = 1; n <= nMax; n++)
901 {
902 n2 = n*n;
903
904 cofn = std::exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
905
906 chny = std::cosh(n*y);
907 shny = std::sinh(n*y);
908
909 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
910 gn = twoxsin2xy*chny + n*cos2xy*shny;
911
912 fn *= cofn;
913 gn *= cofn;
914
915 outRe += fn;
916 outIm += gn;
917 }
918 outRe *= 2*cof1;
919 outIm *= 2*cof1;
920
921 if(std::abs(x) < 0.0001)
922 {
923 outRe += GetErf(x);
924 outIm += cof1*y;
925 }
926 else
927 {
928 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
929 outIm += cof1*sin2xy/twox;
930 }
931 return G4complex(outRe, outIm);
932}
933
934/////////////////////////////////////////////////////////////////
935//
936//
937
938inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
939{
940 G4int n;
941 G4double a =1., b = 1., tmp;
942 G4complex sum = z, d = z;
943
944 for( n = 1; n <= nMax; n++)
945 {
946 a *= 2.;
947 b *= 2.*n +1.;
948 d *= z*z;
949
950 tmp = a/b;
951
952 sum += tmp*d;
953 }
954 sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
955
956 return sum;
957}
958
959/////////////////////////////////////////////////////////////////////
960
961inline G4double G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
962{
963 G4double result;
964
965 result = std::exp(x*x-fReZ*fReZ);
966 result *= std::cos(2.*x*fReZ);
967 return result;
968}
969
970/////////////////////////////////////////////////////////////////////
971
972inline G4double G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
973{
974 G4double result;
975
976 result = std::exp(x*x-fReZ*fReZ);
977 result *= std::sin(2.*x*fReZ);
978 return result;
979}
980
981
982
983/////////////////////////////////////////////////////////////////
984//
985//
986
987inline G4complex G4NuclNuclDiffuseElastic::GetErfInt(G4complex z) // , G4int nMax)
988{
989 G4double outRe, outIm;
990
991 G4double x = z.real();
992 G4double y = z.imag();
993 fReZ = x;
994
995 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
996
997 outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
998 outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
999
1000 outRe *= 2./sqrt(CLHEP::pi);
1001 outIm *= 2./sqrt(CLHEP::pi);
1002
1003 outRe += GetErf(x);
1004
1005 return G4complex(outRe, outIm);
1006}
1007
1008/////////////////////////////////////////////////////////////////
1009//
1010//
1011
1012inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
1013{
1014 G4double out;
1015
1016 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1017
1018 out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
1019
1020 return out;
1021}
1022
1023/////////////////////////////////////////////////////////////////
1024//
1025//
1026
1027inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
1028{
1029 G4double out;
1030
1031 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1032
1033 out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
1034
1035 return out;
1036}
1037
1038
1039/////////////////////////////////////////////////////////////////
1040//
1041//
1042
1043inline G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
1044{
1045 G4complex ca;
1046
1047 G4double sinHalfTheta = std::sin(0.5*theta);
1048 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
1049 sinHalfTheta2 += fAm;
1050
1051 G4double order = 2.*fCoulombPhase0 - fZommerfeld*std::log(sinHalfTheta2);
1052 G4complex z = G4complex(0., order);
1053 ca = std::exp(z);
1054
1055 ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
1056
1057 return ca;
1058}
1059
1060/////////////////////////////////////////////////////////////////
1061//
1062//
1063
1064inline G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
1065{
1066 G4complex ca = CoulombAmplitude(theta);
1067 G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
1068
1069 return out;
1070}
1071
1072/////////////////////////////////////////////////////////////////
1073//
1074//
1075
1076
1077inline void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
1078{
1079 G4complex z = G4complex(1,fZommerfeld);
1080 // G4complex gammalog = GammaLogarithm(z);
1081 G4complex gammalog = GammaLogB2n(z);
1082 fCoulombPhase0 = gammalog.imag();
1083}
1084
1085/////////////////////////////////////////////////////////////////
1086//
1087//
1088
1089
1090inline G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
1091{
1092 G4complex z = G4complex(1. + n, fZommerfeld);
1093 // G4complex gammalog = GammaLogarithm(z);
1094 G4complex gammalog = GammaLogB2n(z);
1095 return gammalog.imag();
1096}
1097
1098
1099/////////////////////////////////////////////////////////////////
1100//
1101//
1102
1103
1104inline void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
1105{
1106 fHalfRutThetaTg = fZommerfeld/fProfileLambda; // (fWaveVector*fNuclearRadius);
1107 fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
1108 fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
1109 // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
1110
1111}
1112
1113/////////////////////////////////////////////////////////////////
1114//
1115//
1116
1117inline G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
1118{
1119 G4double dTheta = fRutherfordTheta - theta;
1120 G4double result = 0., argument = 0.;
1121
1122 if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
1123 else
1124 {
1125 argument = fProfileDelta*dTheta;
1126 result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1127 result /= std::sinh(CLHEP::pi*argument);
1128 result -= 1.;
1129 result /= dTheta;
1130 }
1131 return result;
1132}
1133
1134/////////////////////////////////////////////////////////////////
1135//
1136//
1137
1138inline G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
1139{
1140 G4double dTheta = fRutherfordTheta + theta;
1141 G4double argument = fProfileDelta*dTheta;
1142
1143 G4double result = CLHEP::pi*argument*std::exp(fProfileAlpha*argument);
1144 result /= std::sinh(CLHEP::pi*argument);
1145 result /= dTheta;
1146
1147 return result;
1148}
1149
1150/////////////////////////////////////////////////////////////////
1151//
1152//
1153
1154inline G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
1155{
1156 G4double dTheta = fRutherfordTheta - theta;
1157 G4double result = 0., argument = 0.;
1158
1159 if(std::abs(dTheta) < 0.001) result = 1.;
1160 else
1161 {
1162 argument = fProfileDelta*dTheta;
1163 result = CLHEP::pi*argument;
1164 result /= std::sinh(CLHEP::pi*argument);
1165 }
1166 return result;
1167}
1168
1169/////////////////////////////////////////////////////////////////
1170//
1171//
1172
1173inline G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
1174{
1175 G4double twosigma = 2.*fCoulombPhase0;
1176 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1177 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1178 twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
1179
1180 twosigma *= fCofPhase;
1181
1182 G4complex z = G4complex(0., twosigma);
1183
1184 return std::exp(z);
1185}
1186
1187/////////////////////////////////////////////////////////////////
1188//
1189//
1190
1191inline G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
1192{
1193 G4double twosigma = 2.*fCoulombPhase0;
1194 twosigma -= fZommerfeld*std::log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
1195 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
1196 twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
1197
1198 twosigma *= fCofPhase;
1199
1200 G4complex z = G4complex(0., twosigma);
1201
1202 return std::exp(z);
1203}
1204
1205/////////////////////////////////////////////////////////////////
1206//
1207//
1208
1209
1210inline G4complex G4NuclNuclDiffuseElastic::GammaLess(G4double theta)
1211{
1212 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1213 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1214
1215 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1216 G4double kappa = u/std::sqrt(CLHEP::pi);
1217 G4double dTheta = theta - fRutherfordTheta;
1218 u *= dTheta;
1219 G4double u2 = u*u;
1220 G4double u2m2p3 = u2*2./3.;
1221
1222 G4complex im = G4complex(0.,1.);
1223 G4complex order = G4complex(u,u);
1224 order /= std::sqrt(2.);
1225
1226 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1227 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1228 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1229 G4complex out = gamma*(1. - a1*dTheta) - a0;
1230
1231 return out;
1232}
1233
1234/////////////////////////////////////////////////////////////////
1235//
1236//
1237
1238inline G4complex G4NuclNuclDiffuseElastic::GammaMore(G4double theta)
1239{
1240 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1241 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1242
1243 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1244 G4double kappa = u/std::sqrt(CLHEP::pi);
1245 G4double dTheta = theta - fRutherfordTheta;
1246 u *= dTheta;
1247 G4double u2 = u*u;
1248 G4double u2m2p3 = u2*2./3.;
1249
1250 G4complex im = G4complex(0.,1.);
1251 G4complex order = G4complex(u,u);
1252 order /= std::sqrt(2.);
1253 G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1254 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1255 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1256 G4complex out = -gamma*(1. - a1*dTheta) - a0;
1257
1258 return out;
1259}
1260
1261/////////////////////////////////////////////////////////////////
1262//
1263//
1264
1265inline G4complex G4NuclNuclDiffuseElastic::AmplitudeNear(G4double theta)
1266{
1267 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1268 G4complex out = G4complex(kappa/fWaveVector,0.);
1269
1270 out *= PhaseNear(theta);
1271
1272 if( theta <= fRutherfordTheta )
1273 {
1274 out *= GammaLess(theta) + ProfileNear(theta);
1275 // out *= GammaMore(theta) + ProfileNear(theta);
1276 out += CoulombAmplitude(theta);
1277 }
1278 else
1279 {
1280 out *= GammaMore(theta) + ProfileNear(theta);
1281 // out *= GammaLess(theta) + ProfileNear(theta);
1282 }
1283 return out;
1284}
1285
1286/////////////////////////////////////////////////////////////////
1287//
1288//
1289
1290inline G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
1291{
1292 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1293 G4complex out = G4complex(kappa/fWaveVector,0.);
1294 out *= ProfileFar(theta);
1295 out *= PhaseFar(theta);
1296 return out;
1297}
1298
1299
1300/////////////////////////////////////////////////////////////////
1301//
1302//
1303
1304inline G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
1305{
1306
1307 G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
1308 // G4complex out = AmplitudeNear(theta);
1309 // G4complex out = AmplitudeFar(theta);
1310 return out;
1311}
1312
1313/////////////////////////////////////////////////////////////////
1314//
1315//
1316
1317inline G4double G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
1318{
1319 G4complex out = Amplitude(theta);
1320 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1321 return mod2;
1322}
1323
1324/////////////////////////////////////////////////////////////////
1325//
1326//
1327
1328inline G4complex G4NuclNuclDiffuseElastic::AmplitudeSim(G4double theta)
1329{
1330 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1331 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1332 G4double sindTheta = std::sin(dTheta);
1333 G4double persqrt2 = std::sqrt(0.5);
1334
1335 G4complex order = G4complex(persqrt2,persqrt2);
1336 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1337 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1338
1339 G4complex out;
1340
1341 if ( theta <= fRutherfordTheta )
1342 {
1343 out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1344 }
1345 else
1346 {
1347 out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1348 }
1349
1350 out *= CoulombAmplitude(theta);
1351
1352 return out;
1353}
1354
1355/////////////////////////////////////////////////////////////////
1356//
1357//
1358
1359inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
1360{
1361 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1362 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1363 G4double sindTheta = std::sin(dTheta);
1364
1365 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1366 // G4cout<<"order = "<<order<<G4endl;
1367 G4double cosFresnel = 0.5 - GetCint(order);
1368 G4double sinFresnel = 0.5 - GetSint(order);
1369
1370 G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1371
1372 return out;
1373}
1374
1375/////////////////////////////////////////////////////////////////
1376//
1377// The ratio el/ruth for Fresnel smooth nucleus profile
1378
1379inline G4double G4NuclNuclDiffuseElastic::GetRatioGen(G4double theta)
1380{
1381 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1382 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1383 G4double sindTheta = std::sin(dTheta);
1384
1385 G4double prof = Profile(theta);
1386 G4double prof2 = prof*prof;
1387 // G4double profmod = std::abs(prof);
1388 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1389
1390 order = std::abs(order); // since sin changes sign!
1391 // G4cout<<"order = "<<order<<G4endl;
1392
1393 G4double cosFresnel = GetCint(order);
1394 G4double sinFresnel = GetSint(order);
1395
1396 G4double out;
1397
1398 if ( theta <= fRutherfordTheta )
1399 {
1400 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1401 out += ( cosFresnel + sinFresnel - 1. )*prof;
1402 }
1403 else
1404 {
1405 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1406 }
1407
1408 return out;
1409}
1410
1411/////////////////////////////////////////////////////////////////
1412//
1413// The xsc for Fresnel smooth nucleus profile
1414
1415inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
1416{
1417 G4double ratio = GetRatioGen(theta);
1418 G4double ruthXsc = GetRutherfordXsc(theta);
1419 G4double xsc = ratio*ruthXsc;
1420 return xsc;
1421}
1422
1423/////////////////////////////////////////////////////////////////
1424//
1425// The xsc for Fresnel smooth nucleus profile for integration
1426
1427inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
1428{
1429 G4double theta = std::sqrt(alpha);
1430 G4double xsc = GetFresnelDiffuseXsc(theta);
1431 return xsc;
1432}
1433
1434
1435
1436/////////////////////////////////////////////////////////////////
1437//
1438//
1439
1440inline G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
1441{
1442 G4complex out = AmplitudeSim(theta);
1443 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1444 return mod2;
1445}
1446
1447/////////////////////////////////////////////////////////////////
1448//
1449//
1450
1451inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGla(G4double theta)
1452{
1453 G4int n;
1454 G4double T12b, b, b2; // cosTheta = std::cos(theta);
1455 G4complex out = G4complex(0.,0.), shiftC, shiftN;
1456 G4complex im = G4complex(0.,1.);
1457
1458 for( n = 0; n < fMaxL; n++)
1459 {
1460 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1461 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1462 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1463 b2 = b*b;
1464 T12b = fSumSigma*std::exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1465 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1466 out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1467 }
1468 out /= 2.*im*fWaveVector;
1469 out += CoulombAmplitude(theta);
1470 return out;
1471}
1472
1473/////////////////////////////////////////////////////////////////
1474//
1475//
1476
1477inline G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
1478{
1479 G4complex out = AmplitudeGla(theta);
1480 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1481 return mod2;
1482}
1483
1484
1485/////////////////////////////////////////////////////////////////
1486//
1487//
1488
1489inline G4complex G4NuclNuclDiffuseElastic::AmplitudeGG(G4double theta)
1490{
1491 G4int n;
1492 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1493 G4double sinThetaH2 = sinThetaH*sinThetaH;
1494 G4complex out = G4complex(0.,0.);
1495 G4complex im = G4complex(0.,1.);
1496
1497 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1498 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1499
1500 aTemp = a;
1501
1502 for( n = 1; n < fMaxL; n++)
1503 {
1504 T12b = aTemp*std::exp(-b2/n)/n;
1505 aTemp *= a;
1506 out += T12b;
1507 G4cout<<"out = "<<out<<G4endl;
1508 }
1509 out *= -4.*im*fWaveVector/CLHEP::pi;
1510 out += CoulombAmplitude(theta);
1511 return out;
1512}
1513
1514/////////////////////////////////////////////////////////////////
1515//
1516//
1517
1518inline G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
1519{
1520 G4complex out = AmplitudeGG(theta);
1521 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1522 return mod2;
1523}
1524
1525
1526///////////////////////////////////////////////////////////////////////////////
1527//
1528// Test for given particle and element table of momentum, angle probability.
1529// For the partMom in CMS.
1530
1531inline void G4NuclNuclDiffuseElastic::InitParameters(const G4ParticleDefinition* theParticle,
1532 G4double partMom, G4double Z, G4double A)
1533{
1534 fAtomicNumber = Z; // atomic number
1535 fAtomicWeight = A; // number of nucleons
1536
1537 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1538 G4double A1 = G4double( theParticle->GetBaryonNumber() );
1539 fNuclearRadius1 = CalculateNuclearRad(A1);
1540 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1541 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1542
1543 G4double a = 0.;
1544 G4double z = theParticle->GetPDGCharge();
1545 G4double m1 = theParticle->GetPDGMass();
1546
1547 fWaveVector = partMom/CLHEP::hbarc;
1548
1549 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1550 G4cout<<"kR = "<<lambda<<G4endl;
1551
1552 if( z )
1553 {
1554 a = partMom/m1; // beta*gamma for m1
1555 fBeta = a/std::sqrt(1+a*a);
1556 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1557 fRutherfordRatio = fZommerfeld/fWaveVector;
1558 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1559 }
1560 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1561 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1562 G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1563 fProfileDelta = fCofDelta*fProfileLambda;
1564 fProfileAlpha = fCofAlpha*fProfileLambda;
1565
1566 CalculateCoulombPhaseZero();
1567 CalculateRutherfordAnglePar();
1568
1569 return;
1570}
1571///////////////////////////////////////////////////////////////////////////////
1572//
1573// Test for given particle and element table of momentum, angle probability.
1574// For the partMom in CMS.
1575
1576inline void G4NuclNuclDiffuseElastic::InitDynParameters(const G4ParticleDefinition* theParticle,
1577 G4double partMom)
1578{
1579 G4double a = 0.;
1580 G4double z = theParticle->GetPDGCharge();
1581 G4double m1 = theParticle->GetPDGMass();
1582
1583 fWaveVector = partMom/CLHEP::hbarc;
1584
1585 G4double lambda = fCofLambda*fWaveVector*fNuclearRadius;
1586
1587 if( z )
1588 {
1589 a = partMom/m1; // beta*gamma for m1
1590 fBeta = a/std::sqrt(1+a*a);
1591 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1592 fRutherfordRatio = fZommerfeld/fWaveVector;
1593 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1594 }
1595 fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1596 fProfileDelta = fCofDelta*fProfileLambda;
1597 fProfileAlpha = fCofAlpha*fProfileLambda;
1598
1599 CalculateCoulombPhaseZero();
1600 CalculateRutherfordAnglePar();
1601
1602 return;
1603}
1604
1605
1606///////////////////////////////////////////////////////////////////////////////
1607//
1608// Test for given particle and element table of momentum, angle probability.
1609// For the partMom in CMS.
1610
1611inline void G4NuclNuclDiffuseElastic::InitParametersGla(const G4DynamicParticle* aParticle,
1612 G4double partMom, G4double Z, G4double A)
1613{
1614 fAtomicNumber = Z; // target atomic number
1615 fAtomicWeight = A; // target number of nucleons
1616
1617 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1618 G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1619 fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1620 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1621
1622
1623 G4double a = 0., kR12;
1624 G4double z = aParticle->GetDefinition()->GetPDGCharge();
1625 G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1626
1627 fWaveVector = partMom/CLHEP::hbarc;
1628
1629 G4double pN = A1 - z;
1630 if( pN < 0. ) pN = 0.;
1631
1632 G4double tN = A - Z;
1633 if( tN < 0. ) tN = 0.;
1634
1635 G4double pTkin = aParticle->GetKineticEnergy();
1636 pTkin /= A1;
1637
1638
1639 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1640 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1641
1642 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1643 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiusSquare/CLHEP::millibarn<<" mb"<<G4endl;
1644 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1645 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1646 fMaxL = (G4int(kR12)+1)*4;
1647 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1648
1649 if( z )
1650 {
1651 a = partMom/m1; // beta*gamma for m1
1652 fBeta = a/std::sqrt(1+a*a);
1653 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1654 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1655 }
1656
1657 CalculateCoulombPhaseZero();
1658
1659
1660 return;
1661}
1662
1663
1664/////////////////////////////////////////////////////////////////////////////////////
1665//
1666// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1667// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1668// projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1669
1670inline G4double
1671G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
1672 G4double pTkin,
1673 G4ParticleDefinition* tParticle)
1674{
1675 G4double xsection(0), /*Delta,*/ A0, B0;
1676 G4double hpXsc(0);
1677 G4double hnXsc(0);
1678
1679
1680 G4double targ_mass = tParticle->GetPDGMass();
1681 G4double proj_mass = pParticle->GetPDGMass();
1682
1683 G4double proj_energy = proj_mass + pTkin;
1684 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1685
1686 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1687
1688 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1689 proj_momentum /= CLHEP::GeV;
1690 proj_energy /= CLHEP::GeV;
1691 proj_mass /= CLHEP::GeV;
1692 G4double logS = std::log(sMand);
1693
1694 // General PDG fit constants
1695
1696
1697 // fEtaRatio=Re[f(0)]/Im[f(0)]
1698
1699 if( proj_momentum >= 1.2 )
1700 {
1701 fEtaRatio = 0.13*(logS - 5.8579332)*std::pow(sMand,-0.18);
1702 }
1703 else if( proj_momentum >= 0.6 )
1704 {
1705 fEtaRatio = -75.5*(std::pow(proj_momentum,0.25)-0.95)/
1706 (std::pow(3*proj_momentum,2.2)+1);
1707 }
1708 else
1709 {
1710 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1711 }
1712 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1713
1714 // xsc
1715
1716 if( proj_momentum >= 10. ) // high energy: pp = nn = np
1717 // if( proj_momentum >= 2.)
1718 {
1719 //Delta = 1.;
1720
1721 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1722
1723 if( proj_momentum >= 10.)
1724 {
1725 B0 = 7.5;
1726 A0 = 100. - B0*std::log(3.0e7);
1727
1728 xsection = A0 + B0*std::log(proj_energy) - 11
1729 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
1730 0.93827*0.93827,-0.165); // mb
1731 }
1732 }
1733 else // low energy pp = nn != np
1734 {
1735 if(pParticle == tParticle) // pp or nn // nn to be pp
1736 {
1737 if( proj_momentum < 0.73 )
1738 {
1739 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
1740 }
1741 else if( proj_momentum < 1.05 )
1742 {
1743 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
1744 (std::log(proj_momentum/0.73));
1745 }
1746 else // if( proj_momentum < 10. )
1747 {
1748 hnXsc = 39.0 +
1749 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
1750 }
1751 xsection = hnXsc;
1752 }
1753 else // pn to be np
1754 {
1755 if( proj_momentum < 0.8 )
1756 {
1757 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
1758 }
1759 else if( proj_momentum < 1.4 )
1760 {
1761 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
1762 }
1763 else // if( proj_momentum < 10. )
1764 {
1765 hpXsc = 33.3+
1766 20.8*(std::pow(proj_momentum,2.0)-1.35)/
1767 (std::pow(proj_momentum,2.50)+0.95);
1768 }
1769 xsection = hpXsc;
1770 }
1771 }
1772 xsection *= CLHEP::millibarn; // parametrised in mb
1773 G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1774 return xsection;
1775}
1776
1777////////////////////////////////////////////////////////////////////////////////////
1778//
1779//
1780
1781inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
1782 const G4double mt ,
1783 const G4double Plab )
1784{
1785 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1786 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1787
1788 return sMand;
1789}
1790
1791
1792
1793#endif
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4complex
std::complex< G4double > G4complex
Definition: G4Types.hh:69
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4NuclNuclDiffuseElastic::CalcMandelstamS
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
Definition: G4NuclNuclDiffuseElastic.hh:1781
G4NuclNuclDiffuseElastic::GetRatioSim
G4double GetRatioSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1359
G4NuclNuclDiffuseElastic::GetInvElasticXsc
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:203
G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
Definition: G4NuclNuclDiffuseElastic.cc:1227
G4NuclNuclDiffuseElastic::CalculateZommerfeld
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4NuclNuclDiffuseElastic.hh:568
G4NuclNuclDiffuseElastic::SampleThetaLab
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:1078
G4NuclNuclDiffuseElastic::CalculateParticleBeta
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4NuclNuclDiffuseElastic.hh:554
G4NuclNuclDiffuseElastic::TestAngleTable
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:1286
G4NuclNuclDiffuseElastic::GetCoulombTotalXsc
G4double GetCoulombTotalXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:665
G4NuclNuclDiffuseElastic::GetRatioGen
G4double GetRatioGen(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1379
G4NuclNuclDiffuseElastic::Profile
G4double Profile(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1154
G4NuclNuclDiffuseElastic::AmplitudeSim
G4complex AmplitudeSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1328
G4NuclNuclDiffuseElastic::GetDiffElasticSumProb
G4double GetDiffElasticSumProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:466
G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:340
G4NuclNuclDiffuseElastic::InitialiseOnFly
void InitialiseOnFly(G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:929
G4NuclNuclDiffuseElastic::GetCoulombElasticXsc
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:625
G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit
void SetRecoilKinEnergyLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:375
G4NuclNuclDiffuseElastic::AmplitudeFar
G4complex AmplitudeFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1290
G4NuclNuclDiffuseElastic::AmplitudeGla
G4complex AmplitudeGla(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1451
G4NuclNuclDiffuseElastic::GetErfComp
G4complex GetErfComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:878
G4NuclNuclDiffuseElastic::GammaLogB2n
G4complex GammaLogB2n(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.hh:746
G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc
G4double GetCoulombIntegralXsc(const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
Definition: G4NuclNuclDiffuseElastic.hh:690
G4NuclNuclDiffuseElastic::TestErfcInt
G4complex TestErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:866
G4NuclNuclDiffuseElastic::Amplitude
G4complex Amplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1304
G4NuclNuclDiffuseElastic::AmplitudeSimMod2
G4double AmplitudeSimMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1440
G4NuclNuclDiffuseElastic::GetInvElasticSumXsc
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:288
G4NuclNuclDiffuseElastic::GetLegendrePol
G4double GetLegendrePol(G4int n, G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:814
G4NuclNuclDiffuseElastic::GetIntegrandFunction
G4double GetIntegrandFunction(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:657
G4NuclNuclDiffuseElastic::AmplitudeNear
G4complex AmplitudeNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1265
G4NuclNuclDiffuseElastic::AmplitudeGG
G4complex AmplitudeGG(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1489
G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
Definition: G4NuclNuclDiffuseElastic.hh:1671
G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2
G4double CoulombAmplitudeMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1064
G4NuclNuclDiffuseElastic::IntegralElasticProb
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:673
G4NuclNuclDiffuseElastic::ProfileNear
G4double ProfileNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1117
G4NuclNuclDiffuseElastic::SampleTableThetaCMS
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:806
G4NuclNuclDiffuseElastic::GetErfInt
G4complex GetErfInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:987
G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc
G4double GetFresnelDiffuseXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1415
G4NuclNuclDiffuseElastic::SampleT
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:698
G4NuclNuclDiffuseElastic::InitParametersGla
void InitParametersGla(const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:1611
G4NuclNuclDiffuseElastic::PhaseNear
G4complex PhaseNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1173
G4NuclNuclDiffuseElastic::SetHEModelLowLimit
void SetHEModelLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:385
G4NuclNuclDiffuseElastic::GammaLogarithm
G4complex GammaLogarithm(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.hh:722
G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:253
G4NuclNuclDiffuseElastic::GetErfcSer
G4complex GetErfcSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:798
G4NuclNuclDiffuseElastic::CoulombAmplitude
G4complex CoulombAmplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1043
G4NuclNuclDiffuseElastic::GetScatteringAngle
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
Definition: G4NuclNuclDiffuseElastic.cc:1035
G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA
G4double GetDiffElasticSumProbA(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.cc:561
G4NuclNuclDiffuseElastic::SetPlabLowLimit
void SetPlabLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:380
G4NuclNuclDiffuseElastic::AmplitudeMod2
G4double AmplitudeMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1317
G4NuclNuclDiffuseElastic::AmplitudeGlaMod2
G4double AmplitudeGlaMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1477
G4NuclNuclDiffuseElastic::TestErfcSer
G4complex TestErfcSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:854
G4NuclNuclDiffuseElastic::GetErfcComp
G4complex GetErfcComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:788
G4NuclNuclDiffuseElastic::GetRutherfordXsc
G4double GetRutherfordXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:648
G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
Definition: G4NuclNuclDiffuseElastic.cc:1166
G4NuclNuclDiffuseElastic::SampleTableT
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:791
G4NuclNuclDiffuseElastic::SampleThetaCMS
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:712
G4NuclNuclDiffuseElastic::SetQModelLowLimit
void SetQModelLowLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:390
G4NuclNuclDiffuseElastic::GetErfSer
G4complex GetErfSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:938
G4NuclNuclDiffuseElastic::SetLowestEnergyLimit
void SetLowestEnergyLimit(G4double value)
Definition: G4NuclNuclDiffuseElastic.hh:395
G4NuclNuclDiffuseElastic::TestErfcComp
G4complex TestErfcComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:842
G4NuclNuclDiffuseElastic::AmplitudeGGMod2
G4double AmplitudeGGMod2(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1518
G4NuclNuclDiffuseElastic::InitDynParameters
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
Definition: G4NuclNuclDiffuseElastic.hh:1576
G4NuclNuclDiffuseElastic::CalculateNuclearRad
G4double CalculateNuclearRad(G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:594
G4NuclNuclDiffuseElastic::InitParameters
void InitParameters(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:1531
G4NuclNuclDiffuseElastic::GetErfcInt
G4complex GetErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:808
G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:182
G4NuclNuclDiffuseElastic::SampleInvariantT
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
Definition: G4NuclNuclDiffuseElastic.cc:766
G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc
G4double GetFresnelIntegrandXsc(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.hh:1427
G4NuclNuclDiffuseElastic::PhaseFar
G4complex PhaseFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1191
G4NuclNuclDiffuseElastic::ProfileFar
G4double ProfileFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1138
G4NuclNuclDiffuseElastic::GammaMore
G4complex GammaMore(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1238
G4NuclNuclDiffuseElastic::GetDiffElasticProb
G4double GetDiffElasticProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:389
G4NuclNuclDiffuseElastic::CalculateAm
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:579
G4NuclNuclDiffuseElastic::GammaLess
G4complex GammaLess(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1210
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119