73 lowEnergyRecoilLimit = 100.*
keV;
74 lowEnergyLimitQ = 0.0*
GeV;
75 lowEnergyLimitHE = 0.0*
GeV;
76 lowestEnergyLimit= 0.0*
keV;
77 plabLowLimit = 20.0*
MeV;
109 if(fEnergyVector)
delete fEnergyVector;
131 for(jEl = 0 ; jEl < numOfEl; ++jEl)
133 fAtomicNumber = (*theElementTable)[jEl]->GetZ();
134 fAtomicWeight = (*theElementTable)[jEl]->GetN();
139 G4cout<<
"G4DiffuseElastic::Initialise() the element: "
140 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
142 fElementNumberVector.push_back(fAtomicNumber);
143 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
146 fAngleBank.push_back(fAngleTable);
161 fParticle = particle;
162 fWaveVector = momentum/
hbarc;
189 if (iZ == 1 && iA == 1) theDef = theProton;
190 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
192 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
193 else if (iZ == 2 && iA == 4) theDef = theAlpha;
207 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
209 if( cost >= 1.0 ) cost = 1.0;
210 else if( cost <= -1.0) cost = -1.0;
212 G4double thetaCMS = std::acos(cost);
232 fParticle = particle;
233 fWaveVector = momentum/
hbarc;
241 G4double kRt = fWaveVector*fNuclearRadius*theta;
244 if( z && (kRt > kRtC) )
249 fAm =
CalculateAm( momentum, fZommerfeld, fAtomicNumber);
276 if (iZ == 1 && iA == 1) theDef = theProton;
277 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
279 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
280 else if (iZ == 2 && iA == 4) theDef = theAlpha;
294 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
296 if( cost >= 1.0 ) cost = 1.0;
297 else if( cost <= -1.0) cost = -1.0;
299 G4double thetaCMS = std::acos(cost);
326 if (iZ == 1 && iA == 1) theDef = theProton;
327 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
329 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
330 else if (iZ == 2 && iA == 4) theDef = theAlpha;
344 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
346 if( cost >= 1.0 ) cost = 1.0;
347 else if( cost <= -1.0) cost = -1.0;
349 G4double thetaCMS = std::acos(cost);
369 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
377 G4double kr = fWaveVector*fNuclearRadius;
382 bzero2 = bzero*bzero;
386 bonebyarg2 = bonebyarg*bonebyarg;
388 if (fParticle == theProton)
390 diffuse = 0.63*
fermi;
398 diffuse = 0.63*
fermi;
408 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
415 G4double pikdt = lambda*(1.-std::exp(-
pi*fWaveVector*diffuse*theta/lambda));
420 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
421 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
427 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
428 sigma += kr2*bonebyarg2;
446 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
454 G4double kr = fWaveVector*fNuclearRadius;
459 bzero2 = bzero*bzero;
463 bonebyarg2 = bonebyarg*bonebyarg;
465 if (fParticle == theProton)
467 diffuse = 0.63*
fermi;
476 diffuse = 0.63*
fermi;
484 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
490 G4double sinHalfTheta = std::sin(0.5*theta);
491 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
493 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm);
504 G4double pikdt = lambda*(1.-std::exp(-
pi*fWaveVector*diffuse*theta/lambda));
511 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
512 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
517 sigma += mode2k2*bone2;
518 sigma += e2dk3t*bzero*bone;
521 sigma += kr2*bonebyarg2;
539 theta = std::sqrt(alpha);
543 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
551 G4double kr = fWaveVector*fNuclearRadius;
556 bzero2 = bzero*bzero;
560 bonebyarg2 = bonebyarg*bonebyarg;
562 if (fParticle == theProton)
564 diffuse = 0.63*
fermi;
573 diffuse = 0.63*
fermi;
581 G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));
588 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
590 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm);
601 G4double pikdt = lambda*(1.-std::exp(-
pi*fWaveVector*diffuse*theta/lambda));
608 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
609 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
614 sigma += mode2k2*bone2;
615 sigma += e2dk3t*bzero*bone;
618 sigma += kr2*bonebyarg2;
653 fParticle = particle;
654 fWaveVector = momentum/
hbarc;
675 G4double t = 2*p*p*( 1 - std::cos(theta) );
691 fParticle = particle;
692 fWaveVector = momentum/
hbarc;
697 thetaMax = 10.174/fWaveVector/fNuclearRadius;
699 if (thetaMax >
pi) thetaMax =
pi;
708 for(i = 1; i <= iMax; i++)
710 theta1 = (i-1)*thetaMax/iMax;
711 theta2 = i*thetaMax/iMax;
716 result = 0.5*(theta1 + theta2);
720 if (i > iMax ) result = 0.5*(theta1 + theta2);
726 if(result < 0.) result = 0.;
727 if(result > thetaMax) result = thetaMax;
741 fParticle = aParticle;
743 G4double totElab = std::sqrt(m1*m1+p*p);
782 G4int iMomentum, iAngle;
786 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
788 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5)
break;
790 if ( iElement == fElementNumberVector.size() )
802 fAngleTable = fAngleBank[iElement];
804 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
806 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
808 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
810 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
811 if ( iMomentum < 0 ) iMomentum = 0;
815 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
817 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*
G4UniformRand();
821 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
823 if( position < (*(*fAngleTable)(iMomentum))(iAngle) )
break;
825 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
840 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
843 if( position > (*(*fAngleTable)(iMomentum))(iAngle) )
break;
845 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
862 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
865 if( position > (*(*fAngleTable)(iMomentum))(iAngle) )
break;
867 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
881 randAngle = W1*theta1 + W2*theta2;
905 G4cout<<
"G4DiffuseElastic::Initialise() the element with Z = "
906 <<Z<<
"; and A = "<<A<<
G4endl;
908 fElementNumberVector.push_back(fAtomicNumber);
912 fAngleBank.push_back(fAngleTable);
926 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
932 for( i = 0; i < fEnergyBin; i++)
935 partMom = std::sqrt( kinE*(kinE + 2*m1) );
937 fWaveVector = partMom/
hbarc;
939 G4double kR = fWaveVector*fNuclearRadius;
946 alphaMax = kRmax*kRmax/kR2;
948 if (alphaMax > 4.) alphaMax = 4.;
950 alphaCoulomb = kRcoul*kRcoul/kR2;
955 fBeta = a/std::sqrt(1+a*a);
957 fAm =
CalculateAm( partMom, fZommerfeld, fAtomicNumber);
963 G4double delth = alphaMax/fAngleBin;
971 for(j = fAngleBin-1; j >= 1; j--)
979 alpha1 = delth*(j-1);
981 alpha2 = alpha1 + delth;
984 if( ( alpha1 < alphaCoulomb ) &&
z ) fAddCoulomb =
false;
991 angleVector->
PutValue( j-1 , alpha1, sum );
994 fAngleTable->
insertAt(i,angleVector);
1013 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1018 if ( iAngle >=
G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1020 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1022 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1023 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1025 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1026 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1028 if ( x1 == x2 ) randAngle =
x2;
1034 randAngle = x1 + ( position -
y1 )*( x2 - x1 )/( y2 -
y1 );
1073 t =
SampleT( theParticle, ptot, A);
1076 if(!(t < 0.0 || t >= 0.0))
1080 G4cout <<
"G4DiffuseElastic:WARNING: A = " << A
1081 <<
" mom(GeV)= " << plab/
GeV
1082 <<
" S-wave will be sampled"
1089 G4cout <<
" t= " << t <<
" tmax= " << tmax
1090 <<
" ptot= " << ptot <<
G4endl;
1103 else if( cost <= -1.0)
1110 sint = std::sqrt((1.0-cost)*(1.0+cost));
1114 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1116 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1157 G4double cost = std::cos(thetaCMS);
1165 else if( cost <= -1.0)
1172 sint = std::sqrt((1.0-cost)*(1.0+cost));
1176 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1178 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1217 G4double cost = std::cos(thetaLab);
1225 else if( cost <= -1.0)
1232 sint = std::sqrt((1.0-cost)*(1.0+cost));
1236 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1238 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1266 G4cout<<
"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1267 <<Z<<
"; and A = "<<A<<
G4endl;
1269 fElementNumberVector.push_back(fAtomicNumber);
1276 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1277 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1278 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1285 fWaveVector = partMom/
hbarc;
1287 G4double kR = fWaveVector*fNuclearRadius;
1292 alphaMax = kRmax*kRmax/kR2;
1294 if (alphaMax > 4.) alphaMax = 4.;
1296 alphaCoulomb = kRcoul*kRcoul/kR2;
1301 fBeta = a/std::sqrt(1+a*a);
1303 fAm =
CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1310 fAddCoulomb =
false;
1312 for(j = 1; j < fAngleBin; j++)
1317 alpha1 = alphaMax*(j-1)/fAngleBin;
1318 alpha2 = alphaMax*( j )/fAngleBin;
1320 if( ( alpha2 > alphaCoulomb ) &&
z ) fAddCoulomb =
true;
1325 alpha1, alpha2,epsilon);
1335 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1337 angleVector->
PutValue( j-1 , alpha1, sumL10 );
1339 fAngleTable->
insertAt(i,angleVector);
1340 fAngleBank.push_back(fAngleTable);
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
ThreeVector shoot(const G4int Ap, const G4int Af)
Hep3Vector boostVector() const
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double CalculateNuclearRad(G4double A)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
G4double BesselJzero(G4double z)
void InitialiseOnFly(G4double Z, G4double A)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4double G4NeutronHPJENDLHEData::G4double result
G4ParticleDefinition * GetDefinition() const
G4double GetDiffElasticSumProb(G4double theta)
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double GetLowEdgeEnergy(size_t binNumber) const
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double BesselJone(G4double z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
G4double GetTotalMomentum() const
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
const G4ParticleDefinition * GetDefinition() const
virtual ~G4DiffuseElastic()
G4double GetDiffElasticProb(G4double theta)
HepLorentzVector & boost(double, double, double)
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
static G4Neutron * Neutron()
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
G4LorentzVector Get4Momentum() const
G4double DampFactor(G4double z)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetDiffElasticSumProbA(G4double alpha)
void insertAt(size_t, G4PhysicsVector *)
void SetMaxEnergy(const G4double anEnergy)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double BesselOneByArg(G4double z)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
std::vector< G4Element * > G4ElementTable
static G4ElementTable * GetElementTable()
G4double GetPDGCharge() const
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetIntegrandFunction(G4double theta)
G4double GetTotalMomentum() const