ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclNuclDiffuseElastic.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclNuclDiffuseElastic.hh
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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