124 for ( std::vector<G4PhysicsTable*>::iterator it =
fAngleBank.begin();
127 if ( (*it) ) (*it)->clearAndDestroy();
148 for( jEl = 0; jEl < numOfEl; ++jEl)
156 G4cout<<
"G4DiffuseElastic::Initialise() the element: "
157 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
206 if (iZ == 1 && iA == 1) theDef =
theProton;
209 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
210 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
224 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
226 if( cost >= 1.0 ) cost = 1.0;
227 else if( cost <= -1.0) cost = -1.0;
229 G4double thetaCMS = std::acos(cost);
261 if( z && (kRt > kRtC) )
293 if (iZ == 1 && iA == 1) theDef =
theProton;
296 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
297 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
311 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
313 if( cost >= 1.0 ) cost = 1.0;
314 else if( cost <= -1.0) cost = -1.0;
316 G4double thetaCMS = std::acos(cost);
343 if (iZ == 1 && iA == 1) theDef =
theProton;
346 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
347 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
361 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
363 if( cost >= 1.0 ) cost = 1.0;
364 else if( cost <= -1.0) cost = -1.0;
366 G4double thetaCMS = std::acos(cost);
386 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
396 diffuse = 0.63*
fermi;
404 diffuse = 0.63*
fermi;
415 diffuse = 0.63*
fermi;
426 bzero2 = bzero*bzero;
430 bonebyarg2 = bonebyarg*bonebyarg;
455 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
456 sigma += kr2*bonebyarg2;
474 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
487 bzero2 = bzero*bzero;
491 bonebyarg2 = bonebyarg*bonebyarg;
495 diffuse = 0.63*
fermi;
504 diffuse = 0.63*
fermi;
515 diffuse = 0.63*
fermi;
529 G4double sinHalfTheta = std::sin(0.5*theta);
530 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
556 sigma += mode2k2*bone2;
557 sigma += e2dk3t*bzero*bone;
560 sigma += kr2*bonebyarg2;
578 theta = std::sqrt(
alpha);
582 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
595 bzero2 = bzero*bzero;
599 bonebyarg2 = bonebyarg*bonebyarg;
603 diffuse = 0.63*
fermi;
612 diffuse = 0.63*
fermi;
623 diffuse = 0.63*
fermi;
638 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
663 sigma += mode2k2*bone2;
664 sigma += e2dk3t*bzero*bone;
667 sigma += kr2*bonebyarg2;
724 G4double t = 2*p*p*( 1 - std::cos(theta) );
738 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
748 if (thetaMax >
pi) thetaMax =
pi;
757 for(i = 1; i <= iMax; i++)
759 theta1 = (i-1)*thetaMax/iMax;
760 theta2 = i*thetaMax/iMax;
765 result = 0.5*(theta1 + theta2);
769 if (i > iMax ) result = 0.5*(theta1 + theta2);
775 if(result < 0.) result = 0.;
776 if(result > thetaMax) result = thetaMax;
792 G4double totElab = std::sqrt(m1*m1+p*p);
807 G4double pCMS2 = momentumCMS*momentumCMS;
808 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
858 G4int iMomentum, iAngle;
880 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882 for( iMomentum = 0; iMomentum <
fEnergyBin; iMomentum++)
884 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
887 if ( iMomentum < 0 ) iMomentum = 0;
891 if (iMomentum ==
fEnergyBin -1 || iMomentum == 0 )
897 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
916 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
938 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
957 randAngle = W1*theta1 + W2*theta2;
965 if(randAngle < 0.) randAngle = 0.;
983 G4cout<<
"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
1004 G4double alpha1,
alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1013 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1024 alphaMax = kRmax*kRmax/kR2;
1036 alphaCoulomb = kRcoul*kRcoul/kR2;
1041 fBeta = a/std::sqrt(1+a*a);
1065 alpha1 = delth*(j-1);
1070 if( ( alpha1 < alphaCoulomb ) && z )
fAddCoulomb =
false;
1077 angleVector->
PutValue( j-1 , alpha1, sum );
1095 G4double x1, x2, y1, y2, randAngle;
1099 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1106 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1108 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1109 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1111 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1112 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1114 if ( x1 == x2 ) randAngle = x2;
1117 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*
G4UniformRand();
1120 randAngle = x1 + (
position - y1 )*( x2 - x1 )/( y2 - y1 );
1159 t =
SampleT( theParticle, ptot,
A);
1162 if(!(t < 0.0 || t >= 0.0))
1166 G4cout <<
"G4DiffuseElastic:WARNING: A = " <<
A
1167 <<
" mom(GeV)= " << plab/
GeV
1168 <<
" S-wave will be sampled"
1175 G4cout <<
" t= " << t <<
" tmax= " << tmax
1176 <<
" ptot= " << ptot <<
G4endl;
1189 else if( cost <= -1.0)
1196 sint = std::sqrt((1.0-cost)*(1.0+cost));
1200 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1202 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1243 G4double cost = std::cos(thetaCMS);
1251 else if( cost <= -1.0)
1258 sint = std::sqrt((1.0-cost)*(1.0+cost));
1262 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1264 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1303 G4double cost = std::cos(thetaLab);
1311 else if( cost <= -1.0)
1318 sint = std::sqrt((1.0-cost)*(1.0+cost));
1322 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1324 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1352 G4cout<<
"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1363 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1364 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1378 alphaMax = kRmax*kRmax/kR2;
1380 if (alphaMax > 4.) alphaMax = 4.;
1382 alphaCoulomb = kRcoul*kRcoul/kR2;
1387 fBeta = a/std::sqrt(1+a*a);
1421 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1423 angleVector->
PutValue( j-1 , alpha1, sumL10 );
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
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
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
HepLorentzVector & boost(double, double, double)
static G4Deuteron * Deuteron()
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4ParticleDefinition * theProton
G4ParticleDefinition * theNeutron
G4double GetDiffElasticSumProbA(G4double alpha)
G4double BesselOneByArg(G4double z)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
std::vector< G4double > fElementNumberVector
G4double lowestEnergyLimit
G4ParticleDefinition * theDeuteron
const G4ParticleDefinition * thePionMinus
G4ParticleDefinition * theAlpha
G4double lowEnergyRecoilLimit
G4PhysicsLogVector * fEnergyVector
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
const G4ParticleDefinition * fParticle
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
const G4ParticleDefinition * thePionPlus
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double BesselJzero(G4double z)
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4PhysicsTable * fAngleTable
G4double GetIntegrandFunction(G4double theta)
G4double DampFactor(G4double z)
std::vector< G4String > fElementNameVector
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double lowEnergyLimitHE
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double GetDiffElasticSumProb(G4double theta)
G4double GetDiffElasticProb(G4double theta)
std::vector< G4PhysicsTable * > fAngleBank
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double CalculateNuclearRad(G4double A)
virtual ~G4DiffuseElastic()
void InitialiseOnFly(G4double Z, G4double A)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double BesselJone(G4double z)
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double NeutronTuniform(G4int Z)
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() 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
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetPDGMass() const
G4double GetPDGCharge() 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 G4Proton * Proton()
static G4Triton * Triton()
static constexpr double pi
ThreeVector shoot(const G4int Ap, const G4int Af)