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