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) );
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;
843 G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );
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;
886 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;
887 if ( iMomentum < 0 ) iMomentum = 0;
891 if (iMomentum == fEnergyBin -1 || iMomentum == 0 )
897 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
899 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
901 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
916 for(iAngle = 0; iAngle <
fAngleBin-1; iAngle++)
919 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
921 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
938 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
941 if( position > (*(*
fAngleTable)(iMomentum))(iAngle) )
break;
943 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
957 randAngle = W1*theta1 + W2*theta2;
965 if(randAngle < 0.) randAngle = 0.;
983 G4cout<<
"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
984 <<Z<<
"; and A = "<<A<<
G4endl;
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);
1067 alpha2 = alpha1 + delth;
1070 if( ( alpha1 < alphaCoulomb ) &&
z )
fAddCoulomb =
false;
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;
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 = "
1353 <<Z<<
"; and A = "<<A<<
G4endl;
1362 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
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);
1403 alpha1 = alphaMax*(j-1)/fAngleBin;
1404 alpha2 = alphaMax*( j )/fAngleBin;
1406 if( ( alpha2 > alphaCoulomb ) &&
z )
fAddCoulomb =
true;
1411 alpha1, alpha2,epsilon);
1421 <<sumL10<<
"\t"<<sumL96<<
"\t"<<sumAG<<
G4endl;
1423 angleVector->
PutValue( j-1 , alpha1, sumL10 );
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double BesselOneByArg(G4double z)
static G4PionMinus * PionMinus()
G4double GetDiffElasticProb(G4double theta)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
void SetMinEnergy(G4double anEnergy)
std::vector< ExP01TrackerHit * > a
void SetMaxEnergy(const G4double anEnergy)
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
G4double GetLowEdgeEnergy(size_t binNumber) const
Float_t y1[n_points_granero]
static constexpr double keV
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
Float_t x1[n_points_granero]
G4PhysicsLogVector * fEnergyVector
static constexpr double hbarc
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
static G4PionPlus * PionPlus()
G4IonTable * GetIonTable() const
G4double GetPDGCharge() const
G4double GetTotalMomentum() const
G4PhysicsTable * fAngleTable
G4double lowEnergyRecoilLimit
const G4ParticleDefinition * fParticle
G4double NeutronTuniform(G4int Z)
std::vector< G4String > fElementNameVector
static G4Proton * Proton()
void insertAt(size_t, G4PhysicsVector *)
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
static G4Deuteron * Deuteron()
virtual ~G4DiffuseElastic()
Float_t y2[n_points_geant4]
G4double lowestEnergyLimit
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
static constexpr double TeV
static constexpr double fermi
G4ParticleDefinition * GetDefinition() const
G4double GetDiffElasticSumProb(G4double theta)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4double GetIntegrandFunction(G4double theta)
G4double GetDiffElasticSumProbA(G4double alpha)
std::vector< G4Element * > G4ElementTable
static const G4double alpha
static constexpr double twopi
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4LorentzVector Get4Momentum() const
double A(double temperature)
G4ParticleDefinition * theDeuteron
static G4double GetNuclearMass(const G4double A, const G4double Z)
std::vector< G4PhysicsTable * > fAngleBank
ThreeVector shoot(const G4int Ap, const G4int Af)
const G4LorentzVector & Get4Momentum() const
const G4ParticleDefinition * thePionPlus
G4double G4ParticleHPJENDLHEData::G4double result
double epsilon(double density, double temperature)
const G4ParticleDefinition * GetDefinition() const
static G4Triton * Triton()
static G4Neutron * Neutron()
G4ParticleDefinition * theNeutron
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double BesselJone(G4double z)
G4double BesselJzero(G4double z)
void InitialiseOnFly(G4double Z, G4double A)
static constexpr double degree
std::vector< G4double > fElementNumberVector
G4double DampFactor(G4double z)
G4GLOB_DLL std::ostream G4cout
G4double lowEnergyLimitHE
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
Hep3Vector boostVector() const
static size_t GetNumberOfElements()
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetTotalMomentum() const
static constexpr double pi
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
G4ParticleDefinition * theProton
G4double CalculateNuclearRad(G4double A)
G4ParticleDefinition * theAlpha
Float_t x2[n_points_geant4]
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
static constexpr double GeV
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
static constexpr double pi
HepLorentzVector & boost(double, double, double)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
const G4ParticleDefinition * thePionMinus
static G4NistManager * Instance()
static G4ElementTable * GetElementTable()