117 else if (
Z == 2 &&
A == 4) theTargetDef =
theAlpha;
136 fTmax = 4.0*ptot*ptot;
146 else { PrNucleonMass = NucleonMass; }
147 G4double energyPerN = std::sqrt(
sqr(PlabPerN) +
sqr(PrNucleonMass));
148 energyPerN -= PrNucleonMass;
160 const G4double mevToBarn = 0.38938e+6;
168 G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHadronic);
176 G4double par1 = 2.*(1.+Am)/(1.+ctet1);
180 G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
184 G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
185 G4double PtX= PtProjCMS * std::cos(phi);
186 G4double PtY= PtProjCMS * std::sin(phi);
189 Fproj.
setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
190 T = -(Pproj-Fproj).mag2();
208 if((
A>=12.) && (
A<27) )
fRa=
fRa*0.85;
209 if((
A>=27.) && (
A<48) )
fRa=
fRa*0.90;
210 if((
A>=48.) && (
A<65) )
fRa=
fRa*0.95;
222 { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
223 13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
224 if((Plab < 5500.)&&(Plab >= 610.) )
226 if((Plab >= 5500.)&&(Plab < 12300.) )
229 { rho = 0.135-2.26/(std::sqrt(
S)) ;}
230 Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(
S-4.*0.88))+0.04*
G4Log(
S) ;
231 ceff2 = 0.375 - 2./
S + 0.44/(
sqr(
S-4.)+1.5) ;
238 Ref2 =
fRa*
fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
239 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*
G4Exp(-0.03*sig_pbarp);
243 Ref2 =
fRa*
fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
244 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
248 Ref2 =
fRa*
fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
249 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
253 Ref2 =
fRa*
fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
254 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*
G4Exp(-0.03*sig_pbarp);
258 Ref2 =
fRa*
fRa +2.48*0.01*sig_pbarp*
fRa - 2.23e-6*sig_pbarp*sig_pbarp*
fRa*
fRa;
259 ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*
G4Exp(-0.03*sig_pbarp);
267 ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*
G4Exp(-0.03*sig_pbarp);
271 ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 *
G4Exp(-0.03*sig_pbarp);
275 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 *
G4Exp(-0.02*sig_pbarp);
279 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 *
G4Exp(-0.02*sig_pbarp);
283 ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 *
G4Exp(-0.03*sig_pbarp);
291 ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*
G4Exp(-0.03*sig_pbarp);
295 ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 *
G4Exp(-0.02*sig_pbarp);
299 ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 *
G4Exp(-0.02*sig_pbarp);
303 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 *
G4Exp(-0.03*sig_pbarp);
307 ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*
G4Exp(-0.03*sig_pbarp);
315 ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*
G4Exp(-0.03*sig_pbarp);
319 ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 *
G4Exp(-0.02*sig_pbarp);
323 ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 *
G4Exp(-0.03*sig_pbarp);
327 ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 *
G4Exp(-0.03*sig_pbarp);
331 ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 *
G4Exp(-0.03*sig_pbarp);
335 fRef=std::sqrt(Ref2);
336 fceff = std::sqrt(ceff2);
341 const G4int maxNumberOfLoops = 10000;
342 G4int loopCounter = 0;
349 BracFunct = BracFunct *
Q;
352 ++loopCounter < maxNumberOfLoops );
353 if ( loopCounter >= maxNumberOfLoops ) {
375 if(!(T < 0.0 || T >= 0.0))
379 G4cout <<
"G4DiffuseElastic:WARNING: A = " <<
A
380 <<
" mom(GeV)= " << plab/
GeV
381 <<
" S-wave will be sampled"
391 if(cosTet > 1.0 ) cosTet= 1.;
392 if(cosTet < -1.0 ) cosTet=-1.;
411 if(!(T < 0.0 || T >= 0.0))
415 G4cout <<
"G4DiffuseElastic:WARNING: A = " <<
A
416 <<
" mom(GeV)= " << plab/
GeV
417 <<
" S-wave will be sampled"
434 else if( cost <= -1.0)
441 sint = std::sqrt((1.0-cost)*(1.0+cost));
465 if( std::fabs(x) < 0.01 )
485 fBeta = a/std::sqrt(1+a*a);
521 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
523 modvalue = std::fabs(value);
525 if ( value < 8.0 && value > -8.0 )
527 value2 = value*value;
529 fact1 = 57568490574.0 + value2*(-13362590354.0
530 + value2*( 651619640.7
531 + value2*(-11214424.18
532 + value2*( 77392.33017
533 + value2*(-184.9052456 ) ) ) ) );
535 fact2 = 57568490411.0 + value2*( 1029532985.0
536 + value2*( 9494680.718
537 + value2*(59272.64853
538 + value2*(267.8532712
539 + value2*1.0 ) ) ) );
541 bessel = fact1/fact2;
549 shift = modvalue-0.785398164;
551 fact1 = 1.0 + value2*(-0.1098628627e-2
552 + value2*(0.2734510407e-4
553 + value2*(-0.2073370639e-5
554 + value2*0.2093887211e-6 ) ) );
555 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
556 + value2*(-0.6911147651e-5
557 + value2*(0.7621095161e-6
558 - value2*0.934945152e-7 ) ) );
560 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
572 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
574 modvalue = std::fabs(value);
576 if ( modvalue < 8.0 )
578 value2 = value*value;
579 fact1 = value*(72362614232.0 + value2*(-7895059235.0
580 + value2*( 242396853.1
581 + value2*(-2972611.439
582 + value2*( 15704.48260
583 + value2*(-30.16036606 ) ) ) ) ) );
585 fact2 = 144725228442.0 + value2*(2300535178.0
586 + value2*(18583304.74
587 + value2*(99447.43394
588 + value2*(376.9991397
589 + value2*1.0 ) ) ) );
590 bessel = fact1/fact2;
597 shift = modvalue - 2.356194491;
599 fact1 = 1.0 + value2*( 0.183105e-2
600 + value2*(-0.3516396496e-4
601 + value2*(0.2457520174e-5
602 + value2*(-0.240337019e-6 ) ) ) );
604 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
605 + value2*( 0.8449199096e-5
606 + value2*(-0.88228987e-6
607 + value2*0.105787412e-6 ) ) );
609 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
610 if (value < 0.0) bessel = -bessel;
621 if( std::fabs(x) < 0.01 )
625 result = (2.- x2 + x2*x2/6.)/4.;
644 if(cteta1 < -1.) cteta1 = -1.0;
G4double S(G4double temp)
static const G4double reg
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
static constexpr double twopi
static constexpr double millibarn
static constexpr double fermi
static constexpr double GeV
static constexpr double MeV
static constexpr double pi
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
static G4AntiAlpha * AntiAlpha()
static G4AntiDeuteron * AntiDeuteron()
static G4AntiHe3 * AntiHe3()
static G4AntiNeutron * AntiNeutron()
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) override
~G4AntiNuclElastic() override
G4ParticleDefinition * theAProton
G4double BesselOneByArg(G4double z)
G4ParticleDefinition * theAAlpha
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4ParticleDefinition * theNeutron
G4ComponentAntiNuclNuclearXS * cs
G4ParticleDefinition * theProton
G4double SampleThetaLab(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4ParticleDefinition * theADeuteron
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4ParticleDefinition * theANeutron
G4ParticleDefinition * theATriton
G4ParticleDefinition * theAHe3
G4double DampFactor(G4double z)
G4double BesselJone(G4double z)
G4double GetcosTeta1(G4double plab, G4int A)
G4ParticleDefinition * theDeuteron
G4double BesselJzero(G4double z)
G4double SampleThetaCMS(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4ParticleDefinition * theAlpha
const G4ParticleDefinition * fParticle
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
static G4AntiProton * AntiProton()
static G4AntiTriton * AntiTriton()
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double A)
G4double GetAntiHadronNucleonTotCrSc(const G4ParticleDefinition *aParticle, G4double kinEnergy)
static G4CrossSectionDataSetRegistry * Instance()
static G4Deuteron * Deuteron()
static G4Neutron * Neutron()
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetPDGMass() const
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
static G4Pow * GetInstance()
G4double A13(G4double A) const
G4double Z13(G4int Z) const
G4double Z23(G4int Z) const
static G4Proton * Proton()
static G4Triton * Triton()
G4double energy(const ThreeVector &p, const G4double m)
static const G4double Z1[5]