Loading...
Searching...
No Matches
G4NuclNuclDiffuseElastic Class Reference
#include <G4NuclNuclDiffuseElastic.hh>
Inheritance diagram for G4NuclNuclDiffuseElastic:
Additional Inherited Members | |
| Protected Member Functions inherited from G4HadronicInteraction | |
| void | SetModelName (const G4String &nam) |
| G4bool | IsBlocked () const |
| void | Block () |
| Protected Attributes inherited from G4HadronElastic | |
| G4double | pLocalTmax |
| Protected Attributes inherited from G4HadronicInteraction | |
| G4HadFinalState | theParticleChange |
| G4int | verboseLevel |
| G4double | theMinEnergy |
| G4double | theMaxEnergy |
| G4bool | isBlocked |
Detailed Description
Definition at line 61 of file G4NuclNuclDiffuseElastic.hh.
Constructor & Destructor Documentation
◆ G4NuclNuclDiffuseElastic()
| G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic | ( | ) |
Definition at line 68 of file G4NuclNuclDiffuseElastic.cc.
70{
71 SetMinEnergy( 50*MeV );
73 verboseLevel = 0;
74 lowEnergyRecoilLimit = 100.*keV;
75 lowEnergyLimitQ = 0.0*GeV;
76 lowEnergyLimitHE = 0.0*GeV;
77 lowestEnergyLimit= 0.0*keV;
78 plabLowLimit = 20.0*MeV;
79
80 theProton = G4Proton::Proton();
81 theNeutron = G4Neutron::Neutron();
82 theDeuteron = G4Deuteron::Deuteron();
83 theAlpha = G4Alpha::Alpha();
84 thePionPlus = G4PionPlus::PionPlus();
85 thePionMinus= G4PionMinus::PionMinus();
86
87 fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV
88 fAngleBin = 200;
89
91 fAngleTable = 0;
92
93 fParticle = 0;
94 fWaveVector = 0.;
95 fAtomicWeight = 0.;
96 fAtomicNumber = 0.;
97 fNuclearRadius = 0.;
98 fBeta = 0.;
99 fZommerfeld = 0.;
100 fAm = 0.;
101 fAddCoulomb = false;
102 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
103
104 // Empirical parameters
105
106 fCofAlphaMax = 1.5;
107 fCofAlphaCoulomb = 0.5;
108
109 fProfileDelta = 1.;
110 fProfileAlpha = 0.5;
111
112 fCofLambda = 1.0;
113 fCofDelta = 0.04;
114 fCofAlpha = 0.095;
115
116 fNuclearRadius1 = fNuclearRadius2 = fNuclearRadiusSquare
117 = fRutherfordRatio = fCoulombPhase0 = fHalfRutThetaTg = fHalfRutThetaTg2
118 = fRutherfordTheta = fProfileLambda = fCofPhase = fCofFar
119 = fSumSigma = fEtaRatio = fReZ = 0.0;
120 fMaxL = 0;
121
122 fNuclearRadiusCof = 1.0;
123 fCoulombMuC = 0.0;
124}
Definition: G4HadronElastic.hh:49
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:89
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:102
G4double GetMaxEnergy() const
Definition: G4HadronicInteraction.hh:96
static G4HadronicParameters * Instance()
Definition: G4HadronicParameters.cc:52
Definition: G4PhysicsLogVector.hh:48
Referenced by BuildAngleTable(), GetCint(), GetErfInt(), GetSint(), IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().
◆ ~G4NuclNuclDiffuseElastic()
|
virtual |
Definition at line 130 of file G4NuclNuclDiffuseElastic.cc.
131{
132 if ( fEnergyVector ) {
133 delete fEnergyVector;
134 fEnergyVector = 0;
135 }
136
137 for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
138 it != fAngleBank.end(); ++it ) {
139 if ( (*it) ) (*it)->clearAndDestroy();
140 delete *it;
141 *it = 0;
142 }
143 fAngleTable = 0;
144}
Member Function Documentation
◆ Amplitude()
Definition at line 978 of file G4NuclNuclDiffuseElastic.hh.
979{
980
982 // G4complex out = AmplitudeNear(theta);
983 // G4complex out = AmplitudeFar(theta);
984 return out;
985}
G4complex AmplitudeFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:964
G4complex AmplitudeNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1608
Referenced by AmplitudeMod2().
◆ AmplitudeFar()
Definition at line 964 of file G4NuclNuclDiffuseElastic.hh.
965{
966 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
968 out *= ProfileFar(theta);
969 out *= PhaseFar(theta);
970 return out;
971}
G4complex PhaseFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:945
G4double ProfileFar(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:892
Referenced by Amplitude().
◆ AmplitudeGG()
Definition at line 1691 of file G4NuclNuclDiffuseElastic.cc.
1692{
1694 G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1695 G4double sinThetaH2 = sinThetaH*sinThetaH;
1698
1699 a = -fSumSigma/CLHEP::twopi/fNuclearRadiusSquare;
1700 b2 = fWaveVector*fWaveVector*fNuclearRadiusSquare*sinThetaH2;
1701
1702 aTemp = a;
1703
1705 {
1707 aTemp *= a;
1708 out += T12b;
1710 }
1711 out *= -4.*im*fWaveVector/CLHEP::pi;
1712 out += CoulombAmplitude(theta);
1713 return out;
1714}
G4GLOB_DLL std::ostream G4cout
G4complex CoulombAmplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:797
Referenced by AmplitudeGGMod2().
◆ AmplitudeGGMod2()
Definition at line 1069 of file G4NuclNuclDiffuseElastic.hh.
1070{
1072 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1073 return mod2;
1074}
G4complex AmplitudeGG(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1691
◆ AmplitudeGla()
Definition at line 1664 of file G4NuclNuclDiffuseElastic.cc.
1665{
1670
1672 {
1673 shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1674 // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1675 b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1676 b2 = b*b;
1677 T12b = fSumSigma*G4Exp(-b2/fNuclearRadiusSquare)/CLHEP::pi/fNuclearRadiusSquare;
1678 shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1680 }
1681 out /= 2.*im*fWaveVector;
1682 out += CoulombAmplitude(theta);
1683 return out;
1684}
G4double CalculateCoulombPhase(G4int n)
Definition: G4NuclNuclDiffuseElastic.hh:844
G4double GetLegendrePol(G4int n, G4double x)
Definition: G4NuclNuclDiffuseElastic.cc:1440
Referenced by AmplitudeGlaMod2().
◆ AmplitudeGlaMod2()
Definition at line 1058 of file G4NuclNuclDiffuseElastic.hh.
1059{
1061 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1062 return mod2;
1063}
G4complex AmplitudeGla(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1664
◆ AmplitudeMod2()
Definition at line 991 of file G4NuclNuclDiffuseElastic.hh.
992{
994 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
995 return mod2;
996}
G4complex Amplitude(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:978
◆ AmplitudeNear()
Definition at line 1608 of file G4NuclNuclDiffuseElastic.cc.
1609{
1610 G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1612
1613 out *= PhaseNear(theta);
1614
1615 if( theta <= fRutherfordTheta )
1616 {
1618 // out *= GammaMore(theta) + ProfileNear(theta);
1619 out += CoulombAmplitude(theta);
1620 }
1621 else
1622 {
1624 // out *= GammaLess(theta) + ProfileNear(theta);
1625 }
1626 return out;
1627}
G4double ProfileNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:871
G4complex PhaseNear(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:927
G4complex GammaMore(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1581
G4complex GammaLess(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1553
Referenced by Amplitude().
◆ AmplitudeSim()
Definition at line 1633 of file G4NuclNuclDiffuseElastic.cc.
1634{
1635 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1636 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1637 G4double sindTheta = std::sin(dTheta);
1638 G4double persqrt2 = std::sqrt(0.5);
1639
1641 order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1642 // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1643
1644 G4complex out;
1645
1646 if ( theta <= fRutherfordTheta )
1647 {
1649 }
1650 else
1651 {
1653 }
1654
1655 out *= CoulombAmplitude(theta);
1656
1657 return out;
1658}
G4complex GetErfcInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.hh:672
Referenced by AmplitudeSimMod2().
◆ AmplitudeSimMod2()
Definition at line 1046 of file G4NuclNuclDiffuseElastic.hh.
1047{
1049 G4double mod2 = out.real()*out.real() + out.imag()*out.imag();
1050 return mod2;
1051}
G4complex AmplitudeSim(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1633
◆ BesselJone()
Definition at line 2085 of file G4NuclNuclDiffuseElastic.cc.
2086{
2087 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2088
2089 modvalue = std::fabs(value);
2090
2091 if ( modvalue < 8.0 )
2092 {
2093 value2 = value*value;
2094
2095 fact1 = value*(72362614232.0 + value2*(-7895059235.0
2096 + value2*( 242396853.1
2097 + value2*(-2972611.439
2098 + value2*( 15704.48260
2099 + value2*(-30.16036606 ) ) ) ) ) );
2100
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;
2107 }
2108 else
2109 {
2110 arg = 8.0/modvalue;
2111
2112 value2 = arg*arg;
2113
2114 shift = modvalue - 2.356194491;
2115
2116 fact1 = 1.0 + value2*( 0.183105e-2
2117 + value2*(-0.3516396496e-4
2118 + value2*(0.2457520174e-5
2119 + value2*(-0.240337019e-6 ) ) ) );
2120
2121 fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2122 + value2*( 0.8449199096e-5
2123 + value2*(-0.88228987e-6
2124 + value2*0.105787412e-6 ) ) );
2125
2126 bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2127
2128 if (value < 0.0) bessel = -bessel;
2129 }
2130 return bessel;
2131}
Referenced by BesselOneByArg(), GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BesselJzero()
Definition at line 2033 of file G4NuclNuclDiffuseElastic.cc.
2034{
2035 G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2036
2037 modvalue = std::fabs(value);
2038
2039 if ( value < 8.0 && value > -8.0 )
2040 {
2041 value2 = value*value;
2042
2043 fact1 = 57568490574.0 + value2*(-13362590354.0
2044 + value2*( 651619640.7
2045 + value2*(-11214424.18
2046 + value2*( 77392.33017
2047 + value2*(-184.9052456 ) ) ) ) );
2048
2049 fact2 = 57568490411.0 + value2*( 1029532985.0
2050 + value2*( 9494680.718
2051 + value2*(59272.64853
2052 + value2*(267.8532712
2053 + value2*1.0 ) ) ) );
2054
2055 bessel = fact1/fact2;
2056 }
2057 else
2058 {
2059 arg = 8.0/modvalue;
2060
2061 value2 = arg*arg;
2062
2063 shift = modvalue-0.785398164;
2064
2065 fact1 = 1.0 + value2*(-0.1098628627e-2
2066 + value2*(0.2734510407e-4
2067 + value2*(-0.2073370639e-5
2068 + value2*0.2093887211e-6 ) ) );
2069
2070 fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2071 + value2*(-0.6911147651e-5
2072 + value2*(0.7621095161e-6
2073 - value2*0.934945152e-7 ) ) );
2074
2075 bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2076 }
2077 return bessel;
2078}
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BesselOneByArg()
Definition at line 419 of file G4NuclNuclDiffuseElastic.hh.
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}
G4double BesselJone(G4double z)
Definition: G4NuclNuclDiffuseElastic.cc:2085
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ BuildAngleTable()
| void G4NuclNuclDiffuseElastic::BuildAngleTable | ( | ) |
Definition at line 1007 of file G4NuclNuclDiffuseElastic.cc.
1008{
1009 G4int i, j;
1011 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1012
1013 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
1014
1015 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1016
1018
1019 for( i = 0; i < fEnergyBin; i++)
1020 {
1021 kinE = fEnergyVector->GetLowEdgeEnergy(i);
1022
1023 // G4cout<<G4endl;
1024 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
1025
1026 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1027
1028 InitDynParameters(fParticle, partMom);
1029
1030 alphaMax = fRutherfordTheta*fCofAlphaMax;
1031
1033
1034 // VI: Coverity complain
1035 //alphaMax = pi2;
1036
1037 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
1038
1039 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1040
1042
1043 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1044
1045 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1046
1047 sum = 0.;
1048
1049 // fAddCoulomb = false;
1050 fAddCoulomb = true;
1051
1052 // for(j = 1; j < fAngleBin; j++)
1053 for(j = fAngleBin-1; j >= 1; j--)
1054 {
1055 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1056 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1057
1058 // alpha1 = alphaMax*(j-1)/fAngleBin;
1059 // alpha2 = alphaMax*( j )/fAngleBin;
1060
1061 alpha1 = alphaCoulomb + delth*(j-1);
1062 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1063 alpha2 = alpha1 + delth;
1064
1065 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1066 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1067
1068 sum += delta;
1069
1071 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1072 }
1073 fAngleTable->insertAt(i,angleVector);
1074
1075 // delete[] angleVector;
1076 // delete[] angleBins;
1077 }
1078 return;
1079}
Definition: G4Integrator.hh:45
G4NuclNuclDiffuseElastic()
Definition: G4NuclNuclDiffuseElastic.cc:68
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
Definition: G4NuclNuclDiffuseElastic.cc:1767
G4double GetFresnelIntegrandXsc(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.hh:1035
Definition: G4PhysicsFreeVector.hh:52
void PutValue(std::size_t index, G4double energy, G4double dValue)
Definition: G4PhysicsFreeVector.hh:77
Definition: G4PhysicsTable.hh:56
G4double GetLowEdgeEnergy(std::size_t binNumber) const
Definition: G4PhysicsVector.cc:136
Referenced by Initialise(), and InitialiseOnFly().
◆ CalcMandelstamS()
|
inline |
Definition at line 1080 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by GetHadronNucleonXscNS().
◆ CalculateAm()
Definition at line 468 of file G4NuclNuclDiffuseElastic.hh.
469{
470 G4double k = momentum/CLHEP::hbarc;
473 G4double zn2 = zn*zn;
474 fAm = ch/zn2;
475
476 return fAm;
477}
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().
◆ CalculateCoulombPhase()
Definition at line 844 of file G4NuclNuclDiffuseElastic.hh.
845{
847 // G4complex gammalog = GammaLogarithm(z);
849 return gammalog.imag();
850}
G4complex GammaLogB2n(G4complex xx)
Definition: G4NuclNuclDiffuseElastic.hh:610
Referenced by AmplitudeGla().
◆ CalculateCoulombPhaseZero()
|
inline |
Definition at line 831 of file G4NuclNuclDiffuseElastic.hh.
832{
834 // G4complex gammalog = GammaLogarithm(z);
836 fCoulombPhase0 = gammalog.imag();
837}
Referenced by InitDynParameters(), InitParameters(), and InitParametersGla().
◆ CalculateNuclearRad()
Definition at line 483 of file G4NuclNuclDiffuseElastic.hh.
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 */
506
507 return radius;
508}
Referenced by GetDiffuseElasticSumXsc(), GetDiffuseElasticXsc(), Initialise(), InitialiseOnFly(), InitParameters(), InitParametersGla(), IntegralElasticProb(), SampleCoulombMuCMS(), SampleThetaCMS(), and TestAngleTable().
◆ CalculateParticleBeta()
|
inline |
Definition at line 441 of file G4NuclNuclDiffuseElastic.hh.
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), and GetDiffuseElasticSumXsc().
◆ CalculateRutherfordAnglePar()
|
inline |
Definition at line 858 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by InitDynParameters(), and InitParameters().
◆ CalculateZommerfeld()
|
inline |
Definition at line 456 of file G4NuclNuclDiffuseElastic.hh.
457{
458 fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
459
460 return fZommerfeld;
461}
Referenced by GetCoulombElasticXsc(), GetCoulombIntegralXsc(), GetCoulombTotalXsc(), GetDiffuseElasticSumXsc(), InitDynParameters(), InitParameters(), InitParametersGla(), and TestAngleTable().
◆ CoulombAmplitude()
Definition at line 797 of file G4NuclNuclDiffuseElastic.hh.
798{
799 G4complex ca;
800
801 G4double sinHalfTheta = std::sin(0.5*theta);
802 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
803 sinHalfTheta2 += fAm;
804
807 ca = std::exp(z);
808
809 ca *= -fZommerfeld/(2.*fWaveVector*sinHalfTheta2);
810
811 return ca;
812}
Referenced by AmplitudeGG(), AmplitudeGla(), AmplitudeNear(), AmplitudeSim(), and CoulombAmplitudeMod2().
◆ CoulombAmplitudeMod2()
Definition at line 818 of file G4NuclNuclDiffuseElastic.hh.
◆ DampFactor()
Definition at line 396 of file G4NuclNuclDiffuseElastic.hh.
397{
398 G4double df;
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}
Referenced by GetDiffElasticProb(), GetDiffElasticSumProb(), and GetDiffElasticSumProbA().
◆ GammaLess()
Definition at line 1553 of file G4NuclNuclDiffuseElastic.cc.
1554{
1555 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1556 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1557
1558 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1559 G4double kappa = u/std::sqrt(CLHEP::pi);
1560 G4double dTheta = theta - fRutherfordTheta;
1561 u *= dTheta;
1562 G4double u2 = u*u;
1563 G4double u2m2p3 = u2*2./3.;
1564
1567 order /= std::sqrt(2.);
1568
1571 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1573
1574 return out;
1575}
Referenced by AmplitudeNear().
◆ GammaLogarithm()
Definition at line 2009 of file G4NuclNuclDiffuseElastic.cc.
2010{
2012 24.01409824083091, -1.231739572450155,
2013 0.1208650973866179e-2, -0.5395239384953e-5 } ;
2014 G4int j;
2015 G4complex z = zz - 1.0;
2016 G4complex tmp = z + 5.5;
2017 tmp -= (z + 0.5) * std::log(tmp);
2019
2020 for ( j = 0; j <= 5; j++ )
2021 {
2022 z += 1.0;
2023 ser += cof[j]/z;
2024 }
2025 return -tmp + std::log(2.5066282746310005*ser);
2026}
◆ GammaLogB2n()
Definition at line 610 of file G4NuclNuclDiffuseElastic.hh.
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
623 result += 1./z1 - 1./z3 +1./z5 -1./z7;
624 return result;
625}
Referenced by CalculateCoulombPhase(), and CalculateCoulombPhaseZero().
◆ GammaMore()
Definition at line 1581 of file G4NuclNuclDiffuseElastic.cc.
1582{
1583 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1584 G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1585
1586 G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1587 G4double kappa = u/std::sqrt(CLHEP::pi);
1588 G4double dTheta = theta - fRutherfordTheta;
1589 u *= dTheta;
1590 G4double u2 = u*u;
1591 G4double u2m2p3 = u2*2./3.;
1592
1595 order /= std::sqrt(2.);
1598 G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1600
1601 return out;
1602}
Referenced by AmplitudeNear().
◆ GetCint()
Definition at line 767 of file G4NuclNuclDiffuseElastic.hh.
768{
769 G4double out;
770
771 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
772
774
775 return out;
776}
G4double GetCosHaPit2(G4double t)
Definition: G4NuclNuclDiffuseElastic.hh:195
Referenced by GetRatioGen(), and GetRatioSim().
◆ GetCofAlphaCoulomb()
|
inline |
Definition at line 293 of file G4NuclNuclDiffuseElastic.hh.
293{return fCofAlphaCoulomb;};
◆ GetCofAlphaMax()
|
inline |
Definition at line 292 of file G4NuclNuclDiffuseElastic.hh.
292{return fCofAlphaMax;};
◆ GetCosHaPit2()
Definition at line 195 of file G4NuclNuclDiffuseElastic.hh.
195{return std::cos(CLHEP::halfpi*t*t);};
Referenced by GetCint().
◆ GetCoulombElasticXsc()
|
inline |
Definition at line 514 of file G4NuclNuclDiffuseElastic.hh.
518{
519 G4double sinHalfTheta = std::sin(0.5*theta);
520 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
525 G4double k = momentum/CLHEP::hbarc;
527 G4double ch2 = ch*ch;
528 G4double xsc = ch2/((sinHalfTheta2+am)*(sinHalfTheta2+am));
529
530 return xsc;
531}
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Definition: G4NuclNuclDiffuseElastic.hh:456
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
Definition: G4NuclNuclDiffuseElastic.hh:441
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:468
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:112
Referenced by GetInvCoulombElasticXsc().
◆ GetCoulombIntegralXsc()
|
inline |
Definition at line 579 of file G4NuclNuclDiffuseElastic.hh.
582{
583 G4double c1 = std::cos(theta1);
584 //G4cout<<"c1 = "<<c1<<G4endl;
585 G4double c2 = std::cos(theta2);
586 // G4cout<<"c2 = "<<c2<<G4endl;
588 // G4cout<<"beta = "<<beta<<G4endl;
591 // G4cout<<"fZomerfeld = "<<n<<G4endl;
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;
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}
◆ GetCoulombTotalXsc()
|
inline |
Definition at line 553 of file G4NuclNuclDiffuseElastic.hh.
◆ GetDiffElasticProb()
Definition at line 397 of file G4NuclNuclDiffuseElastic.cc.
402{
403 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
404 G4double delta, diffuse, gamma;
405 G4double e1, e2, bone, bone2;
406
407 // G4double wavek = momentum/hbarc; // wave vector
408 // G4double r0 = 1.08*fermi;
409 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
410
412 G4double kr2 = kr*kr;
413 G4double krt = kr*theta;
414
415 bzero = BesselJzero(krt);
416 bzero2 = bzero*bzero;
417 bone = BesselJone(krt);
418 bone2 = bone*bone;
419 bonebyarg = BesselOneByArg(krt);
420 bonebyarg2 = bonebyarg*bonebyarg;
421
422 // VI - Coverity complains
423 /*
424 if (fParticle == theProton)
425 {
426 diffuse = 0.63*fermi;
427 gamma = 0.3*fermi;
428 delta = 0.1*fermi*fermi;
429 e1 = 0.3*fermi;
430 e2 = 0.35*fermi;
431 }
432 else // as proton, if were not defined
433 {
434 */
435 diffuse = 0.63*fermi;
436 gamma = 0.3*fermi;
437 delta = 0.1*fermi*fermi;
438 e1 = 0.3*fermi;
439 e2 = 0.35*fermi;
440 //}
442
443 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
444
446 G4double kgamma2 = kgamma*kgamma;
447
448 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
449 // G4double dk2t2 = dk2t*dk2t;
450 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
451
453
454 damp = DampFactor(pikdt);
455 damp2 = damp*damp;
456
457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
458 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
459
460
461 sigma = kgamma2;
462 // sigma += dk2t2;
463 sigma *= bzero2;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2;
466 sigma *= damp2; // *rad*rad;
467
468 return sigma;
469}
G4double BesselJzero(G4double z)
Definition: G4NuclNuclDiffuseElastic.cc:2033
G4double BesselOneByArg(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:419
G4double DampFactor(G4double z)
Definition: G4NuclNuclDiffuseElastic.hh:396
Referenced by GetDiffuseElasticXsc().
◆ GetDiffElasticSumProb()
Definition at line 477 of file G4NuclNuclDiffuseElastic.cc.
482{
483 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
484 G4double delta, diffuse, gamma;
485 G4double e1, e2, bone, bone2;
486
487 // G4double wavek = momentum/hbarc; // wave vector
488 // G4double r0 = 1.08*fermi;
489 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
490
492 G4double kr2 = kr*kr;
493 G4double krt = kr*theta;
494
495 bzero = BesselJzero(krt);
496 bzero2 = bzero*bzero;
497 bone = BesselJone(krt);
498 bone2 = bone*bone;
499 bonebyarg = BesselOneByArg(krt);
500 bonebyarg2 = bonebyarg*bonebyarg;
501
502 if (fParticle == theProton)
503 {
504 diffuse = 0.63*fermi;
505 // diffuse = 0.6*fermi;
506 gamma = 0.3*fermi;
507 delta = 0.1*fermi*fermi;
508 e1 = 0.3*fermi;
509 e2 = 0.35*fermi;
510 }
511 else // as proton, if were not defined
512 {
513 diffuse = 0.63*fermi;
514 gamma = 0.3*fermi;
515 delta = 0.1*fermi*fermi;
516 e1 = 0.3*fermi;
517 e2 = 0.35*fermi;
518 }
520 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
522
523 // G4cout<<"kgamma = "<<kgamma<<G4endl;
524
525 if(fAddCoulomb) // add Coulomb correction
526 {
527 G4double sinHalfTheta = std::sin(0.5*theta);
528 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
529
530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
532 }
533
534 G4double kgamma2 = kgamma*kgamma;
535
536 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
537 // G4cout<<"dk2t = "<<dk2t<<G4endl;
538 // G4double dk2t2 = dk2t*dk2t;
539 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
540
542
543 // G4cout<<"pikdt = "<<pikdt<<G4endl;
544
545 damp = DampFactor(pikdt);
546 damp2 = damp*damp;
547
548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
549 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
550
551 sigma = kgamma2;
552 // sigma += dk2t2;
553 sigma *= bzero2;
554 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone;
556
557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
558 sigma += kr2*bonebyarg2; // correction at J1()/()
559
560 sigma *= damp2; // *rad*rad;
561
562 return sigma;
563}
Referenced by GetDiffuseElasticSumXsc().
◆ GetDiffElasticSumProbA()
Definition at line 572 of file G4NuclNuclDiffuseElastic.cc.
573{
574 G4double theta;
575
576 theta = std::sqrt(alpha);
577
578 // theta = std::acos( 1 - alpha/2. );
579
580 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
581 G4double delta, diffuse, gamma;
582 G4double e1, e2, bone, bone2;
583
584 // G4double wavek = momentum/hbarc; // wave vector
585 // G4double r0 = 1.08*fermi;
586 // G4double rad = r0*G4Pow::GetInstance()->A13(A);
587
589 G4double kr2 = kr*kr;
590 G4double krt = kr*theta;
591
592 bzero = BesselJzero(krt);
593 bzero2 = bzero*bzero;
594 bone = BesselJone(krt);
595 bone2 = bone*bone;
596 bonebyarg = BesselOneByArg(krt);
597 bonebyarg2 = bonebyarg*bonebyarg;
598
599 if (fParticle == theProton)
600 {
601 diffuse = 0.63*fermi;
602 // diffuse = 0.6*fermi;
603 gamma = 0.3*fermi;
604 delta = 0.1*fermi*fermi;
605 e1 = 0.3*fermi;
606 e2 = 0.35*fermi;
607 }
608 else // as proton, if were not defined
609 {
610 diffuse = 0.63*fermi;
611 gamma = 0.3*fermi;
612 delta = 0.1*fermi*fermi;
613 e1 = 0.3*fermi;
614 e2 = 0.35*fermi;
615 }
617 // G4double kgamma = fWaveVector*gamma; // wavek*delta;
619
620 // G4cout<<"kgamma = "<<kgamma<<G4endl;
621
622 if(fAddCoulomb) // add Coulomb correction
623 {
625 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
626
627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
629 }
630
631 G4double kgamma2 = kgamma*kgamma;
632
633 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
634 // G4cout<<"dk2t = "<<dk2t<<G4endl;
635 // G4double dk2t2 = dk2t*dk2t;
636 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
637
639
640 // G4cout<<"pikdt = "<<pikdt<<G4endl;
641
642 damp = DampFactor(pikdt);
643 damp2 = damp*damp;
644
645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
646 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
647
648 sigma = kgamma2;
649 // sigma += dk2t2;
650 sigma *= bzero2;
651 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone;
653
654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
655 sigma += kr2*bonebyarg2; // correction at J1()/()
656
657 sigma *= damp2; // *rad*rad;
658
659 return sigma;
660}
Referenced by GetIntegrandFunction().
◆ GetDiffuseElasticSumXsc()
| G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 261 of file G4NuclNuclDiffuseElastic.cc.
265{
266 fParticle = particle;
267 fWaveVector = momentum/hbarc;
268 fAtomicWeight = A;
269 fAtomicNumber = Z;
270 fNuclearRadius = CalculateNuclearRad(A);
271 fAddCoulomb = false;
272
274
275 G4double kRt = fWaveVector*fNuclearRadius*theta;
276 G4double kRtC = 1.9;
277
278 if( z && (kRt > kRtC) )
279 {
280 fAddCoulomb = true;
281 fBeta = CalculateParticleBeta( particle, momentum);
282 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
283 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
284 }
286
287 return sigma;
288}
G4double GetDiffElasticSumProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:477
G4double CalculateNuclearRad(G4double A)
Definition: G4NuclNuclDiffuseElastic.hh:483
Referenced by GetInvElasticSumXsc().
◆ GetDiffuseElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A | ||
| ) |
Definition at line 190 of file G4NuclNuclDiffuseElastic.cc.
194{
195 fParticle = particle;
196 fWaveVector = momentum/hbarc;
197 fAtomicWeight = A;
198 fAddCoulomb = false;
199 fNuclearRadius = CalculateNuclearRad(A);
200
202
203 return sigma;
204}
G4double GetDiffElasticProb(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:397
Referenced by GetInvElasticXsc().
◆ GetErf()
Definition at line 631 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by GetErfComp(), and GetErfInt().
◆ GetErfcComp()
Definition at line 652 of file G4NuclNuclDiffuseElastic.hh.
653{
655 return erfcz;
656}
G4complex GetErfComp(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.cc:1466
◆ GetErfcInt()
Definition at line 672 of file G4NuclNuclDiffuseElastic.hh.
673{
675 return erfcz;
676}
G4complex GetErfInt(G4complex z)
Definition: G4NuclNuclDiffuseElastic.cc:1527
Referenced by AmplitudeSim(), GammaLess(), and GammaMore().
◆ GetErfComp()
Definition at line 1466 of file G4NuclNuclDiffuseElastic.cc.
1467{
1469 G4double n2, cofn, shny, chny, fn, gn;
1470
1471 G4double x = z.real();
1472 G4double y = z.imag();
1473
1474 G4double outRe = 0., outIm = 0.;
1475
1476 G4double twox = 2.*x;
1477 G4double twoxy = twox*y;
1478 G4double twox2 = twox*twox;
1479
1481
1482 G4double cos2xy = std::cos(twoxy);
1483 G4double sin2xy = std::sin(twoxy);
1484
1485 G4double twoxcos2xy = twox*cos2xy;
1486 G4double twoxsin2xy = twox*sin2xy;
1487
1489 {
1491
1493
1494 chny = std::cosh(n*y);
1495 shny = std::sinh(n*y);
1496
1497 fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1498 gn = twoxsin2xy*chny + n*cos2xy*shny;
1499
1500 fn *= cofn;
1501 gn *= cofn;
1502
1503 outRe += fn;
1504 outIm += gn;
1505 }
1506 outRe *= 2*cof1;
1507 outIm *= 2*cof1;
1508
1509 if(std::abs(x) < 0.0001)
1510 {
1511 outRe += GetErf(x);
1512 outIm += cof1*y;
1513 }
1514 else
1515 {
1516 outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
1517 outIm += cof1*sin2xy/twox;
1518 }
1520}
G4double GetErf(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:631
Referenced by GetErfcComp(), and TestErfcComp().
◆ GetErfcSer()
Definition at line 662 of file G4NuclNuclDiffuseElastic.hh.
663{
665 return erfcz;
666}
G4complex GetErfSer(G4complex z, G4int nMax)
Definition: G4NuclNuclDiffuseElastic.hh:719
◆ GetErfInt()
Definition at line 1527 of file G4NuclNuclDiffuseElastic.cc.
1528{
1529 G4double outRe, outIm;
1530
1531 G4double x = z.real();
1532 G4double y = z.imag();
1533 fReZ = x;
1534
1535 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1536
1539
1540 outRe *= 2./std::sqrt(CLHEP::pi);
1541 outIm *= 2./std::sqrt(CLHEP::pi);
1542
1543 outRe += GetErf(x);
1544
1546}
G4double GetExpSin(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:753
G4double GetExpCos(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:742
Referenced by GetErfcInt(), and TestErfcInt().
◆ GetErfSer()
Definition at line 719 of file G4NuclNuclDiffuseElastic.hh.
720{
722 G4double a =1., b = 1., tmp;
723 G4complex sum = z, d = z;
724
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}
Referenced by GetErfcSer(), and TestErfcSer().
◆ GetExpCos()
Definition at line 742 of file G4NuclNuclDiffuseElastic.hh.
743{
744 G4double result;
745
746 result = G4Exp(x*x-fReZ*fReZ);
747 result *= std::cos(2.*x*fReZ);
748 return result;
749}
Referenced by GetErfInt().
◆ GetExpSin()
Definition at line 753 of file G4NuclNuclDiffuseElastic.hh.
754{
755 G4double result;
756
757 result = G4Exp(x*x-fReZ*fReZ);
758 result *= std::sin(2.*x*fReZ);
759 return result;
760}
Referenced by GetErfInt().
◆ GetFresnelDiffuseXsc()
Definition at line 1023 of file G4NuclNuclDiffuseElastic.hh.
1024{
1027 G4double xsc = ratio*ruthXsc;
1028 return xsc;
1029}
G4double GetRatioGen(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:1972
G4double GetRutherfordXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:537
Referenced by GetFresnelIntegrandXsc().
◆ GetFresnelIntegrandXsc()
Definition at line 1035 of file G4NuclNuclDiffuseElastic.hh.
1036{
1037 G4double theta = std::sqrt(alpha);
1039 return xsc;
1040}
G4double GetFresnelDiffuseXsc(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:1023
Referenced by BuildAngleTable().
◆ GetHadronNucleonXscNS()
| G4double G4NuclNuclDiffuseElastic::GetHadronNucleonXscNS | ( | G4ParticleDefinition * | pParticle, |
| G4double | pTkin, | ||
| G4ParticleDefinition * | tParticle | ||
| ) |
Definition at line 1862 of file G4NuclNuclDiffuseElastic.cc.
1865{
1867 G4double hpXsc(0);
1868 G4double hnXsc(0);
1869
1870
1873
1874 G4double proj_energy = proj_mass + pTkin;
1875 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1876
1878
1879 sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1880 proj_momentum /= CLHEP::GeV;
1881 proj_energy /= CLHEP::GeV;
1882 proj_mass /= CLHEP::GeV;
1884
1885 // General PDG fit constants
1886
1887
1888 // fEtaRatio=Re[f(0)]/Im[f(0)]
1889
1890 if( proj_momentum >= 1.2 )
1891 {
1893 }
1894 else if( proj_momentum >= 0.6 )
1895 {
1898 }
1899 else
1900 {
1901 fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1902 }
1904
1905 // xsc
1906
1907 if( proj_momentum >= 10. ) // high energy: pp = nn = np
1908 // if( proj_momentum >= 2.)
1909 {
1910 //Delta = 1.;
1911
1912 //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1913
1914 //AR-12Aug2016 if( proj_momentum >= 10.)
1915 {
1916 B0 = 7.5;
1917 A0 = 100. - B0*G4Log(3.0e7);
1918
1919 xsection = A0 + B0*G4Log(proj_energy) - 11
1921 0.93827*0.93827,-0.165); // mb
1922 }
1923 }
1924 else // low energy pp = nn != np
1925 {
1926 if(pParticle == tParticle) // pp or nn // nn to be pp
1927 {
1928 if( proj_momentum < 0.73 )
1929 {
1931 }
1932 else if( proj_momentum < 1.05 )
1933 {
1934 hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))*
1935 (G4Log(proj_momentum/0.73));
1936 }
1937 else // if( proj_momentum < 10. )
1938 {
1939 hnXsc = 39.0 +
1941 }
1942 xsection = hnXsc;
1943 }
1944 else // pn to be np
1945 {
1946 if( proj_momentum < 0.8 )
1947 {
1949 }
1950 else if( proj_momentum < 1.4 )
1951 {
1953 }
1954 else // if( proj_momentum < 10. )
1955 {
1956 hpXsc = 33.3+
1959 }
1960 xsection = hpXsc;
1961 }
1962 }
1963 xsection *= CLHEP::millibarn; // parametrised in mb
1965 return xsection;
1966}
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
Definition: G4NuclNuclDiffuseElastic.hh:1080
Referenced by InitParametersGla().
◆ GetIntegrandFunction()
Definition at line 668 of file G4NuclNuclDiffuseElastic.cc.
669{
670 G4double result;
671
672 result = GetDiffElasticSumProbA(alpha);
673
674 // result *= 2*pi*std::sin(theta);
675
676 return result;
677}
G4double GetDiffElasticSumProbA(G4double alpha)
Definition: G4NuclNuclDiffuseElastic.cc:572
Referenced by IntegralElasticProb(), SampleThetaCMS(), and TestAngleTable().
◆ GetInvCoulombElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | tMand, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 348 of file G4NuclNuclDiffuseElastic.cc.
352{
354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
355
358 G4ParticleDefinition * theDef = 0;
359
360 if (iZ == 1 && iA == 1) theDef = theProton;
361 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
364 else if (iZ == 2 && iA == 4) theDef = theAlpha;
366
368
369 G4LorentzVector lv(0.0,0.0,0.0,tmass);
370 lv += lv1;
371
372 G4ThreeVector bst = lv.boostVector();
373 lv1.boost(-bst);
374
375 G4ThreeVector p1 = lv1.vect();
377 G4double ptot2 = ptot*ptot;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
379
380 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0;
382
383 G4double thetaCMS = std::acos(cost);
384
386
387 sigma *= pi/ptot2;
388
389 return sigma;
390}
Definition: ThreeVector.h:36
double mag() const
Definition: LorentzVector.h:67
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
Definition: G4NuclNuclDiffuseElastic.hh:514
Definition: G4ParticleDefinition.hh:60
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
◆ GetInvElasticSumXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvElasticSumXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | tMand, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 296 of file G4NuclNuclDiffuseElastic.cc.
300{
302
303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
304
307
308 G4ParticleDefinition* theDef = 0;
309
310 if (iZ == 1 && iA == 1) theDef = theProton;
311 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
314 else if (iZ == 2 && iA == 4) theDef = theAlpha;
316
318
319 G4LorentzVector lv(0.0,0.0,0.0,tmass);
320 lv += lv1;
321
322 G4ThreeVector bst = lv.boostVector();
323 lv1.boost(-bst);
324
325 G4ThreeVector p1 = lv1.vect();
327 G4double ptot2 = ptot*ptot;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
329
330 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0;
332
333 G4double thetaCMS = std::acos(cost);
334
336
337 sigma *= pi/ptot2;
338
339 return sigma;
340}
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
Definition: G4NuclNuclDiffuseElastic.cc:261
◆ GetInvElasticXsc()
| G4double G4NuclNuclDiffuseElastic::GetInvElasticXsc | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A, | ||
| G4double | Z | ||
| ) |
Definition at line 211 of file G4NuclNuclDiffuseElastic.cc.
215{
217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
218
221 G4ParticleDefinition * theDef = 0;
222
223 if (iZ == 1 && iA == 1) theDef = theProton;
224 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
227 else if (iZ == 2 && iA == 4) theDef = theAlpha;
229
231
232 G4LorentzVector lv(0.0,0.0,0.0,tmass);
233 lv += lv1;
234
235 G4ThreeVector bst = lv.boostVector();
236 lv1.boost(-bst);
237
238 G4ThreeVector p1 = lv1.vect();
240 G4double ptot2 = ptot*ptot;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
242
243 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0;
245
246 G4double thetaCMS = std::acos(cost);
247
249
250 sigma *= pi/ptot2;
251
252 return sigma;
253}
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:190
◆ GetLegendrePol()
Definition at line 1440 of file G4NuclNuclDiffuseElastic.cc.
1441{
1443 G4double x = std::cos(theta);
1444
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.;
1453 else
1454 {
1455 // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
1456
1457 legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
1458 }
1459 return legPol;
1460}
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96
Referenced by AmplitudeGla().
◆ GetNuclearRadius()
|
inline |
Definition at line 185 of file G4NuclNuclDiffuseElastic.hh.
185{return fNuclearRadius;};
◆ GetProfileLambda()
|
inline |
Definition at line 274 of file G4NuclNuclDiffuseElastic.hh.
274{return fProfileLambda;};
◆ GetRatioGen()
Definition at line 1972 of file G4NuclNuclDiffuseElastic.cc.
1973{
1974 G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1975 G4double dTheta = 0.5*(theta - fRutherfordTheta);
1976 G4double sindTheta = std::sin(dTheta);
1977
1979 G4double prof2 = prof*prof;
1980 // G4double profmod = std::abs(prof);
1981 G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1982
1983 order = std::abs(order); // since sin changes sign!
1984 // G4cout<<"order = "<<order<<G4endl;
1985
1988
1989 G4double out;
1990
1991 if ( theta <= fRutherfordTheta )
1992 {
1993 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1994 out += ( cosFresnel + sinFresnel - 1. )*prof;
1995 }
1996 else
1997 {
1998 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1999 }
2000
2001 return out;
2002}
G4double Profile(G4double theta)
Definition: G4NuclNuclDiffuseElastic.hh:908
G4double GetSint(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:782
G4double GetCint(G4double x)
Definition: G4NuclNuclDiffuseElastic.hh:767
Referenced by GetFresnelDiffuseXsc().
◆ GetRatioSim()
Definition at line 1003 of file G4NuclNuclDiffuseElastic.hh.
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;
1013
1014 G4double out = 0.5*( cosFresnel*cosFresnel + sinFresnel*sinFresnel );
1015
1016 return out;
1017}
◆ GetRutherfordXsc()
Definition at line 537 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by GetFresnelDiffuseXsc().
◆ GetScatteringAngle()
| G4double G4NuclNuclDiffuseElastic::GetScatteringAngle | ( | G4int | iMomentum, |
| G4int | iAngle, | ||
| G4double | position | ||
| ) |
Definition at line 1086 of file G4NuclNuclDiffuseElastic.cc.
1087{
1088 G4double x1, x2, y1, y2, randAngle;
1089
1090 if( iAngle == 0 )
1091 {
1092 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1093 // iAngle++;
1094 }
1095 else
1096 {
1098 {
1099 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1100 }
1101 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1102 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1103
1104 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1105 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1106
1107 if ( x1 == x2 ) randAngle = x2;
1108 else
1109 {
1111 else
1112 {
1113 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1114 }
1115 }
1116 }
1117 return randAngle;
1118}
Definition: xmltok.h:110
Referenced by SampleTableThetaCMS().
◆ GetSinHaPit2()
Definition at line 196 of file G4NuclNuclDiffuseElastic.hh.
196{return std::sin(CLHEP::halfpi*t*t);};
Referenced by GetSint().
◆ GetSint()
Definition at line 782 of file G4NuclNuclDiffuseElastic.hh.
783{
784 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
785
786 G4double out =
788
789 return out;
790}
G4double GetSinHaPit2(G4double t)
Definition: G4NuclNuclDiffuseElastic.hh:196
Referenced by GetRatioGen(), and GetRatioSim().
◆ InitDynParameters()
| void G4NuclNuclDiffuseElastic::InitDynParameters | ( | const G4ParticleDefinition * | theParticle, |
| G4double | partMom | ||
| ) |
Definition at line 1767 of file G4NuclNuclDiffuseElastic.cc.
1769{
1770 G4double a = 0.;
1773
1774 fWaveVector = partMom/CLHEP::hbarc;
1775
1777
1778 if( z )
1779 {
1780 a = partMom/m1; // beta*gamma for m1
1781 fBeta = a/std::sqrt(1+a*a);
1782 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1783 fRutherfordRatio = fZommerfeld/fWaveVector;
1784 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1785 }
1787 fProfileDelta = fCofDelta*fProfileLambda;
1788 fProfileAlpha = fCofAlpha*fProfileLambda;
1789
1790 CalculateCoulombPhaseZero();
1791 CalculateRutherfordAnglePar();
1792
1793 return;
1794}
void CalculateRutherfordAnglePar()
Definition: G4NuclNuclDiffuseElastic.hh:858
void CalculateCoulombPhaseZero()
Definition: G4NuclNuclDiffuseElastic.hh:831
Referenced by BuildAngleTable(), and SampleCoulombMuCMS().
◆ Initialise()
| void G4NuclNuclDiffuseElastic::Initialise | ( | ) |
Definition at line 150 of file G4NuclNuclDiffuseElastic.cc.
151{
152
153 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
154
157
158 // projectile radius
159
162
163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
164 {
165 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
166 fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
167
168 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
169 fNuclearRadius += R1;
170
172 {
174 <<(*theElementTable)[jEl]->GetName()<<G4endl;
175 }
176 fElementNumberVector.push_back(fAtomicNumber);
177 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
178
179 BuildAngleTable();
180 fAngleBank.push_back(fAngleTable);
181 }
182}
G4double GetAtomicMassAmu(const G4String &symb) const
Definition: G4NistManager.hh:376
void BuildAngleTable()
Definition: G4NuclNuclDiffuseElastic.cc:1007
G4int GetBaryonNumber() const
Definition: G4ParticleDefinition.hh:132
◆ InitialiseOnFly()
Definition at line 978 of file G4NuclNuclDiffuseElastic.cc.
979{
980 fAtomicNumber = Z; // atomic number
982
985
986 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
987
989 {
992 }
993 fElementNumberVector.push_back(fAtomicNumber);
994
995 BuildAngleTable();
996
997 fAngleBank.push_back(fAngleTable);
998
999 return;
1000}
Referenced by SampleTableThetaCMS().
◆ InitParameters()
| void G4NuclNuclDiffuseElastic::InitParameters | ( | const G4ParticleDefinition * | theParticle, |
| G4double | partMom, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 1722 of file G4NuclNuclDiffuseElastic.cc.
1724{
1725 fAtomicNumber = Z; // atomic number
1727
1728 fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight);
1730 fNuclearRadius1 = CalculateNuclearRad(A1);
1731 // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1732 fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1733
1734 G4double a = 0.;
1737
1738 fWaveVector = partMom/CLHEP::hbarc;
1739
1742
1743 if( z )
1744 {
1745 a = partMom/m1; // beta*gamma for m1
1746 fBeta = a/std::sqrt(1+a*a);
1747 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1748 fRutherfordRatio = fZommerfeld/fWaveVector;
1749 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1750 }
1754 fProfileDelta = fCofDelta*fProfileLambda;
1755 fProfileAlpha = fCofAlpha*fProfileLambda;
1756
1757 CalculateCoulombPhaseZero();
1758 CalculateRutherfordAnglePar();
1759
1760 return;
1761}
◆ InitParametersGla()
| void G4NuclNuclDiffuseElastic::InitParametersGla | ( | const G4DynamicParticle * | aParticle, |
| G4double | partMom, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 1802 of file G4NuclNuclDiffuseElastic.cc.
1804{
1805 fAtomicNumber = Z; // target atomic number
1807
1811 fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1812
1813
1814 G4double a = 0., kR12;
1817
1818 fWaveVector = partMom/CLHEP::hbarc;
1819
1820 G4double pN = A1 - z;
1821 if( pN < 0. ) pN = 0.;
1822
1824 if( tN < 0. ) tN = 0.;
1825
1827 pTkin /= A1;
1828
1829
1830 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1831 (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1832
1835 kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1837 fMaxL = (G4int(kR12)+1)*4;
1839
1840 if( z )
1841 {
1842 a = partMom/m1; // beta*gamma for m1
1843 fBeta = a/std::sqrt(1+a*a);
1844 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1845 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1846 }
1847
1848 CalculateCoulombPhaseZero();
1849
1850
1851 return;
1852}
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
Definition: G4NuclNuclDiffuseElastic.cc:1862
◆ IntegralElasticProb()
| G4double G4NuclNuclDiffuseElastic::IntegralElasticProb | ( | const G4ParticleDefinition * | particle, |
| G4double | theta, | ||
| G4double | momentum, | ||
| G4double | A | ||
| ) |
Definition at line 684 of file G4NuclNuclDiffuseElastic.cc.
688{
689 G4double result;
690 fParticle = particle;
691 fWaveVector = momentum/hbarc;
692 fAtomicWeight = A;
693
694 fNuclearRadius = CalculateNuclearRad(A);
695
696
697 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
698
699 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
701
702 return result;
703}
G4double GetIntegrandFunction(G4double theta)
Definition: G4NuclNuclDiffuseElastic.cc:668
◆ PhaseFar()
Definition at line 945 of file G4NuclNuclDiffuseElastic.hh.
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
955
956 return std::exp(z);
957}
Referenced by AmplitudeFar().
◆ PhaseNear()
Definition at line 927 of file G4NuclNuclDiffuseElastic.hh.
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
937
938 return std::exp(z);
939}
Referenced by AmplitudeNear().
◆ Profile()
Definition at line 908 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by GetRatioGen().
◆ ProfileFar()
Definition at line 892 of file G4NuclNuclDiffuseElastic.hh.
893{
894 G4double dTheta = fRutherfordTheta + theta;
895 G4double argument = fProfileDelta*dTheta;
896
898 result /= std::sinh(CLHEP::pi*argument);
899 result /= dTheta;
900
901 return result;
902}
Referenced by AmplitudeFar().
◆ ProfileNear()
Definition at line 871 of file G4NuclNuclDiffuseElastic.hh.
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}
Referenced by AmplitudeNear(), and AmplitudeSim().
◆ SampleCoulombMuCMS()
| G4double G4NuclNuclDiffuseElastic::SampleCoulombMuCMS | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p | ||
| ) |
Definition at line 809 of file G4NuclNuclDiffuseElastic.cc.
810{
811 G4double t(0.), rand(0.), mu(0.);
812
815
816 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
817 fNuclearRadius += R1;
818
819 InitDynParameters(fParticle, p);
820
821 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutThetaTg2);
822
823 rand = G4UniformRand();
824
825 // sample (1-cosTheta) in 0 < Theta < ThetaC as Coulomb scattering
826
827 mu = fCoulombMuC*rand*fAm;
828 mu /= fAm + fCoulombMuC*( 1. - rand );
829
830 t = 4.*p*p*mu;
831
832 return t;
833}
Referenced by SampleInvariantT().
◆ SampleInvariantT()
|
virtual |
Reimplemented from G4HadronElastic.
Definition at line 777 of file G4NuclNuclDiffuseElastic.cc.
779{
780 fParticle = aParticle;
781 fAtomicNumber = G4double(Z);
782 fAtomicWeight = G4double(A);
783
785 G4double totElab = std::sqrt(m1*m1+p*p);
787 G4LorentzVector lv1(p,0.0,0.0,totElab);
788 G4LorentzVector lv(0.0,0.0,0.0,mass2);
789 lv += lv1;
790
791 G4ThreeVector bst = lv.boostVector();
792 lv1.boost(-bst);
793
794 G4ThreeVector p1 = lv1.vect();
796
797 // t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
798
800
801
802 return t;
803}
G4double SampleCoulombMuCMS(const G4ParticleDefinition *aParticle, G4double p)
Definition: G4NuclNuclDiffuseElastic.cc:809
static G4double GetNuclearMass(const G4double A, const G4double Z)
Definition: G4NucleiProperties.cc:50
◆ SampleT()
| G4double G4NuclNuclDiffuseElastic::SampleT | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | A | ||
| ) |
Definition at line 709 of file G4NuclNuclDiffuseElastic.cc.
711{
714 return t;
715}
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:723
Referenced by SampleThetaLab().
◆ SampleTableT()
| G4double G4NuclNuclDiffuseElastic::SampleTableT | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 840 of file G4NuclNuclDiffuseElastic.cc.
842{
844 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
846 return t;
847}
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:855
◆ SampleTableThetaCMS()
| G4double G4NuclNuclDiffuseElastic::SampleTableThetaCMS | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 855 of file G4NuclNuclDiffuseElastic.cc.
857{
858 size_t iElement;
859 G4int iMomentum, iAngle;
861 G4double m1 = particle->GetPDGMass();
862
863 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
864 {
865 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
866 }
867 if ( iElement == fElementNumberVector.size() )
868 {
870
871 // iElement--;
872
873 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
874 // << " is not found, return zero angle" << G4endl;
875 // return 0.; // no table for this element
876 }
877 // G4cout<<"iElement = "<<iElement<<G4endl;
878
879 fAngleTable = fAngleBank[iElement];
880
881 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882
883 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
884 {
885 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
886
887 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
888 }
889 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
890
891
892 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
893 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
894
895
896 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
897 {
899
900 // G4cout<<"position = "<<position<<G4endl;
901
902 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
903 {
905 }
906 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
907
908 // G4cout<<"iAngle = "<<iAngle<<G4endl;
909
911
912 // G4cout<<"randAngle = "<<randAngle<<G4endl;
913 }
914 else // kinE inside between energy table edges
915 {
916 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
918
919 // G4cout<<"position = "<<position<<G4endl;
920
921 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
922 {
923 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
925 }
926 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
927
928 // G4cout<<"iAngle = "<<iAngle<<G4endl;
929
931
932 // G4cout<<"theta2 = "<<theta2<<G4endl;
933
934 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
935
936 // G4cout<<"E2 = "<<E2<<G4endl;
937
938 iMomentum--;
939
940 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
941
942 // G4cout<<"position = "<<position<<G4endl;
943
944 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
945 {
946 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
948 }
949 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
950
952
953 // G4cout<<"theta1 = "<<theta1<<G4endl;
954
955 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
956
957 // G4cout<<"E1 = "<<E1<<G4endl;
958
959 W = 1.0/(E2 - E1);
960 W1 = (E2 - kinE)*W;
961 W2 = (kinE - E1)*W;
962
963 randAngle = W1*theta1 + W2*theta2;
964
965 // randAngle = theta2;
966 }
967 // G4double angle = randAngle;
968 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
969 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
970
971 return randAngle;
972}
void InitialiseOnFly(G4double Z, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:978
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
Definition: G4NuclNuclDiffuseElastic.cc:1086
Referenced by SampleTableT().
◆ SampleThetaCMS()
| G4double G4NuclNuclDiffuseElastic::SampleThetaCMS | ( | const G4ParticleDefinition * | aParticle, |
| G4double | p, | ||
| G4double | A | ||
| ) |
Definition at line 723 of file G4NuclNuclDiffuseElastic.cc.
725{
726 G4int i, iMax = 100;
727 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
728
729 fParticle = particle;
730 fWaveVector = momentum/hbarc;
731 fAtomicWeight = A;
732
733 fNuclearRadius = CalculateNuclearRad(A);
734
735 thetaMax = 10.174/fWaveVector/fNuclearRadius;
736
738
739 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
740
741 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
742 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
743
744 norm *= G4UniformRand();
745
746 for(i = 1; i <= iMax; i++)
747 {
748 theta1 = (i-1)*thetaMax/iMax;
749 theta2 = i*thetaMax/iMax;
750 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
751
752 if ( sum >= norm )
753 {
754 result = 0.5*(theta1 + theta2);
755 break;
756 }
757 }
758 if (i > iMax ) result = 0.5*(theta1 + theta2);
759
761
762 result += G4RandGauss::shoot(0.,sigma);
763
764 if(result < 0.) result = 0.;
765 if(result > thetaMax) result = thetaMax;
766
767 return result;
768}
Referenced by SampleT().
◆ SampleThetaLab()
| G4double G4NuclNuclDiffuseElastic::SampleThetaLab | ( | const G4HadProjectile * | aParticle, |
| G4double | tmass, | ||
| G4double | A | ||
| ) |
Definition at line 1129 of file G4NuclNuclDiffuseElastic.cc.
1131{
1136 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1137 lv += lv1;
1138
1139 G4ThreeVector bst = lv.boostVector();
1140 lv1.boost(-bst);
1141
1144 G4double tmax = 4.0*ptot*ptot;
1145 G4double t = 0.0;
1146
1147
1148 //
1149 // Sample t
1150 //
1151
1152 t = SampleT( theParticle, ptot, A);
1153
1154 // NaN finder
1155 if(!(t < 0.0 || t >= 0.0))
1156 {
1158 {
1160 << " mom(GeV)= " << plab/GeV
1161 << " S-wave will be sampled"
1162 << G4endl;
1163 }
1164 t = G4UniformRand()*tmax;
1165 }
1167 {
1170 }
1171 // Sampling of angles in CM system
1172
1174 G4double cost = 1. - 2.0*t/tmax;
1175 G4double sint;
1176
1177 if( cost >= 1.0 )
1178 {
1179 cost = 1.0;
1180 sint = 0.0;
1181 }
1182 else if( cost <= -1.0)
1183 {
1184 cost = -1.0;
1185 sint = 0.0;
1186 }
1187 else
1188 {
1189 sint = std::sqrt((1.0-cost)*(1.0+cost));
1190 }
1192 {
1194 }
1195 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1196 v1 *= ptot;
1197 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1198
1199 nlv1.boost(bst);
1200
1201 G4ThreeVector np1 = nlv1.vect();
1202
1203 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1204
1206
1207 return theta;
1208}
double theta() const
HepLorentzVector & boost(double, double, double)
Definition: LorentzVector.cc:55
Hep3Vector vect() const
G4double GetTotalMomentum() const
Definition: G4HadProjectile.hh:111
const G4ParticleDefinition * GetDefinition() const
Definition: G4HadProjectile.hh:86
const G4LorentzVector & Get4Momentum() const
Definition: G4HadProjectile.hh:91
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Definition: G4NuclNuclDiffuseElastic.cc:709
◆ SetCofAlpha()
|
inline |
Definition at line 281 of file G4NuclNuclDiffuseElastic.hh.
281{fCofAlpha = pa;};
◆ SetCofAlphaCoulomb()
|
inline |
Definition at line 283 of file G4NuclNuclDiffuseElastic.hh.
283{fCofAlphaCoulomb = pa;};
◆ SetCofAlphaMax()
|
inline |
Definition at line 282 of file G4NuclNuclDiffuseElastic.hh.
282{fCofAlphaMax = pa;};
◆ SetCofDelta()
|
inline |
Definition at line 285 of file G4NuclNuclDiffuseElastic.hh.
285{fCofDelta = pa;};
◆ SetCofFar()
|
inline |
Definition at line 287 of file G4NuclNuclDiffuseElastic.hh.
287{fCofFar = pa;};
◆ SetCofLambda()
|
inline |
Definition at line 279 of file G4NuclNuclDiffuseElastic.hh.
279{fCofLambda = pa;};
◆ SetCofPhase()
|
inline |
Definition at line 286 of file G4NuclNuclDiffuseElastic.hh.
286{fCofPhase = pa;};
◆ SetEtaRatio()
|
inline |
Definition at line 288 of file G4NuclNuclDiffuseElastic.hh.
288{fEtaRatio = pa;};
◆ SetHEModelLowLimit()
|
inline |
Definition at line 377 of file G4NuclNuclDiffuseElastic.hh.
378{
379 lowEnergyLimitHE = value;
380}
◆ SetLowestEnergyLimit()
|
inline |
Definition at line 387 of file G4NuclNuclDiffuseElastic.hh.
388{
389 lowestEnergyLimit = value;
390}
◆ SetMaxL()
|
inline |
Definition at line 289 of file G4NuclNuclDiffuseElastic.hh.
289{fMaxL = l;};
◆ SetNuclearRadiusCof()
|
inline |
Definition at line 290 of file G4NuclNuclDiffuseElastic.hh.
290{fNuclearRadiusCof = r;};
◆ SetPlabLowLimit()
|
inline |
Definition at line 372 of file G4NuclNuclDiffuseElastic.hh.
373{
374 plabLowLimit = value;
375}
◆ SetProfileAlpha()
|
inline |
Definition at line 278 of file G4NuclNuclDiffuseElastic.hh.
278{fProfileAlpha = pa;};
◆ SetProfileDelta()
|
inline |
Definition at line 277 of file G4NuclNuclDiffuseElastic.hh.
277{fProfileDelta = pd;};
◆ SetProfileLambda()
|
inline |
Definition at line 276 of file G4NuclNuclDiffuseElastic.hh.
276{fProfileLambda = pl;};
◆ SetQModelLowLimit()
|
inline |
Definition at line 382 of file G4NuclNuclDiffuseElastic.hh.
383{
384 lowEnergyLimitQ = value;
385}
◆ SetRecoilKinEnergyLimit()
|
inline |
Definition at line 367 of file G4NuclNuclDiffuseElastic.hh.
368{
369 lowEnergyRecoilLimit = value;
370}
◆ TestAngleTable()
| void G4NuclNuclDiffuseElastic::TestAngleTable | ( | const G4ParticleDefinition * | theParticle, |
| G4double | partMom, | ||
| G4double | Z, | ||
| G4double | A | ||
| ) |
Definition at line 1337 of file G4NuclNuclDiffuseElastic.cc.
1339{
1340 fAtomicNumber = Z; // atomic number
1342 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1343
1344
1345
1348
1349 fElementNumberVector.push_back(fAtomicNumber);
1350
1351
1352
1353
1354 G4int i=0, j;
1356 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1357 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1358 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1360
1361 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1362
1364
1365 fWaveVector = partMom/hbarc;
1366
1367 G4double kR = fWaveVector*fNuclearRadius;
1368 G4double kR2 = kR*kR;
1371
1372 alphaMax = kRmax*kRmax/kR2;
1373
1374 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1375
1376 alphaCoulomb = kRcoul*kRcoul/kR2;
1377
1378 if( z )
1379 {
1380 a = partMom/m1; // beta*gamma for m1
1381 fBeta = a/std::sqrt(1+a*a);
1382 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1383 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1384 }
1386
1387 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1388
1389
1390 fAddCoulomb = false;
1391
1392 for(j = 1; j < fAngleBin; j++)
1393 {
1394 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1395 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1396
1397 alpha1 = alphaMax*(j-1)/fAngleBin;
1398 alpha2 = alphaMax*( j )/fAngleBin;
1399
1400 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1401
1402 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1403 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1405 alpha1, alpha2,epsilon);
1406
1407 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1408 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1409
1410 sumL10 += deltaL10;
1411 sumL96 += deltaL96;
1412 sumAG += deltaAG;
1413
1416
1418 }
1419 fAngleTable->insertAt(i,angleVector);
1420 fAngleBank.push_back(fAngleTable);
1421
1422 /*
1423 // Integral over all angle range - Bad accuracy !!!
1424
1425 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1426 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1427 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1428 0., alpha2,epsilon);
1429 G4cout<<G4endl;
1430 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1431 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1432 */
1433 return;
1434}
◆ TestErfcComp()
Definition at line 683 of file G4NuclNuclDiffuseElastic.hh.
◆ TestErfcInt()
◆ TestErfcSer()
◆ ThetaCMStoThetaLab()
| G4double G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab | ( | const G4DynamicParticle * | aParticle, |
| G4double | tmass, | ||
| G4double | thetaCMS | ||
| ) |
Definition at line 1217 of file G4NuclNuclDiffuseElastic.cc.
1219{
1222 // G4double plab = aParticle->GetTotalMomentum();
1224 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1225
1226 lv += lv1;
1227
1228 G4ThreeVector bst = lv.boostVector();
1229
1230 lv1.boost(-bst);
1231
1234
1236 G4double cost = std::cos(thetaCMS);
1237 G4double sint;
1238
1239 if( cost >= 1.0 )
1240 {
1241 cost = 1.0;
1242 sint = 0.0;
1243 }
1244 else if( cost <= -1.0)
1245 {
1246 cost = -1.0;
1247 sint = 0.0;
1248 }
1249 else
1250 {
1251 sint = std::sqrt((1.0-cost)*(1.0+cost));
1252 }
1254 {
1256 }
1257 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1258 v1 *= ptot;
1259 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1260
1261 nlv1.boost(bst);
1262
1263 G4ThreeVector np1 = nlv1.vect();
1264
1265
1267
1268 return thetaLab;
1269}
G4LorentzVector Get4Momentum() const
◆ ThetaLabToThetaCMS()
| G4double G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS | ( | const G4DynamicParticle * | aParticle, |
| G4double | tmass, | ||
| G4double | thetaLab | ||
| ) |
Definition at line 1278 of file G4NuclNuclDiffuseElastic.cc.
1280{
1285 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1286
1287 lv += lv1;
1288
1289 G4ThreeVector bst = lv.boostVector();
1290
1291 // lv1.boost(-bst);
1292
1293 // G4ThreeVector p1 = lv1.vect();
1294 // G4double ptot = p1.mag();
1295
1297 G4double cost = std::cos(thetaLab);
1298 G4double sint;
1299
1300 if( cost >= 1.0 )
1301 {
1302 cost = 1.0;
1303 sint = 0.0;
1304 }
1305 else if( cost <= -1.0)
1306 {
1307 cost = -1.0;
1308 sint = 0.0;
1309 }
1310 else
1311 {
1312 sint = std::sqrt((1.0-cost)*(1.0+cost));
1313 }
1315 {
1317 }
1318 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1319 v1 *= plab;
1320 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1321
1322 nlv1.boost(-bst);
1323
1324 G4ThreeVector np1 = nlv1.vect();
1325
1326
1328
1329 return thetaCMS;
1330}
G4double GetTotalMomentum() const
The documentation for this class was generated from the following files:
