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