137 for ( std::vector<G4PhysicsTable*>::iterator it =
fAngleBank.begin();
139 if ( (*it) ) (*it)->clearAndDestroy();
163 for(jEl = 0 ; jEl < numOfEl; ++jEl)
173 G4cout<<
"G4NuclNuclDiffuseElastic::Initialise() the element: "
174 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
223 if (iZ == 1 && iA == 1) theDef =
theProton;
226 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
227 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
243 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0;
246 G4double thetaCMS = std::acos(cost);
278 if( z && (kRt > kRtC) )
310 if (iZ == 1 && iA == 1) theDef =
theProton;
313 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
314 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
330 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0;
333 G4double thetaCMS = std::acos(cost);
360 if (iZ == 1 && iA == 1) theDef =
theProton;
363 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
364 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
380 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0;
383 G4double thetaCMS = std::acos(cost);
403 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
416 bzero2 = bzero*bzero;
420 bonebyarg2 = bonebyarg*bonebyarg;
435 diffuse = 0.63*
fermi;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2;
483 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
496 bzero2 = bzero*bzero;
500 bonebyarg2 = bonebyarg*bonebyarg;
504 diffuse = 0.63*
fermi;
513 diffuse = 0.63*
fermi;
527 G4double sinHalfTheta = std::sin(0.5*theta);
528 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
554 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone;
558 sigma += kr2*bonebyarg2;
576 theta = std::sqrt(
alpha);
580 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
593 bzero2 = bzero*bzero;
597 bonebyarg2 = bonebyarg*bonebyarg;
601 diffuse = 0.63*
fermi;
610 diffuse = 0.63*
fermi;
625 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
651 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone;
655 sigma += kr2*bonebyarg2;
713 G4double t = 2*p*p*( 1 - std::cos(theta) );
727 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
737 if (thetaMax >
pi) thetaMax =
pi;
746 for(i = 1; i <= iMax; i++)
748 theta1 = (i-1)*thetaMax/iMax;
749 theta2 = i*thetaMax/iMax;
754 result = 0.5*(theta1 + theta2);
758 if (i > iMax ) result = 0.5*(theta1 + theta2);
764 if(result < 0.) result = 0.;
765 if(result > thetaMax) result = thetaMax;
785 G4double totElab = std::sqrt(m1*m1+p*p);
859 G4int iMomentum, iAngle;
881 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
883 for( iMomentum = 0; iMomentum <
fEnergyBin; iMomentum++)
887 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
893 if ( iMomentum < 0 ) iMomentum = 0;
896 if (iMomentum ==
fEnergyBin -1 || iMomentum == 0 )
902 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
921 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
944 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
963 randAngle = W1*theta1 + W2*theta2;
990 G4cout<<
"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
1011 G4double alpha1,
alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1026 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1032 if(alphaMax >
pi) alphaMax =
pi;
1061 alpha1 = alphaCoulomb + delth*(j-1);
1070 angleVector->
PutValue( j-1 , alpha1, sum );
1088 G4double x1, x2, y1, y2, randAngle;
1092 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1099 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1101 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1102 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1104 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1105 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1107 if ( x1 == x2 ) randAngle = x2;
1110 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*
G4UniformRand();
1113 randAngle = x1 + (
position - y1 )*( x2 - x1 )/( y2 - y1 );
1152 t =
SampleT( theParticle, ptot,
A);
1155 if(!(t < 0.0 || t >= 0.0))
1159 G4cout <<
"G4NuclNuclDiffuseElastic:WARNING: A = " <<
A
1160 <<
" mom(GeV)= " << plab/
GeV
1161 <<
" S-wave will be sampled"
1168 G4cout <<
" t= " << t <<
" tmax= " << tmax
1169 <<
" ptot= " << ptot <<
G4endl;
1182 else if( cost <= -1.0)
1189 sint = std::sqrt((1.0-cost)*(1.0+cost));
1193 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1195 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1236 G4double cost = std::cos(thetaCMS);
1244 else if( cost <= -1.0)
1251 sint = std::sqrt((1.0-cost)*(1.0+cost));
1255 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1257 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1297 G4double cost = std::cos(thetaLab);
1305 else if( cost <= -1.0)
1312 sint = std::sqrt((1.0-cost)*(1.0+cost));
1316 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1318 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1346 G4cout<<
"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1357 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1358 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1372 alphaMax = kRmax*kRmax/kR2;
1374 if (alphaMax > 4.) alphaMax = 4.;
1376 alphaCoulomb = kRcoul*kRcoul/kR2;
1381 fBeta = a/std::sqrt(1+a*a);
1415 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1417 angleVector->
PutValue( j-1 , alpha1, sumL10 );
1445 if (
n < 0 ) legPol = 0.;
1446 else if(
n == 0 ) legPol = 1.;
1447 else if(
n == 1 ) legPol = x;
1448 else if(
n == 2 ) legPol = (3.*x*x-1.)/2.;
1449 else if(
n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
1450 else if(
n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1451 else if(
n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
1452 else if(
n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1469 G4double n2, cofn, shny, chny, fn, gn;
1488 for(
n = 1;
n <= nMax;
n++)
1492 cofn =
G4Exp(-0.5*n2)/(n2+twox2);
1494 chny = std::cosh(
n*y);
1495 shny = std::sinh(
n*y);
1497 fn = twox - twoxcos2xy*chny +
n*sin2xy*shny;
1498 gn = twoxsin2xy*chny +
n*cos2xy*shny;
1509 if(std::abs(x) < 0.0001)
1516 outRe +=
GetErf(x) + cof1*(1-cos2xy)/twox;
1517 outIm += cof1*sin2xy/twox;
1567 order /= std::sqrt(2.);
1570 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1571 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1595 order /= std::sqrt(2.);
1597 G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1598 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1637 G4double sindTheta = std::sin(dTheta);
1638 G4double persqrt2 = std::sqrt(0.5);
1678 shiftN = std::exp( -0.5*(1.-im*
fEtaRatio)*T12b ) - 1.;
1694 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1695 G4double sinThetaH2 = sinThetaH*sinThetaH;
1746 fBeta = a/std::sqrt(1+a*a);
1781 fBeta = a/std::sqrt(1+a*a);
1821 if( pN < 0. ) pN = 0.;
1824 if( tN < 0. ) tN = 0.;
1843 fBeta = a/std::sqrt(1+a*a);
1874 G4double proj_energy = proj_mass + pTkin;
1875 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1890 if( proj_momentum >= 1.2 )
1894 else if( proj_momentum >= 0.6 )
1901 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1907 if( proj_momentum >= 10. )
1917 A0 = 100. - B0*
G4Log(3.0e7);
1919 xsection = A0 + B0*
G4Log(proj_energy) - 11
1921 0.93827*0.93827,-0.165);
1926 if(pParticle == tParticle)
1928 if( proj_momentum < 0.73 )
1932 else if( proj_momentum < 1.05 )
1934 hnXsc = 23 + 40*(
G4Log(proj_momentum/0.73))*
1935 (
G4Log(proj_momentum/0.73));
1946 if( proj_momentum < 0.8 )
1950 else if( proj_momentum < 1.4 )
1976 G4double sindTheta = std::sin(dTheta);
1983 order = std::abs(order);
1993 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1994 out += ( cosFresnel + sinFresnel - 1. )*prof;
1998 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
2011 const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
2012 24.01409824083091, -1.231739572450155,
2013 0.1208650973866179e-2, -0.5395239384953e-5 } ;
2017 tmp -= (z + 0.5) * std::log(tmp);
2020 for ( j = 0; j <= 5; j++ )
2025 return -tmp + std::log(2.5066282746310005*ser);
2035 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2037 modvalue = std::fabs(value);
2039 if ( value < 8.0 && value > -8.0 )
2041 value2 = value*value;
2043 fact1 = 57568490574.0 + value2*(-13362590354.0
2044 + value2*( 651619640.7
2045 + value2*(-11214424.18
2046 + value2*( 77392.33017
2047 + value2*(-184.9052456 ) ) ) ) );
2049 fact2 = 57568490411.0 + value2*( 1029532985.0
2050 + value2*( 9494680.718
2051 + value2*(59272.64853
2052 + value2*(267.8532712
2053 + value2*1.0 ) ) ) );
2055 bessel = fact1/fact2;
2063 shift = modvalue-0.785398164;
2065 fact1 = 1.0 + value2*(-0.1098628627e-2
2066 + value2*(0.2734510407e-4
2067 + value2*(-0.2073370639e-5
2068 + value2*0.2093887211e-6 ) ) );
2070 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2071 + value2*(-0.6911147651e-5
2072 + value2*(0.7621095161e-6
2073 - value2*0.934945152e-7 ) ) );
2075 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2087 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2089 modvalue = std::fabs(value);
2091 if ( modvalue < 8.0 )
2093 value2 = value*value;
2095 fact1 = value*(72362614232.0 + value2*(-7895059235.0
2096 + value2*( 242396853.1
2097 + value2*(-2972611.439
2098 + value2*( 15704.48260
2099 + value2*(-30.16036606 ) ) ) ) ) );
2101 fact2 = 144725228442.0 + value2*(2300535178.0
2102 + value2*(18583304.74
2103 + value2*(99447.43394
2104 + value2*(376.9991397
2105 + value2*1.0 ) ) ) );
2106 bessel = fact1/fact2;
2114 shift = modvalue - 2.356194491;
2116 fact1 = 1.0 + value2*( 0.183105e-2
2117 + value2*(-0.3516396496e-4
2118 + value2*(0.2457520174e-5
2119 + value2*(-0.240337019e-6 ) ) ) );
2121 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2122 + value2*( 0.8449199096e-5
2123 + value2*(-0.88228987e-6
2124 + value2*0.105787412e-6 ) ) );
2126 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2128 if (value < 0.0) bessel = -bessel;
static const G4double e1[44]
static const G4double e2[44]
G4double epsilon(G4double density, G4double temperature)
std::vector< G4Element * > G4ElementTable
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static const G4double alpha
G4double G4Log(G4double x)
static constexpr double twopi
static constexpr double keV
static constexpr double fermi
static constexpr double GeV
static constexpr double MeV
static constexpr double degree
static constexpr double pi
std::complex< G4double > G4complex
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
static G4Deuteron * Deuteron()
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
G4double GetKineticEnergy() const
G4double GetTotalMomentum() const
static G4ElementTable * GetElementTable()
static size_t GetNumberOfElements()
G4double GetTotalMomentum() const
const G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & Get4Momentum() const
void SetMinEnergy(G4double anEnergy)
void SetMaxEnergy(const G4double anEnergy)
G4double GetMaxEnergy() const
static G4HadronicParameters * Instance()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Neutron * Neutron()
static G4NistManager * Instance()
G4double GetAtomicMassAmu(const G4String &symb) const
std::vector< G4String > fElementNameVector
std::vector< G4double > fElementNumberVector
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double fRutherfordRatio
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
const G4ParticleDefinition * thePionPlus
G4double GetRatioGen(G4double theta)
G4double Profile(G4double theta)
G4complex AmplitudeSim(G4double theta)
const G4ParticleDefinition * thePionMinus
G4double GetDiffElasticSumProb(G4double theta)
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
void InitialiseOnFly(G4double Z, G4double A)
G4NuclNuclDiffuseElastic()
G4double BesselJzero(G4double z)
G4ParticleDefinition * theProton
G4double fCofAlphaCoulomb
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4complex AmplitudeGla(G4double theta)
virtual ~G4NuclNuclDiffuseElastic()
G4complex GetErfComp(G4complex z, G4int nMax)
G4ParticleDefinition * theAlpha
G4double CalculateCoulombPhase(G4int n)
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetLegendrePol(G4int n, G4double x)
G4double GetIntegrandFunction(G4double theta)
G4complex AmplitudeNear(G4double theta)
G4complex AmplitudeGG(G4double theta)
G4double GetSint(G4double x)
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double lowEnergyRecoilLimit
G4double ProfileNear(G4double theta)
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4complex GetErfInt(G4complex z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
void CalculateRutherfordAnglePar()
G4double GetCint(G4double x)
void InitParametersGla(const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
G4double BesselOneByArg(G4double z)
G4double fNuclearRadiusCof
G4complex PhaseNear(G4double theta)
G4complex GammaLogarithm(G4complex xx)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4ParticleDefinition * theNeutron
G4complex CoulombAmplitude(G4double theta)
G4double lowestEnergyLimit
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double GetDiffElasticSumProbA(G4double alpha)
G4double lowEnergyLimitHE
void CalculateCoulombPhaseZero()
G4ParticleDefinition * theDeuteron
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double fNuclearRadiusSquare
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
std::vector< G4PhysicsTable * > fAngleBank
G4double GetExpSin(G4double x)
G4double fHalfRutThetaTg2
G4double BesselJone(G4double z)
G4PhysicsTable * fAngleTable
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
G4double CalculateNuclearRad(G4double A)
void InitParameters(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
G4PhysicsLogVector * fEnergyVector
G4complex GetErfcInt(G4complex z)
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double fRutherfordTheta
G4double GetFresnelIntegrandXsc(G4double alpha)
G4double GetExpCos(G4double x)
G4double GetErf(G4double x)
G4double SampleCoulombMuCMS(const G4ParticleDefinition *aParticle, G4double p)
G4double DampFactor(G4double z)
const G4ParticleDefinition * fParticle
G4complex GammaMore(G4double theta)
G4double GetDiffElasticProb(G4double theta)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4complex GammaLess(G4double theta)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetPDGMass() const
G4double GetPDGCharge() const
G4int GetBaryonNumber() const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
void PutValue(const std::size_t index, const G4double e, const G4double value)
void insertAt(std::size_t, G4PhysicsVector *)
G4double GetLowEdgeEnergy(const std::size_t index) const
static G4PionMinus * PionMinus()
static G4PionPlus * PionPlus()
static G4Pow * GetInstance()
G4double powA(G4double A, G4double y) const
static G4Proton * Proton()
static G4Triton * Triton()
static constexpr double millibarn
static constexpr double GeV
static constexpr double twopi
static constexpr double hbarc
static constexpr double pi
ThreeVector shoot(const G4int Ap, const G4int Af)