123 for ( std::vector<G4PhysicsTable*>::iterator it =
fAngleBank.begin();
126 if ( (*it) ) (*it)->clearAndDestroy();
147 for( jEl = 0; jEl < numOfEl; ++jEl)
155 G4cout<<
"G4DiffuseElastic::Initialise() the element: "
156 <<(*theElementTable)[jEl]->GetName()<<
G4endl;
205 if (iZ == 1 && iA == 1) theDef =
theProton;
208 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
209 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
223 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
225 if( cost >= 1.0 ) cost = 1.0;
226 else if( cost <= -1.0) cost = -1.0;
228 G4double thetaCMS = std::acos(cost);
260 if( z && (kRt > kRtC) )
292 if (iZ == 1 && iA == 1) theDef =
theProton;
295 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
296 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
310 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
312 if( cost >= 1.0 ) cost = 1.0;
313 else if( cost <= -1.0) cost = -1.0;
315 G4double thetaCMS = std::acos(cost);
342 if (iZ == 1 && iA == 1) theDef =
theProton;
345 else if (iZ == 2 && iA == 3) theDef =
G4He3::He3();
346 else if (iZ == 2 && iA == 4) theDef =
theAlpha;
360 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
362 if( cost >= 1.0 ) cost = 1.0;
363 else if( cost <= -1.0) cost = -1.0;
365 G4double thetaCMS = std::acos(cost);
385 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
395 diffuse = 0.63*
fermi;
403 diffuse = 0.63*
fermi;
414 diffuse = 0.63*
fermi;
425 bzero2 = bzero*bzero;
429 bonebyarg2 = bonebyarg*bonebyarg;
454 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
455 sigma += kr2*bonebyarg2;
473 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
486 bzero2 = bzero*bzero;
490 bonebyarg2 = bonebyarg*bonebyarg;
494 diffuse = 0.63*
fermi;
503 diffuse = 0.63*
fermi;
514 diffuse = 0.63*
fermi;
528 G4double sinHalfTheta = std::sin(0.5*theta);
529 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
555 sigma += mode2k2*bone2;
556 sigma += e2dk3t*bzero*bone;
559 sigma += kr2*bonebyarg2;
577 theta = std::sqrt(alpha);
581 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
594 bzero2 = bzero*bzero;
598 bonebyarg2 = bonebyarg*bonebyarg;
602 diffuse = 0.63*
fermi;
611 diffuse = 0.63*
fermi;
622 diffuse = 0.63*
fermi;
637 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
662 sigma += mode2k2*bone2;
663 sigma += e2dk3t*bzero*bone;
666 sigma += kr2*bonebyarg2;
723 G4double t = 2*p*p*( 1 - std::cos(theta) );
737 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
747 if (thetaMax >
pi) thetaMax =
pi;
756 for(i = 1; i <= iMax; i++)
758 theta1 = (i-1)*thetaMax/iMax;
759 theta2 = i*thetaMax/iMax;
764 result = 0.5*(theta1 + theta2);
768 if (i > iMax ) result = 0.5*(theta1 + theta2);
774 if(result < 0.) result = 0.;
775 if(result > thetaMax) result = thetaMax;
791 G4double totElab = std::sqrt(m1*m1+p*p);
806 G4double pCMS2 = momentumCMS*momentumCMS;
807 G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
842 G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );
857 G4int iMomentum, iAngle;
879 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
881 for( iMomentum = 0; iMomentum <
fEnergyBin; iMomentum++)
883 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) )
break;
885 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
886 if ( iMomentum < 0 ) iMomentum = 0;
890 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
896 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
898 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
900 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
915 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
918 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
920 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
937 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
940 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
942 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
956 randAngle = W1*theta1 + W2*theta2;
964 if(randAngle < 0.) randAngle = 0.;
982 G4cout<<
"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
983 <<Z<<
"; and A = "<<A<<
G4endl;
1003 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1012 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1023 alphaMax = kRmax*kRmax/kR2;
1035 alphaCoulomb = kRcoul*kRcoul/kR2;
1040 fBeta = a/std::sqrt(1+a*a);
1064 alpha1 = delth*(j-1);
1066 alpha2 = alpha1 + delth;
1069 if( ( alpha1 < alphaCoulomb ) &&
z )
fAddCoulomb =
false;
1076 angleVector->
PutValue( j-1 , alpha1, sum );
1094 G4double x1, x2, y1, y2, randAngle;
1098 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1105 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1107 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1108 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1110 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1111 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1113 if ( x1 == x2 ) randAngle = x2;
1116 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*
G4UniformRand();
1119 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1158 t =
SampleT( theParticle, ptot, A);
1161 if(!(t < 0.0 || t >= 0.0))
1165 G4cout <<
"G4DiffuseElastic:WARNING: A = " << A
1166 <<
" mom(GeV)= " << plab/
GeV
1167 <<
" S-wave will be sampled"
1174 G4cout <<
" t= " << t <<
" tmax= " << tmax
1175 <<
" ptot= " << ptot <<
G4endl;
1188 else if( cost <= -1.0)
1195 sint = std::sqrt((1.0-cost)*(1.0+cost));
1199 G4cout <<
"cos(t)=" << cost <<
" std::sin(t)=" << sint <<
G4endl;
1201 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1203 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1242 G4double cost = std::cos(thetaCMS);
1250 else if( cost <= -1.0)
1257 sint = std::sqrt((1.0-cost)*(1.0+cost));
1261 G4cout <<
"cos(tcms)=" << cost <<
" std::sin(tcms)=" << sint <<
G4endl;
1263 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1265 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1302 G4double cost = std::cos(thetaLab);
1310 else if( cost <= -1.0)
1317 sint = std::sqrt((1.0-cost)*(1.0+cost));
1321 G4cout <<
"cos(tlab)=" << cost <<
" std::sin(tlab)=" << sint <<
G4endl;
1323 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1325 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1351 G4cout<<
"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1352 <<Z<<
"; and A = "<<A<<
G4endl;
1361 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1362 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1363 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1377 alphaMax = kRmax*kRmax/kR2;
1379 if (alphaMax > 4.) alphaMax = 4.;
1381 alphaCoulomb = kRcoul*kRcoul/kR2;
1386 fBeta = a/std::sqrt(1+a*a);
1402 alpha1 = alphaMax*(j-1)/fAngleBin;
1403 alpha2 = alphaMax*( j )/fAngleBin;
1405 if( ( alpha2 > alphaCoulomb ) &&
z )
fAddCoulomb =
true;
1410 alpha1, alpha2,epsilon);
1420 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1422 angleVector->
PutValue( j-1 , alpha1, sumL10 );
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)
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 lowEnergyLimitHE
G4double BesselJzero(G4double z)
CLHEP::Hep3Vector G4ThreeVector
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)
std::vector< G4String > fElementNameVector
G4double NeutronTuniform(G4int Z)
G4ParticleDefinition * theProton
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)
static G4NistManager * Instance()
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
const G4ParticleDefinition * thePionPlus
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
const G4ParticleDefinition * thePionMinus
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
static size_t GetNumberOfElements()
const G4ParticleDefinition * GetDefinition() const
virtual ~G4DiffuseElastic()
G4double GetDiffElasticProb(G4double theta)
std::vector< G4PhysicsTable * > fAngleBank
static const double twopi
static G4Triton * Triton()
static G4Proton * Proton()
static G4PionPlus * PionPlus()
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double lowestEnergyLimit
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
static G4Neutron * Neutron()
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
G4LorentzVector Get4Momentum() const
G4ParticleDefinition * theAlpha
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double DampFactor(G4double z)
G4PhysicsLogVector * fEnergyVector
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double lowEnergyRecoilLimit
static G4PionMinus * PionMinus()
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetDiffElasticSumProbA(G4double alpha)
G4double GetAtomicMassAmu(const G4String &symb) const
static const double degree
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)
G4ParticleDefinition * theDeuteron
G4double BesselOneByArg(G4double z)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
const G4ParticleDefinition * fParticle
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
static const G4double alpha
std::vector< G4double > fElementNumberVector
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4PhysicsTable * fAngleTable
double epsilon(double density, double temperature)
static const double fermi
G4double GetIntegrandFunction(G4double theta)
G4double GetTotalMomentum() const
CLHEP::HepLorentzVector G4LorentzVector
G4ParticleDefinition * theNeutron