| Geant4 Cross Reference |
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 //
28 //
29 // G4 Model: optical elastic scattering with 4-momentum balance
30 //
31 // Class Description
32 // Final state production model for nucleus-nucleus elastic scattering;
33 // Coulomb amplitude is not considered as correction
34 // (as in G4DiffuseElastic)
35 // Class Description - End
36 //
37 //
38 // 17.03.09 V. Grichine implementation for Coulomb elastic scattering
39
40
41 #ifndef G4NuclNuclDiffuseElastic_h
42 #define G4NuclNuclDiffuseElastic_h 1
43
44 #include <complex>
45 #include <CLHEP/Units/PhysicalConstants.h>
46 #include "globals.hh"
47 #include "G4Integrator.hh"
48 #include "G4HadronElastic.hh"
49 #include "G4HadProjectile.hh"
50 #include "G4Nucleus.hh"
51
52 #include "G4Exp.hh"
53 #include "G4Log.hh"
54 #include "G4Pow.hh"
55
56
57 class G4ParticleDefinition;
58 class G4PhysicsTable;
59 class G4PhysicsLogVector;
60
61 class G4NuclNuclDiffuseElastic : public G4HadronElastic // G4HadronicInteraction
62 {
63 public:
64
65 G4NuclNuclDiffuseElastic();
66
67 // G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle);
68
69 virtual ~G4NuclNuclDiffuseElastic();
70
71 void Initialise();
72
73 void InitialiseOnFly(G4double Z, G4double A);
74
75 void BuildAngleTable();
76
77 virtual G4double SampleInvariantT(const G4ParticleDefinition* p,
78 G4double plab,
79 G4int Z, G4int A);
80
81 void SetPlabLowLimit(G4double value);
82
83 void SetHEModelLowLimit(G4double value);
84
85 void SetQModelLowLimit(G4double value);
86
87 void SetLowestEnergyLimit(G4double value);
88
89 void SetRecoilKinEnergyLimit(G4double value);
90
91 G4double SampleT(const G4ParticleDefinition* aParticle,
92 G4double p, G4double A);
93
94 G4double SampleTableT(const G4ParticleDefinition* aParticle,
95 G4double p, G4double Z, G4double A);
96
97 G4double SampleThetaCMS( const G4ParticleDefinition* aParticle, G4double p, G4double A);
98
99 G4double SampleCoulombMuCMS( const G4ParticleDefinition* aParticle, G4double p);
100
101 G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p,
102 G4double Z, G4double A);
103
104 G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
105
106 G4double SampleThetaLab(const G4HadProjectile* aParticle,
107 G4double tmass, G4double A);
108
109 G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
110 G4double theta,
111 G4double momentum,
112 G4double A );
113
114 G4double GetInvElasticXsc( const G4ParticleDefinition* particle,
115 G4double theta,
116 G4double momentum,
117 G4double A, G4double Z );
118
119 G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
120 G4double theta,
121 G4double momentum,
122 G4double A, G4double Z );
123
124 G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle,
125 G4double tMand,
126 G4double momentum,
127 G4double A, G4double Z );
128
129 G4double IntegralElasticProb( const G4ParticleDefinition* particle,
130 G4double theta,
131 G4double momentum,
132 G4double A );
133
134
135 G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle,
136 G4double theta,
137 G4double momentum,
138 G4double Z );
139
140 G4double GetRutherfordXsc( G4double theta );
141
142 G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
143 G4double tMand,
144 G4double momentum,
145 G4double A, G4double Z );
146
147 G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,
148 G4double momentum, G4double Z );
149
150 G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,
151 G4double momentum, G4double Z,
152 G4double theta1, G4double theta2 );
153
154
155 G4double CalculateParticleBeta( const G4ParticleDefinition* particle,
156 G4double momentum );
157
158 G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
159
160 G4double CalculateAm( G4double momentum, G4double n, G4double Z);
161
162 G4double CalculateNuclearRad( G4double A);
163
164 G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle,
165 G4double tmass, G4double thetaCMS);
166
167 G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle,
168 G4double tmass, G4double thetaLab);
169
170 void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
171 G4double Z, G4double A);
172
173
174
175 G4double BesselJzero(G4double z);
176 G4double BesselJone(G4double z);
177 G4double DampFactor(G4double z);
178 G4double BesselOneByArg(G4double z);
179
180 G4double GetDiffElasticProb(G4double theta);
181 G4double GetDiffElasticSumProb(G4double theta);
182 G4double GetDiffElasticSumProbA(G4double alpha);
183 G4double GetIntegrandFunction(G4double theta);
184
185 G4double GetNuclearRadius(){return fNuclearRadius;};
186
187
188 // Technical math functions for strong Coulomb contribution
189
190 G4complex GammaLogarithm(G4complex xx);
191 G4complex GammaLogB2n(G4complex xx);
192
193 G4double GetErf(G4double x);
194
195 G4double GetCosHaPit2(G4double t){return std::cos(CLHEP::halfpi*t*t);};
196 G4double GetSinHaPit2(G4double t){return std::sin(CLHEP::halfpi*t*t);};
197
198 G4double GetCint(G4double x);
199 G4double GetSint(G4double x);
200
201
202 G4complex GetErfcComp(G4complex z, G4int nMax);
203 G4complex GetErfcSer(G4complex z, G4int nMax);
204 G4complex GetErfcInt(G4complex z); // , G4int nMax);
205
206 G4complex GetErfComp(G4complex z, G4int nMax); // AandS algorithm != Ser, Int
207 G4complex GetErfSer(G4complex z, G4int nMax);
208
209 G4double GetExpCos(G4double x);
210 G4double GetExpSin(G4double x);
211 G4complex GetErfInt(G4complex z); // , G4int nMax);
212
213 G4double GetLegendrePol(G4int n, G4double x);
214
215 G4complex TestErfcComp(G4complex z, G4int nMax);
216 G4complex TestErfcSer(G4complex z, G4int nMax);
217 G4complex TestErfcInt(G4complex z); // , G4int nMax);
218
219 G4complex CoulombAmplitude(G4double theta);
220 G4double CoulombAmplitudeMod2(G4double theta);
221
222 void CalculateCoulombPhaseZero();
223 G4double CalculateCoulombPhase(G4int n);
224 void CalculateRutherfordAnglePar();
225
226 G4double ProfileNear(G4double theta);
227 G4double ProfileFar(G4double theta);
228 G4double Profile(G4double theta);
229
230 G4complex PhaseNear(G4double theta);
231 G4complex PhaseFar(G4double theta);
232
233 G4complex GammaLess(G4double theta);
234 G4complex GammaMore(G4double theta);
235
236 G4complex AmplitudeNear(G4double theta);
237 G4complex AmplitudeFar(G4double theta);
238
239 G4complex Amplitude(G4double theta);
240 G4double AmplitudeMod2(G4double theta);
241
242 G4complex AmplitudeSim(G4double theta);
243 G4double AmplitudeSimMod2(G4double theta);
244
245 G4double GetRatioSim(G4double theta);
246 G4double GetRatioGen(G4double theta);
247
248 G4double GetFresnelDiffuseXsc(G4double theta);
249 G4double GetFresnelIntegrandXsc(G4double alpha);
250
251
252 G4complex AmplitudeGla(G4double theta);
253 G4double AmplitudeGlaMod2(G4double theta);
254
255 G4complex AmplitudeGG(G4double theta);
256 G4double AmplitudeGGMod2(G4double theta);
257
258 void InitParameters(const G4ParticleDefinition* theParticle,
259 G4double partMom, G4double Z, G4double A);
260
261 void InitDynParameters(const G4ParticleDefinition* theParticle,
262 G4double partMom);
263
264 void InitParametersGla(const G4DynamicParticle* aParticle,
265 G4double partMom, G4double Z, G4double A);
266
267 G4double GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
268 G4double pTkin,
269 G4ParticleDefinition* tParticle);
270
271 G4double CalcMandelstamS( const G4double mp , const G4double mt ,
272 const G4double Plab );
273
274 G4double GetProfileLambda(){return fProfileLambda;};
275
276 inline void SetProfileLambda(G4double pl) {fProfileLambda = pl;};
277 inline void SetProfileDelta(G4double pd) {fProfileDelta = pd;};
278 inline void SetProfileAlpha(G4double pa){fProfileAlpha = pa;};
279 inline void SetCofLambda(G4double pa){fCofLambda = pa;};
280
281 inline void SetCofAlpha(G4double pa){fCofAlpha = pa;};
282 inline void SetCofAlphaMax(G4double pa){fCofAlphaMax = pa;};
283 inline void SetCofAlphaCoulomb(G4double pa){fCofAlphaCoulomb = pa;};
284
285 inline void SetCofDelta(G4double pa){fCofDelta = pa;};
286 inline void SetCofPhase(G4double pa){fCofPhase = pa;};
287 inline void SetCofFar(G4double pa){fCofFar = pa;};
288 inline void SetEtaRatio(G4double pa){fEtaRatio = pa;};
289 inline void SetMaxL(G4int l){fMaxL = l;};
290 inline void SetNuclearRadiusCof(G4double r){fNuclearRadiusCof = r;};
291
292 inline G4double GetCofAlphaMax(){return fCofAlphaMax;};
293 inline G4double GetCofAlphaCoulomb(){return fCofAlphaCoulomb;};
294
295 private:
296
297 G4ParticleDefinition* theProton;
298 G4ParticleDefinition* theNeutron;
299 G4ParticleDefinition* theDeuteron;
300 G4ParticleDefinition* theAlpha;
301
302 const G4ParticleDefinition* thePionPlus;
303 const G4ParticleDefinition* thePionMinus;
304
305 G4double lowEnergyRecoilLimit;
306 G4double lowEnergyLimitHE;
307 G4double lowEnergyLimitQ;
308 G4double lowestEnergyLimit;
309 G4double plabLowLimit;
310
311 G4int fEnergyBin;
312 G4int fAngleBin;
313
314 G4PhysicsLogVector* fEnergyVector;
315 G4PhysicsTable* fAngleTable;
316 std::vector<G4PhysicsTable*> fAngleBank;
317
318 std::vector<G4double> fElementNumberVector;
319 std::vector<G4String> fElementNameVector;
320
321 const G4ParticleDefinition* fParticle;
322
323 G4double fWaveVector;
324 G4double fAtomicWeight;
325 G4double fAtomicNumber;
326
327 G4double fNuclearRadius1;
328 G4double fNuclearRadius2;
329 G4double fNuclearRadius;
330 G4double fNuclearRadiusSquare;
331 G4double fNuclearRadiusCof;
332
333 G4double fBeta;
334 G4double fZommerfeld;
335 G4double fRutherfordRatio;
336 G4double fAm;
337 G4bool fAddCoulomb;
338
339 G4double fCoulombPhase0;
340 G4double fHalfRutThetaTg;
341 G4double fHalfRutThetaTg2;
342 G4double fRutherfordTheta;
343
344 G4double fProfileLambda;
345 G4double fProfileDelta;
346 G4double fProfileAlpha;
347
348 G4double fCofLambda;
349 G4double fCofAlpha;
350 G4double fCofDelta;
351 G4double fCofPhase;
352 G4double fCofFar;
353
354 G4double fCofAlphaMax;
355 G4double fCofAlphaCoulomb;
356
357 G4int fMaxL;
358 G4double fSumSigma;
359 G4double fEtaRatio;
360
361 G4double fReZ;
362 G4double fCoulombMuC;
363
364 };
365
366
367 inline void G4NuclNuclDiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
368 {
369 lowEnergyRecoilLimit = value;
370 }
371
372 inline void G4NuclNuclDiffuseElastic::SetPlabLowLimit(G4double value)
373 {
374 plabLowLimit = value;
375 }
376
377 inline void G4NuclNuclDiffuseElastic::SetHEModelLowLimit(G4double value)
378 {
379 lowEnergyLimitHE = value;
380 }
381
382 inline void G4NuclNuclDiffuseElastic::SetQModelLowLimit(G4double value)
383 {
384 lowEnergyLimitQ = value;
385 }
386
387 inline void G4NuclNuclDiffuseElastic::SetLowestEnergyLimit(G4double value)
388 {
389 lowestEnergyLimit = value;
390 }
391
392 ////////////////////////////////////////////////////////////////////
393 //
394 // damp factor in diffraction x/sh(x), x was already *pi
395
396 inline G4double G4NuclNuclDiffuseElastic::DampFactor(G4double x)
397 {
398 G4double df;
399 G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
400
401 // x *= pi;
402
403 if( std::fabs(x) < 0.01 )
404 {
405 df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
406 }
407 else
408 {
409 df = x/std::sinh(x);
410 }
411 return df;
412 }
413
414
415 ////////////////////////////////////////////////////////////////////
416 //
417 // return J1(x)/x with special case for small x
418
419 inline G4double G4NuclNuclDiffuseElastic::BesselOneByArg(G4double x)
420 {
421 G4double x2, result;
422
423 if( std::fabs(x) < 0.01 )
424 {
425 x *= 0.5;
426 x2 = x*x;
427 result = 2. - x2 + x2*x2/6.;
428 }
429 else
430 {
431 result = BesselJone(x)/x;
432 }
433 return result;
434 }
435
436 ////////////////////////////////////////////////////////////////////
437 //
438 // return particle beta
439
440 inline G4double
441 G4NuclNuclDiffuseElastic::CalculateParticleBeta(const G4ParticleDefinition* particle,
442 G4double momentum)
443 {
444 G4double mass = particle->GetPDGMass();
445 G4double a = momentum/mass;
446 fBeta = a/std::sqrt(1+a*a);
447
448 return fBeta;
449 }
450
451 ////////////////////////////////////////////////////////////////////
452 //
453 // return Zommerfeld parameter for Coulomb scattering
454
455 inline G4double
456 G4NuclNuclDiffuseElastic::CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2 )
457 {
458 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
459
460 return fZommerfeld;
461 }
462
463 ////////////////////////////////////////////////////////////////////
464 //
465 // return Wentzel correction for Coulomb scattering
466
467 inline G4double
468 G4NuclNuclDiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
469 {
470 G4double k = momentum/CLHEP::hbarc;
471 G4double ch = 1.13 + 3.76*n*n;
472 G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
473 G4double zn2 = zn*zn;
474 fAm = ch/zn2;
475
476 return fAm;
477 }
478
479 ////////////////////////////////////////////////////////////////////
480 //
481 // calculate nuclear radius for different atomic weights using different approximations
482
483 inline G4double G4NuclNuclDiffuseElastic::CalculateNuclearRad( G4double A)
484 {
485 G4double r0 = 1.*CLHEP::fermi, radius;
486 // r0 *= 1.12;
487 // r0 *= 1.44;
488 r0 *= fNuclearRadiusCof;
489
490 /*
491 if( A < 50. )
492 {
493 if( A > 10. ) r0 = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08*fermi;
494 else r0 = 1.1*CLHEP::fermi;
495
496 radius = r0*G4Pow::GetInstance()->A13(A);
497 }
498 else
499 {
500 r0 = 1.7*CLHEP::fermi; // 1.7*fermi;
501
502 radius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
503 }
504 */
505 radius = r0*G4Pow::GetInstance()->A13(A);
506
507 return radius;
508 }
509
510 ////////////////////////////////////////////////////////////////////
511 //
512 // return Coulomb scattering differential xsc with Wentzel correction. Test function
513
514 inline G4double G4NuclNuclDiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle,
515 G4double theta,
516 G4double momentum,
517 G4double Z )
518 {
519 G4double sinHalfTheta = std::sin(0.5*theta);
520 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
521 G4double beta = CalculateParticleBeta( particle, momentum);
522 G4double z = particle->GetPDGCharge();
523 G4double n = CalculateZommerfeld( beta, z, Z );
524 G4double am = CalculateAm( momentum, n, Z);
525 G4double k = momentum/CLHEP::hbarc;
526 G4double ch = 0.5*n/k;
527 G4double ch2 = ch*ch;
528 G4double xsc = ch2/((sinHalfTheta2+am)*(sinHalfTheta2+am));
529
530 return xsc;
531 }
532
533 ////////////////////////////////////////////////////////////////////
534 //
535 // return Rutherford scattering differential xsc with Wentzel correction. For Sampling.
536
537 inline G4double G4NuclNuclDiffuseElastic::GetRutherfordXsc(G4double theta)
538 {
539 G4double sinHalfTheta = std::sin(0.5*theta);
540 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
541
542 G4double ch2 = fRutherfordRatio*fRutherfordRatio;
543
544 G4double xsc = ch2/(sinHalfTheta2+fAm)/(sinHalfTheta2+fAm);
545
546 return xsc;
547 }
548
549 ////////////////////////////////////////////////////////////////////
550 //
551 // return Coulomb scattering total xsc with Wentzel correction
552
553 inline G4double G4NuclNuclDiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,
554 G4double momentum, G4double Z )
555 {
556 G4double beta = CalculateParticleBeta( particle, momentum);
557 G4cout<<"beta = "<<beta<<G4endl;
558 G4double z = particle->GetPDGCharge();
559 G4double n = CalculateZommerfeld( beta, z, Z );
560 G4cout<<"fZomerfeld = "<<n<<G4endl;
561 G4double am = CalculateAm( momentum, n, Z);
562 G4cout<<"cof Am = "<<am<<G4endl;
563 G4double k = momentum/CLHEP::hbarc;
564 G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
565 G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
566 G4double ch = n/k;
567 G4double ch2 = ch*ch;
568 G4double xsc = ch2*CLHEP::pi/(am +am*am);
569
570 return xsc;
571 }
572
573 ////////////////////////////////////////////////////////////////////
574 //
575 // return Coulomb scattering xsc with Wentzel correction integrated between
576 // theta1 and < theta2
577
578 inline G4double
579 G4NuclNuclDiffuseElastic::GetCoulombIntegralXsc(const G4ParticleDefinition* particle,
580 G4double momentum, G4double Z,
581 G4double theta1, G4double theta2 )
582 {
583 G4double c1 = std::cos(theta1);
584 //G4cout<<"c1 = "<<c1<<G4endl;
585 G4double c2 = std::cos(theta2);
586 // G4cout<<"c2 = "<<c2<<G4endl;
587 G4double beta = CalculateParticleBeta( particle, momentum);
588 // G4cout<<"beta = "<<beta<<G4endl;
589 G4double z = particle->GetPDGCharge();
590 G4double n = CalculateZommerfeld( beta, z, Z );
591 // G4cout<<"fZomerfeld = "<<n<<G4endl;
592 G4double am = CalculateAm( momentum, n, Z);
593 // G4cout<<"cof Am = "<<am<<G4endl;
594 G4double k = momentum/CLHEP::hbarc;
595 // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
596 // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
597 G4double ch = n/k;
598 G4double ch2 = ch*ch;
599 am *= 2.;
600 G4double xsc = ch2*CLHEP::twopi*(c1-c2)/((1 - c1 + am)*(1 - c2 + am));
601
602 return xsc;
603 }
604
605 ///////////////////////////////////////////////////////////////////
606 //
607 // For the calculation of arg Gamma(z) one needs complex extension
608 // of ln(Gamma(z)) here is approximate algorithm
609
610 inline G4complex G4NuclNuclDiffuseElastic::GammaLogB2n(G4complex z)
611 {
612 G4complex z1 = 12.*z;
613 G4complex z2 = z*z;
614 G4complex z3 = z2*z;
615 G4complex z5 = z2*z3;
616 G4complex z7 = z2*z5;
617
618 z3 *= 360.;
619 z5 *= 1260.;
620 z7 *= 1680.;
621
622 G4complex result = (z-0.5)*std::log(z)-z+0.5*G4Log(CLHEP::twopi);
623 result += 1./z1 - 1./z3 +1./z5 -1./z7;
624 return result;
625 }
626
627 /////////////////////////////////////////////////////////////////
628 //
629 //
630
631 inline G4double G4NuclNuclDiffuseElastic::GetErf(G4double x)
632 {
633 G4double t, z, tmp, result;
634
635 z = std::fabs(x);
636 t = 1.0/(1.0+0.5*z);
637
638 tmp = t*std::exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
639 t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
640 t*(-0.82215223+t*0.17087277)))))))));
641
642 if( x >= 0.) result = 1. - tmp;
643 else result = 1. + tmp;
644
645 return result;
646 }
647
648 /////////////////////////////////////////////////////////////////
649 //
650 //
651
652 inline G4complex G4NuclNuclDiffuseElastic::GetErfcComp(G4complex z, G4int nMax)
653 {
654 G4complex erfcz = 1. - GetErfComp( z, nMax);
655 return erfcz;
656 }
657
658 /////////////////////////////////////////////////////////////////
659 //
660 //
661
662 inline G4complex G4NuclNuclDiffuseElastic::GetErfcSer(G4complex z, G4int nMax)
663 {
664 G4complex erfcz = 1. - GetErfSer( z, nMax);
665 return erfcz;
666 }
667
668 /////////////////////////////////////////////////////////////////
669 //
670 //
671
672 inline G4complex G4NuclNuclDiffuseElastic::GetErfcInt(G4complex z) // , G4int nMax)
673 {
674 G4complex erfcz = 1. - GetErfInt( z); // , nMax);
675 return erfcz;
676 }
677
678
679 /////////////////////////////////////////////////////////////////
680 //
681 //
682
683 inline G4complex G4NuclNuclDiffuseElastic::TestErfcComp(G4complex z, G4int nMax)
684 {
685 G4complex miz = G4complex( z.imag(), -z.real() );
686 G4complex erfcz = 1. - GetErfComp( miz, nMax);
687 G4complex w = std::exp(-z*z)*erfcz;
688 return w;
689 }
690
691 /////////////////////////////////////////////////////////////////
692 //
693 //
694
695 inline G4complex G4NuclNuclDiffuseElastic::TestErfcSer(G4complex z, G4int nMax)
696 {
697 G4complex miz = G4complex( z.imag(), -z.real() );
698 G4complex erfcz = 1. - GetErfSer( miz, nMax);
699 G4complex w = std::exp(-z*z)*erfcz;
700 return w;
701 }
702
703 /////////////////////////////////////////////////////////////////
704 //
705 //
706
707 inline G4complex G4NuclNuclDiffuseElastic::TestErfcInt(G4complex z) // , G4int nMax)
708 {
709 G4complex miz = G4complex( z.imag(), -z.real() );
710 G4complex erfcz = 1. - GetErfInt( miz); // , nMax);
711 G4complex w = std::exp(-z*z)*erfcz;
712 return w;
713 }
714
715 /////////////////////////////////////////////////////////////////
716 //
717 //
718
719 inline G4complex G4NuclNuclDiffuseElastic::GetErfSer(G4complex z, G4int nMax)
720 {
721 G4int n;
722 G4double a =1., b = 1., tmp;
723 G4complex sum = z, d = z;
724
725 for( n = 1; n <= nMax; n++)
726 {
727 a *= 2.;
728 b *= 2.*n +1.;
729 d *= z*z;
730
731 tmp = a/b;
732
733 sum += tmp*d;
734 }
735 sum *= 2.*std::exp(-z*z)/std::sqrt(CLHEP::pi);
736
737 return sum;
738 }
739
740 /////////////////////////////////////////////////////////////////////
741
742 inline G4double G4NuclNuclDiffuseElastic::GetExpCos(G4double x)
743 {
744 G4double result;
745
746 result = G4Exp(x*x-fReZ*fReZ);
747 result *= std::cos(2.*x*fReZ);
748 return result;
749 }
750
751 /////////////////////////////////////////////////////////////////////
752
753 inline G4double G4NuclNuclDiffuseElastic::GetExpSin(G4double x)
754 {
755 G4double result;
756
757 result = G4Exp(x*x-fReZ*fReZ);
758 result *= std::sin(2.*x*fReZ);
759 return result;
760 }
761
762
763 /////////////////////////////////////////////////////////////////
764 //
765 //
766
767 inline G4double G4NuclNuclDiffuseElastic::GetCint(G4double x)
768 {
769 G4double out;
770
771 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
772
773 out= integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetCosHaPit2, 0., x );
774
775 return out;
776 }
777
778 /////////////////////////////////////////////////////////////////
779 //
780 //
781
782 inline G4double G4NuclNuclDiffuseElastic::GetSint(G4double x)
783 {
784 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
785
786 G4double out =
787 integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetSinHaPit2, 0., x );
788
789 return out;
790 }
791
792
793 /////////////////////////////////////////////////////////////////
794 //
795 //
796
797 inline G4complex G4NuclNuclDiffuseElastic::CoulombAmplitude(G4double theta)
798 {
799 G4complex ca;
800
801 G4double sinHalfTheta = std::sin(0.5*theta);
802 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
803 sinHalfTheta2 += fAm;
804
805 G4double order = 2.*fCoulombPhase0 - fZommerfeld*G4Log(sinHalfTheta2);
806 G4complex z = G4complex(0., order);
807 ca = std::exp(z);
808
809 ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
810
811 return ca;
812 }
813
814 /////////////////////////////////////////////////////////////////
815 //
816 //
817
818 inline G4double G4NuclNuclDiffuseElastic::CoulombAmplitudeMod2(G4double theta)
819 {
820 G4complex ca = CoulombAmplitude(theta);
821 G4double out = ca.real()*ca.real() + ca.imag()*ca.imag();
822
823 return out;
824 }
825
826 /////////////////////////////////////////////////////////////////
827 //
828 //
829
830
831 inline void G4NuclNuclDiffuseElastic::CalculateCoulombPhaseZero()
832 {
833 G4complex z = G4complex(1,fZommerfeld);
834 // G4complex gammalog = GammaLogarithm(z);
835 G4complex gammalog = GammaLogB2n(z);
836 fCoulombPhase0 = gammalog.imag();
837 }
838
839 /////////////////////////////////////////////////////////////////
840 //
841 //
842
843
844 inline G4double G4NuclNuclDiffuseElastic::CalculateCoulombPhase(G4int n)
845 {
846 G4complex z = G4complex(1. + n, fZommerfeld);
847 // G4complex gammalog = GammaLogarithm(z);
848 G4complex gammalog = GammaLogB2n(z);
849 return gammalog.imag();
850 }
851
852
853 /////////////////////////////////////////////////////////////////
854 //
855 //
856
857
858 inline void G4NuclNuclDiffuseElastic::CalculateRutherfordAnglePar()
859 {
860 fHalfRutThetaTg = fZommerfeld/fProfileLambda; // (fWaveVector*fNuclearRadius);
861 fRutherfordTheta = 2.*std::atan(fHalfRutThetaTg);
862 fHalfRutThetaTg2 = fHalfRutThetaTg*fHalfRutThetaTg;
863 // G4cout<<"fRutherfordTheta = "<<fRutherfordTheta/degree<<" degree"<<G4endl;
864
865 }
866
867 /////////////////////////////////////////////////////////////////
868 //
869 //
870
871 inline G4double G4NuclNuclDiffuseElastic::ProfileNear(G4double theta)
872 {
873 G4double dTheta = fRutherfordTheta - theta;
874 G4double result = 0., argument = 0.;
875
876 if(std::abs(dTheta) < 0.001) result = fProfileAlpha*fProfileDelta;
877 else
878 {
879 argument = fProfileDelta*dTheta;
880 result = CLHEP::pi*argument*G4Exp(fProfileAlpha*argument);
881 result /= std::sinh(CLHEP::pi*argument);
882 result -= 1.;
883 result /= dTheta;
884 }
885 return result;
886 }
887
888 /////////////////////////////////////////////////////////////////
889 //
890 //
891
892 inline G4double G4NuclNuclDiffuseElastic::ProfileFar(G4double theta)
893 {
894 G4double dTheta = fRutherfordTheta + theta;
895 G4double argument = fProfileDelta*dTheta;
896
897 G4double result = CLHEP::pi*argument*G4Exp(fProfileAlpha*argument);
898 result /= std::sinh(CLHEP::pi*argument);
899 result /= dTheta;
900
901 return result;
902 }
903
904 /////////////////////////////////////////////////////////////////
905 //
906 //
907
908 inline G4double G4NuclNuclDiffuseElastic::Profile(G4double theta)
909 {
910 G4double dTheta = fRutherfordTheta - theta;
911 G4double result = 0., argument = 0.;
912
913 if(std::abs(dTheta) < 0.001) result = 1.;
914 else
915 {
916 argument = fProfileDelta*dTheta;
917 result = CLHEP::pi*argument;
918 result /= std::sinh(result);
919 }
920 return result;
921 }
922
923 /////////////////////////////////////////////////////////////////
924 //
925 //
926
927 inline G4complex G4NuclNuclDiffuseElastic::PhaseNear(G4double theta)
928 {
929 G4double twosigma = 2.*fCoulombPhase0;
930 twosigma -= fZommerfeld*G4Log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
931 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
932 twosigma -= fProfileLambda*theta - 0.25*CLHEP::pi;
933
934 twosigma *= fCofPhase;
935
936 G4complex z = G4complex(0., twosigma);
937
938 return std::exp(z);
939 }
940
941 /////////////////////////////////////////////////////////////////
942 //
943 //
944
945 inline G4complex G4NuclNuclDiffuseElastic::PhaseFar(G4double theta)
946 {
947 G4double twosigma = 2.*fCoulombPhase0;
948 twosigma -= fZommerfeld*G4Log(fHalfRutThetaTg2/(1.+fHalfRutThetaTg2));
949 twosigma += fRutherfordTheta*fZommerfeld/fHalfRutThetaTg - CLHEP::halfpi;
950 twosigma += fProfileLambda*theta - 0.25*CLHEP::pi;
951
952 twosigma *= fCofPhase;
953
954 G4complex z = G4complex(0., twosigma);
955
956 return std::exp(z);
957 }
958
959
960 /////////////////////////////////////////////////////////////////
961 //
962 //
963
964 inline G4complex G4NuclNuclDiffuseElastic::AmplitudeFar(G4double theta)
965 {
966 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
967 G4complex out = G4complex(kappa/fWaveVector,0.);
968 out *= ProfileFar(theta);
969 out *= PhaseFar(theta);
970 return out;
971 }
972
973
974 /////////////////////////////////////////////////////////////////
975 //
976 //
977
978 inline G4complex G4NuclNuclDiffuseElastic::Amplitude(G4double theta)
979 {
980
981 G4complex out = AmplitudeNear(theta) + fCofFar*AmplitudeFar(theta);
982 // G4complex out = AmplitudeNear(theta);
983 // G4complex out = AmplitudeFar(theta);
984 return out;
985 }
986
987 /////////////////////////////////////////////////////////////////
988 //
989 //
990
991 inline G4double G4NuclNuclDiffuseElastic::AmplitudeMod2(G4double theta)
992 {
993 G4complex out = Amplitude(theta);
994 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
995 return mod2;
996 }
997
998
999 /////////////////////////////////////////////////////////////////
1000 //
1001 //
1002
1003 inline G4double G4NuclNuclDiffuseElastic::GetRatioSim(G4double theta)
1004 {
1005 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1006 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1007 G4double sindTheta = std::sin(dTheta);
1008
1009 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1010 // G4cout<<"order = "<<order<<G4endl;
1011 G4double cosFresnel = 0.5 - GetCint(order);
1012 G4double sinFresnel = 0.5 - GetSint(order);
1013
1014 G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1015
1016 return out;
1017 }
1018
1019 /////////////////////////////////////////////////////////////////
1020 //
1021 // The xsc for Fresnel smooth nucleus profile
1022
1023 inline G4double G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc(G4double theta)
1024 {
1025 G4double ratio = GetRatioGen(theta);
1026 G4double ruthXsc = GetRutherfordXsc(theta);
1027 G4double xsc = ratio*ruthXsc;
1028 return xsc;
1029 }
1030
1031 /////////////////////////////////////////////////////////////////
1032 //
1033 // The xsc for Fresnel smooth nucleus profile for integration
1034
1035 inline G4double G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc(G4double alpha)
1036 {
1037 G4double theta = std::sqrt(alpha);
1038 G4double xsc = GetFresnelDiffuseXsc(theta);
1039 return xsc;
1040 }
1041
1042 /////////////////////////////////////////////////////////////////
1043 //
1044 //
1045
1046 inline G4double G4NuclNuclDiffuseElastic::AmplitudeSimMod2(G4double theta)
1047 {
1048 G4complex out = AmplitudeSim(theta);
1049 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1050 return mod2;
1051 }
1052
1053
1054 /////////////////////////////////////////////////////////////////
1055 //
1056 //
1057
1058 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGlaMod2(G4double theta)
1059 {
1060 G4complex out = AmplitudeGla(theta);
1061 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1062 return mod2;
1063 }
1064
1065 /////////////////////////////////////////////////////////////////
1066 //
1067 //
1068
1069 inline G4double G4NuclNuclDiffuseElastic::AmplitudeGGMod2(G4double theta)
1070 {
1071 G4complex out = AmplitudeGG(theta);
1072 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1073 return mod2;
1074 }
1075
1076 ////////////////////////////////////////////////////////////////////////////////////
1077 //
1078 //
1079
1080 inline G4double G4NuclNuclDiffuseElastic::CalcMandelstamS( const G4double mp ,
1081 const G4double mt ,
1082 const G4double Plab )
1083 {
1084 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1085 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1086
1087 return sMand;
1088 }
1089
1090
1091
1092 #endif
1093