Geant4  10.02.p03
G4DiffuseElastic Class Reference

#include <G4DiffuseElastic.hh>

Inheritance diagram for G4DiffuseElastic:
Collaboration diagram for G4DiffuseElastic:

Public Member Functions

 G4DiffuseElastic ()
 
virtual ~G4DiffuseElastic ()
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
void Initialise ()
 
void InitialiseOnFly (G4double Z, G4double A)
 
void BuildAngleTable ()
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
G4double NeutronTuniform (G4int Z)
 
void SetPlabLowLimit (G4double value)
 
void SetHEModelLowLimit (G4double value)
 
void SetQModelLowLimit (G4double value)
 
void SetLowestEnergyLimit (G4double value)
 
void SetRecoilKinEnergyLimit (G4double value)
 
G4double SampleT (const G4ParticleDefinition *aParticle, G4double p, G4double A)
 
G4double SampleTableT (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
 
G4double SampleThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double A)
 
G4double SampleTableThetaCMS (const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
 
G4double GetScatteringAngle (G4int iMomentum, G4int iAngle, G4double position)
 
G4double SampleThetaLab (const G4HadProjectile *aParticle, G4double tmass, G4double A)
 
G4double GetDiffuseElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
 
G4double GetInvElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
 
G4double GetDiffuseElasticSumXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
 
G4double GetInvElasticSumXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
 
G4double IntegralElasticProb (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
 
G4double GetCoulombElasticXsc (const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
 
G4double GetInvCoulombElasticXsc (const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
 
G4double GetCoulombTotalXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z)
 
G4double GetCoulombIntegralXsc (const G4ParticleDefinition *particle, G4double momentum, G4double Z, G4double theta1, G4double theta2)
 
G4double CalculateParticleBeta (const G4ParticleDefinition *particle, G4double momentum)
 
G4double CalculateZommerfeld (G4double beta, G4double Z1, G4double Z2)
 
G4double CalculateAm (G4double momentum, G4double n, G4double Z)
 
G4double CalculateNuclearRad (G4double A)
 
G4double ThetaCMStoThetaLab (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
 
G4double ThetaLabToThetaCMS (const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
 
void TestAngleTable (const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
 
G4double BesselJzero (G4double z)
 
G4double BesselJone (G4double z)
 
G4double DampFactor (G4double z)
 
G4double BesselOneByArg (G4double z)
 
G4double GetDiffElasticProb (G4double theta)
 
G4double GetDiffElasticSumProb (G4double theta)
 
G4double GetDiffElasticSumProbA (G4double alpha)
 
G4double GetIntegrandFunction (G4double theta)
 
G4double GetNuclearRadius ()
 
- Public Member Functions inherited from G4HadronElastic
 G4HadronElastic (const G4String &name="hElasticLHEP")
 
virtual ~G4HadronElastic ()
 
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
void SetLowestEnergyLimit (G4double value)
 
G4double LowestEnergyLimit () const
 
G4double ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual void ModelDescription (std::ostream &) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteractionGetMyPointer () const
 
virtual G4int GetVerboseLevel () const
 
virtual void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair< G4double, G4doubleGetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4doubleGetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 

Private Attributes

G4ParticleDefinitiontheProton
 
G4ParticleDefinitiontheNeutron
 
G4ParticleDefinitiontheDeuteron
 
G4ParticleDefinitiontheAlpha
 
const G4ParticleDefinitionthePionPlus
 
const G4ParticleDefinitionthePionMinus
 
G4double lowEnergyRecoilLimit
 
G4double lowEnergyLimitHE
 
G4double lowEnergyLimitQ
 
G4double lowestEnergyLimit
 
G4double plabLowLimit
 
G4int fEnergyBin
 
G4int fAngleBin
 
G4PhysicsLogVectorfEnergyVector
 
G4PhysicsTablefAngleTable
 
std::vector< G4PhysicsTable * > fAngleBank
 
std::vector< G4doublefElementNumberVector
 
std::vector< G4StringfElementNameVector
 
const G4ParticleDefinitionfParticle
 
G4double fWaveVector
 
G4double fAtomicWeight
 
G4double fAtomicNumber
 
G4double fNuclearRadius
 
G4double fBeta
 
G4double fZommerfeld
 
G4double fAm
 
G4bool fAddCoulomb
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 59 of file G4DiffuseElastic.hh.

Constructor & Destructor Documentation

◆ G4DiffuseElastic()

G4DiffuseElastic::G4DiffuseElastic ( )

Definition at line 74 of file G4DiffuseElastic.cc.

75  : G4HadronElastic("DiffuseElastic"), fParticle(0)
76 {
77  SetMinEnergy( 0.01*MeV ); // 0.01*GeV );
78  SetMaxEnergy( 1.*TeV );
79 
80  verboseLevel = 0;
81  lowEnergyRecoilLimit = 100.*keV;
82  lowEnergyLimitQ = 0.0*GeV;
83  lowEnergyLimitHE = 0.0*GeV;
84  lowestEnergyLimit = 0.0*keV;
85  plabLowLimit = 20.0*MeV;
86 
93 
94  fEnergyBin = 200;
95  fAngleBin = 200;
96 
98 
99  fAngleTable = 0;
100 
101  fParticle = 0;
102  fWaveVector = 0.;
103  fAtomicWeight = 0.;
104  fAtomicNumber = 0.;
105  fNuclearRadius = 0.;
106  fBeta = 0.;
107  fZommerfeld = 0.;
108  fAm = 0.;
109  fAddCoulomb = false;
110 }
static const double MeV
Definition: G4SIunits.hh:211
G4HadronElastic(const G4String &name="hElasticLHEP")
G4ParticleDefinition * theProton
const G4ParticleDefinition * thePionPlus
const G4ParticleDefinition * thePionMinus
void SetMinEnergy(G4double anEnergy)
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static const double GeV
Definition: G4SIunits.hh:214
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4ParticleDefinition * theAlpha
G4PhysicsLogVector * fEnergyVector
G4double lowEnergyRecoilLimit
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
void SetMaxEnergy(const G4double anEnergy)
G4ParticleDefinition * theDeuteron
static const double TeV
Definition: G4SIunits.hh:215
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static const double keV
Definition: G4SIunits.hh:213
const G4ParticleDefinition * fParticle
G4PhysicsTable * fAngleTable
G4ParticleDefinition * theNeutron
Here is the call graph for this function:

◆ ~G4DiffuseElastic()

G4DiffuseElastic::~G4DiffuseElastic ( )
virtual

Definition at line 116 of file G4DiffuseElastic.cc.

117 {
118  if ( fEnergyVector )
119  {
120  delete fEnergyVector;
121  fEnergyVector = 0;
122  }
123  for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
124  it != fAngleBank.end(); ++it )
125  {
126  if ( (*it) ) (*it)->clearAndDestroy();
127 
128  delete *it;
129  *it = 0;
130  }
131  fAngleTable = 0;
132 }
std::vector< G4PhysicsTable * > fAngleBank
G4PhysicsLogVector * fEnergyVector
G4PhysicsTable * fAngleTable

Member Function Documentation

◆ BesselJone()

G4double G4DiffuseElastic::BesselJone ( G4double  z)
inline

Definition at line 330 of file G4DiffuseElastic.hh.

331 {
332  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
333 
334  modvalue = std::fabs(value);
335 
336  if ( modvalue < 8.0 )
337  {
338  value2 = value*value;
339 
340  fact1 = value*(72362614232.0 + value2*(-7895059235.0
341  + value2*( 242396853.1
342  + value2*(-2972611.439
343  + value2*( 15704.48260
344  + value2*(-30.16036606 ) ) ) ) ) );
345 
346  fact2 = 144725228442.0 + value2*(2300535178.0
347  + value2*(18583304.74
348  + value2*(99447.43394
349  + value2*(376.9991397
350  + value2*1.0 ) ) ) );
351  bessel = fact1/fact2;
352  }
353  else
354  {
355  arg = 8.0/modvalue;
356 
357  value2 = arg*arg;
358 
359  shift = modvalue - 2.356194491;
360 
361  fact1 = 1.0 + value2*( 0.183105e-2
362  + value2*(-0.3516396496e-4
363  + value2*(0.2457520174e-5
364  + value2*(-0.240337019e-6 ) ) ) );
365 
366  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
367  + value2*( 0.8449199096e-5
368  + value2*(-0.88228987e-6
369  + value2*0.105787412e-6 ) ) );
370 
371  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
372 
373  if (value < 0.0) bessel = -bessel;
374  }
375  return bessel;
376 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ BesselJzero()

G4double G4DiffuseElastic::BesselJzero ( G4double  z)
inline

Definition at line 278 of file G4DiffuseElastic.hh.

279 {
280  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
281 
282  modvalue = std::fabs(value);
283 
284  if ( value < 8.0 && value > -8.0 )
285  {
286  value2 = value*value;
287 
288  fact1 = 57568490574.0 + value2*(-13362590354.0
289  + value2*( 651619640.7
290  + value2*(-11214424.18
291  + value2*( 77392.33017
292  + value2*(-184.9052456 ) ) ) ) );
293 
294  fact2 = 57568490411.0 + value2*( 1029532985.0
295  + value2*( 9494680.718
296  + value2*(59272.64853
297  + value2*(267.8532712
298  + value2*1.0 ) ) ) );
299 
300  bessel = fact1/fact2;
301  }
302  else
303  {
304  arg = 8.0/modvalue;
305 
306  value2 = arg*arg;
307 
308  shift = modvalue-0.785398164;
309 
310  fact1 = 1.0 + value2*(-0.1098628627e-2
311  + value2*(0.2734510407e-4
312  + value2*(-0.2073370639e-5
313  + value2*0.2093887211e-6 ) ) );
314 
315  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
316  + value2*(-0.6911147651e-5
317  + value2*(0.7621095161e-6
318  - value2*0.934945152e-7 ) ) );
319 
320  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
321  }
322  return bessel;
323 }
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ BesselOneByArg()

G4double G4DiffuseElastic::BesselOneByArg ( G4double  z)
inline

Definition at line 405 of file G4DiffuseElastic.hh.

406 {
407  G4double x2, result;
408 
409  if( std::fabs(x) < 0.01 )
410  {
411  x *= 0.5;
412  x2 = x*x;
413  result = 2. - x2 + x2*x2/6.;
414  }
415  else
416  {
417  result = BesselJone(x)/x;
418  }
419  return result;
420 }
Double_t x2[nxs]
G4double BesselJone(G4double z)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ BuildAngleTable()

void G4DiffuseElastic::BuildAngleTable ( )

Definition at line 999 of file G4DiffuseElastic.cc.

1000 {
1001  G4int i, j;
1002  G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1003  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1004 
1006 
1008 
1009  for( i = 0; i < fEnergyBin; i++)
1010  {
1011  kinE = fEnergyVector->GetLowEdgeEnergy(i);
1012  partMom = std::sqrt( kinE*(kinE + 2*m1) );
1013 
1014  fWaveVector = partMom/hbarc;
1015 
1017  G4double kR2 = kR*kR;
1018  G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1019  G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
1020  // G4double kRlim = 1.2;
1021  // G4double kRlim2 = kRlim*kRlim/kR2;
1022 
1023  alphaMax = kRmax*kRmax/kR2;
1024 
1025 
1026  // if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1027  // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = 15.; // vmg27.11.14
1028 
1029  // if ( alphaMax > 4. || alphaMax < 1. ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg06.01.15
1030 
1031  // G4cout<<"alphaMax = "<<alphaMax<<", ";
1032 
1033  if ( alphaMax >= CLHEP::pi*CLHEP::pi ) alphaMax = CLHEP::pi*CLHEP::pi; // vmg21.10.15
1034 
1035  alphaCoulomb = kRcoul*kRcoul/kR2;
1036 
1037  if( z )
1038  {
1039  a = partMom/m1; // beta*gamma for m1
1040  fBeta = a/std::sqrt(1+a*a);
1042  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1043  }
1044  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1045 
1046  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1047 
1048  G4double delth = alphaMax/fAngleBin;
1049 
1050  sum = 0.;
1051 
1052  // fAddCoulomb = false;
1053  fAddCoulomb = true;
1054 
1055  // for(j = 1; j < fAngleBin; j++)
1056  for(j = fAngleBin-1; j >= 1; j--)
1057  {
1058  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1059  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1060 
1061  // alpha1 = alphaMax*(j-1)/fAngleBin;
1062  // alpha2 = alphaMax*( j )/fAngleBin;
1063 
1064  alpha1 = delth*(j-1);
1065  // if(alpha1 < kRlim2) alpha1 = kRlim2;
1066  alpha2 = alpha1 + delth;
1067 
1068  // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1069  if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
1070 
1071  delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1072  // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1073 
1074  sum += delta;
1075 
1076  angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1077  // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl;
1078  }
1079  fAngleTable->insertAt(i, angleVector);
1080 
1081  // delete[] angleVector;
1082  // delete[] angleBins;
1083  }
1084  return;
1085 }
G4double Legendre10(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)
int G4int
Definition: G4Types.hh:78
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4double GetLowEdgeEnergy(size_t binNumber) const
float hbarc
Definition: hepunit.py:265
static const double pi
Definition: SystemOfUnits.h:53
G4PhysicsLogVector * fEnergyVector
void insertAt(size_t, G4PhysicsVector *)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
G4PhysicsTable * fAngleTable
G4double GetPDGCharge() const
G4double GetIntegrandFunction(G4double theta)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateAm()

G4double G4DiffuseElastic::CalculateAm ( G4double  momentum,
G4double  n,
G4double  Z 
)
inline

Definition at line 451 of file G4DiffuseElastic.hh.

452 {
453  G4double k = momentum/CLHEP::hbarc;
454  G4double ch = 1.13 + 3.76*n*n;
455  G4double zn = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
456  G4double zn2 = zn*zn;
457  fAm = ch/zn2;
458 
459  return fAm;
460 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double hbarc
static const double Bohr_radius
Char_t n[5]
Float_t Z
G4double A13(G4double A) const
Definition: G4Pow.hh:132
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateNuclearRad()

G4double G4DiffuseElastic::CalculateNuclearRad ( G4double  A)
inline

Definition at line 466 of file G4DiffuseElastic.hh.

467 {
468  G4double R, r0, a11, a12, a13, a2, a3;
469 
470  a11 = 1.26; // 1.08, 1.16
471  a12 = 1.; // 1.08, 1.16
472  a13 = 1.12; // 1.08, 1.16
473  a2 = 1.1;
474  a3 = 1.;
475 
476  // Special rms radii for light nucleii
477 
478  if (A < 50.)
479  {
480  if (std::abs(A-1.) < 0.5) return 0.89*CLHEP::fermi; // p
481  else if(std::abs(A-2.) < 0.5) return 2.13*CLHEP::fermi; // d
482  else if( // std::abs(Z-1.) < 0.5 &&
483 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
484 
485  // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
486  else if( // std::abs(Z-2.) < 0.5 &&
487 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
488 
489  else if( // std::abs(Z-3.) < 0.5
490  std::abs(A-7.) < 0.5 ) return 2.40*CLHEP::fermi; // Li7
491  else if( // std::abs(Z-4.) < 0.5
492 std::abs(A-9.) < 0.5) return 2.51*CLHEP::fermi; // Be9
493 
494  else if( 10. < A && A <= 16. ) r0 = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi;
495  else if( 15. < A && A <= 20. ) r0 = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
496  else if( 20. < A && A <= 30. ) r0 = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
497  else r0 = a2*CLHEP::fermi;
498 
499  R = r0*G4Pow::GetInstance()->A13(A);
500  }
501  else
502  {
503  r0 = a3*CLHEP::fermi;
504 
505  R = r0*G4Pow::GetInstance()->powA(A, 0.27);
506  }
507  fNuclearRadius = R;
508  return R;
509  /*
510  G4double r0;
511  if( A < 50. )
512  {
513  if( A > 10. ) r0 = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi; // 1.08CLHEP::fermi;
514  else r0 = 1.1*CLHEP::fermi;
515  fNuclearRadius = r0*G4Pow::GetInstance()->A13(A);
516  }
517  else
518  {
519  r0 = 1.7*CLHEP::fermi; // 1.7*CLHEP::fermi;
520  fNuclearRadius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
521  }
522  return fNuclearRadius;
523  */
524 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double A23(G4double A) const
Definition: G4Pow.hh:160
double A(double temperature)
static const double fermi
Definition: SystemOfUnits.h:82
static const G4double a3
G4double A13(G4double A) const
Definition: G4Pow.hh:132
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
double G4double
Definition: G4Types.hh:76
static const G4double a2
const G4double r0
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateParticleBeta()

G4double G4DiffuseElastic::CalculateParticleBeta ( const G4ParticleDefinition particle,
G4double  momentum 
)
inline

Definition at line 426 of file G4DiffuseElastic.hh.

428 {
429  G4double mass = particle->GetPDGMass();
430  G4double a = momentum/mass;
431  fBeta = a/std::sqrt(1+a*a);
432 
433  return fBeta;
434 }
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ CalculateZommerfeld()

G4double G4DiffuseElastic::CalculateZommerfeld ( G4double  beta,
G4double  Z1,
G4double  Z2 
)
inline

Definition at line 440 of file G4DiffuseElastic.hh.

441 {
443 
444  return fZommerfeld;
445 }
Double_t Z2
static const double fine_structure_const
Double_t Z1
Here is the caller graph for this function:

◆ DampFactor()

G4double G4DiffuseElastic::DampFactor ( G4double  z)
inline

Definition at line 382 of file G4DiffuseElastic.hh.

383 {
384  G4double df;
385  G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
386 
387  // x *= pi;
388 
389  if( std::fabs(x) < 0.01 )
390  {
391  df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
392  }
393  else
394  {
395  df = x/std::sinh(x);
396  }
397  return df;
398 }
Float_t f4
Float_t f2
Float_t f3
double G4double
Definition: G4Types.hh:76
Here is the caller graph for this function:

◆ GetCoulombElasticXsc()

G4double G4DiffuseElastic::GetCoulombElasticXsc ( const G4ParticleDefinition particle,
G4double  theta,
G4double  momentum,
G4double  Z 
)
inline

Definition at line 530 of file G4DiffuseElastic.hh.

534 {
535  G4double sinHalfTheta = std::sin(0.5*theta);
536  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
537  G4double beta = CalculateParticleBeta( particle, momentum);
538  G4double z = particle->GetPDGCharge();
539  G4double n = CalculateZommerfeld( beta, z, Z );
540  G4double am = CalculateAm( momentum, n, Z);
541  G4double k = momentum/CLHEP::hbarc;
542  G4double ch = 0.5*n/k;
543  G4double ch2 = ch*ch;
544  G4double xsc = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
545 
546  return xsc;
547 }
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
static const double hbarc
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Char_t n[5]
Float_t Z
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetCoulombIntegralXsc()

G4double G4DiffuseElastic::GetCoulombIntegralXsc ( const G4ParticleDefinition particle,
G4double  momentum,
G4double  Z,
G4double  theta1,
G4double  theta2 
)
inline

Definition at line 579 of file G4DiffuseElastic.hh.

582 {
583  G4double c1 = std::cos(theta1);
584  G4cout<<"c1 = "<<c1<<G4endl;
585  G4double c2 = std::cos(theta2);
586  G4cout<<"c2 = "<<c2<<G4endl;
587  G4double beta = CalculateParticleBeta( particle, momentum);
588  // G4cout<<"beta = "<<beta<<G4endl;
589  G4double z = particle->GetPDGCharge();
590  G4double n = CalculateZommerfeld( beta, z, Z );
591  // G4cout<<"fZomerfeld = "<<n<<G4endl;
592  G4double am = CalculateAm( momentum, n, Z);
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;
597  G4double ch = n/k;
598  G4double ch2 = ch*ch;
599  am *= 2.;
600  G4double xsc = ch2*CLHEP::twopi*(c1-c2);
601  xsc /= (1 - c1 + am)*(1 - c2 + am);
602 
603  return xsc;
604 }
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
static const double hbarc
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
Float_t Z
TCanvas * c2
Definition: plot_hist.C:75
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static const double twopi
Definition: SystemOfUnits.h:54
G4double GetPDGCharge() const
Here is the call graph for this function:

◆ GetCoulombTotalXsc()

G4double G4DiffuseElastic::GetCoulombTotalXsc ( const G4ParticleDefinition particle,
G4double  momentum,
G4double  Z 
)
inline

Definition at line 554 of file G4DiffuseElastic.hh.

556 {
557  G4double beta = CalculateParticleBeta( particle, momentum);
558  G4cout<<"beta = "<<beta<<G4endl;
559  G4double z = particle->GetPDGCharge();
560  G4double n = CalculateZommerfeld( beta, z, Z );
561  G4cout<<"fZomerfeld = "<<n<<G4endl;
562  G4double am = CalculateAm( momentum, n, Z);
563  G4cout<<"cof Am = "<<am<<G4endl;
564  G4double k = momentum/CLHEP::hbarc;
565  G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
566  G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
567  G4double ch = n/k;
568  G4double ch2 = ch*ch;
569  G4double xsc = ch2*CLHEP::pi/(am +am*am);
570 
571  return xsc;
572 }
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
static const double hbarc
static const double Bohr_radius
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
Float_t Z
static const double fermi
Definition: SystemOfUnits.h:82
static const double pi
Definition: SystemOfUnits.h:53
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
Here is the call graph for this function:

◆ GetDiffElasticProb()

G4double G4DiffuseElastic::GetDiffElasticProb ( G4double  theta)

Definition at line 379 of file G4DiffuseElastic.cc.

384 {
385  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
386  G4double delta, diffuse, gamma;
387  G4double e1, e2, bone, bone2;
388 
389  // G4double wavek = momentum/hbarc; // wave vector
390  // G4double r0 = 1.08*fermi;
391  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
392 
393  if (fParticle == theProton)
394  {
395  diffuse = 0.63*fermi;
396  gamma = 0.3*fermi;
397  delta = 0.1*fermi*fermi;
398  e1 = 0.3*fermi;
399  e2 = 0.35*fermi;
400  }
401  else if (fParticle == theNeutron)
402  {
403  diffuse = 0.63*fermi; // 1.63*fermi; //
404  G4double k0 = 1*GeV/hbarc;
405  diffuse *= k0/fWaveVector;
406 
407  gamma = 0.3*fermi;
408  delta = 0.1*fermi*fermi;
409  e1 = 0.3*fermi;
410  e2 = 0.35*fermi;
411  }
412  else // as proton, if were not defined
413  {
414  diffuse = 0.63*fermi;
415  gamma = 0.3*fermi;
416  delta = 0.1*fermi*fermi;
417  e1 = 0.3*fermi;
418  e2 = 0.35*fermi;
419  }
420  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
421  G4double kr2 = kr*kr;
422  G4double krt = kr*theta;
423 
424  bzero = BesselJzero(krt);
425  bzero2 = bzero*bzero;
426  bone = BesselJone(krt);
427  bone2 = bone*bone;
428  bonebyarg = BesselOneByArg(krt);
429  bonebyarg2 = bonebyarg*bonebyarg;
430 
431  G4double lambda = 15.; // 15 ok
432 
433  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
434 
435  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
436  G4double kgamma2 = kgamma*kgamma;
437 
438  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
439  // G4double dk2t2 = dk2t*dk2t;
440  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
441 
442  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
443 
444  damp = DampFactor(pikdt);
445  damp2 = damp*damp;
446 
447  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
448  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
449 
450 
451  sigma = kgamma2;
452  // sigma += dk2t2;
453  sigma *= bzero2;
454  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
455  sigma += kr2*bonebyarg2;
456  sigma *= damp2; // *rad*rad;
457 
458  return sigma;
459 }
G4double BesselJzero(G4double z)
static const G4double e2
G4ParticleDefinition * theProton
G4double BesselJone(G4double z)
float hbarc
Definition: hepunit.py:265
static const double GeV
Definition: G4SIunits.hh:214
static const G4double e1
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
static const double pi
Definition: G4SIunits.hh:74
G4double BesselOneByArg(G4double z)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
static const double fermi
Definition: G4SIunits.hh:102
G4ParticleDefinition * theNeutron
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetDiffElasticSumProb()

G4double G4DiffuseElastic::GetDiffElasticSumProb ( G4double  theta)

Definition at line 467 of file G4DiffuseElastic.cc.

472 {
473  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
474  G4double delta, diffuse, gamma;
475  G4double e1, e2, bone, bone2;
476 
477  // G4double wavek = momentum/hbarc; // wave vector
478  // G4double r0 = 1.08*fermi;
479  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
480 
481  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
482  G4double kr2 = kr*kr;
483  G4double krt = kr*theta;
484 
485  bzero = BesselJzero(krt);
486  bzero2 = bzero*bzero;
487  bone = BesselJone(krt);
488  bone2 = bone*bone;
489  bonebyarg = BesselOneByArg(krt);
490  bonebyarg2 = bonebyarg*bonebyarg;
491 
492  if (fParticle == theProton)
493  {
494  diffuse = 0.63*fermi;
495  // diffuse = 0.6*fermi;
496  gamma = 0.3*fermi;
497  delta = 0.1*fermi*fermi;
498  e1 = 0.3*fermi;
499  e2 = 0.35*fermi;
500  }
501  else if (fParticle == theNeutron)
502  {
503  diffuse = 0.63*fermi;
504  // diffuse = 0.6*fermi;
505  G4double k0 = 1*GeV/hbarc;
506  diffuse *= k0/fWaveVector;
507  gamma = 0.3*fermi;
508  delta = 0.1*fermi*fermi;
509  e1 = 0.3*fermi;
510  e2 = 0.35*fermi;
511  }
512  else // as proton, if were not defined
513  {
514  diffuse = 0.63*fermi;
515  gamma = 0.3*fermi;
516  delta = 0.1*fermi*fermi;
517  e1 = 0.3*fermi;
518  e2 = 0.35*fermi;
519  }
520  G4double lambda = 15.; // 15 ok
521  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
522  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
523 
524  // G4cout<<"kgamma = "<<kgamma<<G4endl;
525 
526  if(fAddCoulomb) // add Coulomb correction
527  {
528  G4double sinHalfTheta = std::sin(0.5*theta);
529  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
530 
531  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
532  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
533  }
534 
535  G4double kgamma2 = kgamma*kgamma;
536 
537  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
538  // G4cout<<"dk2t = "<<dk2t<<G4endl;
539  // G4double dk2t2 = dk2t*dk2t;
540  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
541 
542  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
543 
544  // G4cout<<"pikdt = "<<pikdt<<G4endl;
545 
546  damp = DampFactor(pikdt);
547  damp2 = damp*damp;
548 
549  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
550  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
551 
552  sigma = kgamma2;
553  // sigma += dk2t2;
554  sigma *= bzero2;
555  sigma += mode2k2*bone2;
556  sigma += e2dk3t*bzero*bone;
557 
558  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
559  sigma += kr2*bonebyarg2; // correction at J1()/()
560 
561  sigma *= damp2; // *rad*rad;
562 
563  return sigma;
564 }
G4double BesselJzero(G4double z)
static const G4double e2
G4ParticleDefinition * theProton
G4double BesselJone(G4double z)
float hbarc
Definition: hepunit.py:265
static const double GeV
Definition: G4SIunits.hh:214
static const G4double e1
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
static const double pi
Definition: G4SIunits.hh:74
G4double BesselOneByArg(G4double z)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
static const double fermi
Definition: G4SIunits.hh:102
G4ParticleDefinition * theNeutron
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetDiffElasticSumProbA()

G4double G4DiffuseElastic::GetDiffElasticSumProbA ( G4double  alpha)

Definition at line 573 of file G4DiffuseElastic.cc.

574 {
575  G4double theta;
576 
577  theta = std::sqrt(alpha);
578 
579  // theta = std::acos( 1 - alpha/2. );
580 
581  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
582  G4double delta, diffuse, gamma;
583  G4double e1, e2, bone, bone2;
584 
585  // G4double wavek = momentum/hbarc; // wave vector
586  // G4double r0 = 1.08*fermi;
587  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
588 
589  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
590  G4double kr2 = kr*kr;
591  G4double krt = kr*theta;
592 
593  bzero = BesselJzero(krt);
594  bzero2 = bzero*bzero;
595  bone = BesselJone(krt);
596  bone2 = bone*bone;
597  bonebyarg = BesselOneByArg(krt);
598  bonebyarg2 = bonebyarg*bonebyarg;
599 
600  if ( fParticle == theProton )
601  {
602  diffuse = 0.63*fermi;
603  // diffuse = 0.6*fermi;
604  gamma = 0.3*fermi;
605  delta = 0.1*fermi*fermi;
606  e1 = 0.3*fermi;
607  e2 = 0.35*fermi;
608  }
609  else if ( fParticle == theNeutron )
610  {
611  diffuse = 0.63*fermi;
612  // diffuse = 0.6*fermi;
613  // G4double k0 = 0.8*GeV/hbarc;
614  // diffuse *= k0/fWaveVector;
615  gamma = 0.3*fermi;
616  delta = 0.1*fermi*fermi;
617  e1 = 0.3*fermi;
618  e2 = 0.35*fermi;
619  }
620  else // as proton, if were not defined
621  {
622  diffuse = 0.63*fermi;
623  gamma = 0.3*fermi;
624  delta = 0.1*fermi*fermi;
625  e1 = 0.3*fermi;
626  e2 = 0.35*fermi;
627  }
628  G4double lambda = 15.; // 15 ok
629  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
630  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
631 
632  // G4cout<<"kgamma = "<<kgamma<<G4endl;
633 
634  if( fAddCoulomb ) // add Coulomb correction
635  {
636  G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
637  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
638 
639  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
640  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
641  }
642  G4double kgamma2 = kgamma*kgamma;
643 
644  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
645  // G4cout<<"dk2t = "<<dk2t<<G4endl;
646  // G4double dk2t2 = dk2t*dk2t;
647  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
648 
649  G4double pikdt = lambda*(1. - G4Exp( -pi*fWaveVector*diffuse*theta/lambda ) ); // wavek*delta;
650 
651  // G4cout<<"pikdt = "<<pikdt<<G4endl;
652 
653  damp = DampFactor( pikdt );
654  damp2 = damp*damp;
655 
656  G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVector*fWaveVector;
657  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
658 
659  sigma = kgamma2;
660  // sigma += dk2t2;
661  sigma *= bzero2;
662  sigma += mode2k2*bone2;
663  sigma += e2dk3t*bzero*bone;
664 
665  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
666  sigma += kr2*bonebyarg2; // correction at J1()/()
667 
668  sigma *= damp2; // *rad*rad;
669 
670  return sigma;
671 }
G4double BesselJzero(G4double z)
static const G4double e2
G4ParticleDefinition * theProton
G4double BesselJone(G4double z)
static const G4double e1
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
static const double pi
Definition: G4SIunits.hh:74
G4double BesselOneByArg(G4double z)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
static const G4double alpha
static const double fermi
Definition: G4SIunits.hh:102
G4ParticleDefinition * theNeutron
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetDiffuseElasticSumXsc()

G4double G4DiffuseElastic::GetDiffuseElasticSumXsc ( const G4ParticleDefinition particle,
G4double  theta,
G4double  momentum,
G4double  A,
G4double  Z 
)

Definition at line 243 of file G4DiffuseElastic.cc.

247 {
248  fParticle = particle;
249  fWaveVector = momentum/hbarc;
250  fAtomicWeight = A;
251  fAtomicNumber = Z;
253  fAddCoulomb = false;
254 
255  G4double z = particle->GetPDGCharge();
256 
257  G4double kRt = fWaveVector*fNuclearRadius*theta;
258  G4double kRtC = 1.9;
259 
260  if( z && (kRt > kRtC) )
261  {
262  fAddCoulomb = true;
263  fBeta = CalculateParticleBeta( particle, momentum);
265  fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
266  }
268 
269  return sigma;
270 }
G4double CalculateNuclearRad(G4double A)
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double GetDiffElasticSumProb(G4double theta)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
double A(double temperature)
float hbarc
Definition: hepunit.py:265
Float_t Z
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetDiffuseElasticXsc()

G4double G4DiffuseElastic::GetDiffuseElasticXsc ( const G4ParticleDefinition particle,
G4double  theta,
G4double  momentum,
G4double  A 
)

Definition at line 172 of file G4DiffuseElastic.cc.

176 {
177  fParticle = particle;
178  fWaveVector = momentum/hbarc;
179  fAtomicWeight = A;
180  fAddCoulomb = false;
182 
184 
185  return sigma;
186 }
G4double CalculateNuclearRad(G4double A)
double A(double temperature)
float hbarc
Definition: hepunit.py:265
G4double GetDiffElasticProb(G4double theta)
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetIntegrandFunction()

G4double G4DiffuseElastic::GetIntegrandFunction ( G4double  theta)

Definition at line 679 of file G4DiffuseElastic.cc.

680 {
681  G4double result;
682 
683  result = GetDiffElasticSumProbA(alpha);
684 
685  // result *= 2*pi*std::sin(theta);
686 
687  return result;
688 }
G4double GetDiffElasticSumProbA(G4double alpha)
double G4double
Definition: G4Types.hh:76
static const G4double alpha
Here is the call graph for this function:
Here is the caller graph for this function:

◆ GetInvCoulombElasticXsc()

G4double G4DiffuseElastic::GetInvCoulombElasticXsc ( const G4ParticleDefinition particle,
G4double  tMand,
G4double  momentum,
G4double  A,
G4double  Z 
)

Definition at line 330 of file G4DiffuseElastic.cc.

334 {
335  G4double m1 = particle->GetPDGMass();
336  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
337 
338  G4int iZ = static_cast<G4int>(Z+0.5);
339  G4int iA = static_cast<G4int>(A+0.5);
340  G4ParticleDefinition * theDef = 0;
341 
342  if (iZ == 1 && iA == 1) theDef = theProton;
343  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
344  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
345  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
346  else if (iZ == 2 && iA == 4) theDef = theAlpha;
347  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
348 
349  G4double tmass = theDef->GetPDGMass();
350 
351  G4LorentzVector lv(0.0,0.0,0.0,tmass);
352  lv += lv1;
353 
354  G4ThreeVector bst = lv.boostVector();
355  lv1.boost(-bst);
356 
357  G4ThreeVector p1 = lv1.vect();
358  G4double ptot = p1.mag();
359  G4double ptot2 = ptot*ptot;
360  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
361 
362  if( cost >= 1.0 ) cost = 1.0;
363  else if( cost <= -1.0) cost = -1.0;
364 
365  G4double thetaCMS = std::acos(cost);
366 
367  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
368 
369  sigma *= pi/ptot2;
370 
371  return sigma;
372 }
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4ParticleDefinition * theProton
int G4int
Definition: G4Types.hh:78
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
double A(double temperature)
G4IonTable * GetIonTable() const
Float_t Z
double mag() const
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4ParticleDefinition * theAlpha
static const double pi
Definition: G4SIunits.hh:74
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * theDeuteron
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
Here is the call graph for this function:

◆ GetInvElasticSumXsc()

G4double G4DiffuseElastic::GetInvElasticSumXsc ( const G4ParticleDefinition particle,
G4double  tMand,
G4double  momentum,
G4double  A,
G4double  Z 
)

Definition at line 278 of file G4DiffuseElastic.cc.

282 {
283  G4double m1 = particle->GetPDGMass();
284 
285  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
286 
287  G4int iZ = static_cast<G4int>(Z+0.5);
288  G4int iA = static_cast<G4int>(A+0.5);
289 
290  G4ParticleDefinition* theDef = 0;
291 
292  if (iZ == 1 && iA == 1) theDef = theProton;
293  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
294  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
295  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
296  else if (iZ == 2 && iA == 4) theDef = theAlpha;
297  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
298 
299  G4double tmass = theDef->GetPDGMass();
300 
301  G4LorentzVector lv(0.0,0.0,0.0,tmass);
302  lv += lv1;
303 
304  G4ThreeVector bst = lv.boostVector();
305  lv1.boost(-bst);
306 
307  G4ThreeVector p1 = lv1.vect();
308  G4double ptot = p1.mag();
309  G4double ptot2 = ptot*ptot;
310  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
311 
312  if( cost >= 1.0 ) cost = 1.0;
313  else if( cost <= -1.0) cost = -1.0;
314 
315  G4double thetaCMS = std::acos(cost);
316 
317  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
318 
319  sigma *= pi/ptot2;
320 
321  return sigma;
322 }
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4ParticleDefinition * theProton
int G4int
Definition: G4Types.hh:78
double A(double temperature)
G4IonTable * GetIonTable() const
Float_t Z
double mag() const
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4ParticleDefinition * theAlpha
static const double pi
Definition: G4SIunits.hh:74
static G4ParticleTable * GetParticleTable()
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4ParticleDefinition * theDeuteron
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
Here is the call graph for this function:

◆ GetInvElasticXsc()

G4double G4DiffuseElastic::GetInvElasticXsc ( const G4ParticleDefinition particle,
G4double  theta,
G4double  momentum,
G4double  A,
G4double  Z 
)

Definition at line 193 of file G4DiffuseElastic.cc.

197 {
198  G4double m1 = particle->GetPDGMass();
199  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
200 
201  G4int iZ = static_cast<G4int>(Z+0.5);
202  G4int iA = static_cast<G4int>(A+0.5);
203  G4ParticleDefinition * theDef = 0;
204 
205  if (iZ == 1 && iA == 1) theDef = theProton;
206  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
207  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
208  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
209  else if (iZ == 2 && iA == 4) theDef = theAlpha;
210  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
211 
212  G4double tmass = theDef->GetPDGMass();
213 
214  G4LorentzVector lv(0.0,0.0,0.0,tmass);
215  lv += lv1;
216 
217  G4ThreeVector bst = lv.boostVector();
218  lv1.boost(-bst);
219 
220  G4ThreeVector p1 = lv1.vect();
221  G4double ptot = p1.mag();
222  G4double ptot2 = ptot*ptot;
223  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
224 
225  if( cost >= 1.0 ) cost = 1.0;
226  else if( cost <= -1.0) cost = -1.0;
227 
228  G4double thetaCMS = std::acos(cost);
229 
230  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
231 
232  sigma *= pi/ptot2;
233 
234  return sigma;
235 }
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
G4ParticleDefinition * theProton
int G4int
Definition: G4Types.hh:78
double A(double temperature)
G4IonTable * GetIonTable() const
Float_t Z
double mag() const
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4ParticleDefinition * theAlpha
static const double pi
Definition: G4SIunits.hh:74
static G4ParticleTable * GetParticleTable()
G4ParticleDefinition * theDeuteron
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
Here is the call graph for this function:

◆ GetNuclearRadius()

G4double G4DiffuseElastic::GetNuclearRadius ( )
inline

Definition at line 192 of file G4DiffuseElastic.hh.

192 {return fNuclearRadius;};

◆ GetScatteringAngle()

G4double G4DiffuseElastic::GetScatteringAngle ( G4int  iMomentum,
G4int  iAngle,
G4double  position 
)

Definition at line 1092 of file G4DiffuseElastic.cc.

1093 {
1094  G4double x1, x2, y1, y2, randAngle;
1095 
1096  if( iAngle == 0 )
1097  {
1098  randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1099  // iAngle++;
1100  }
1101  else
1102  {
1103  if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1104  {
1105  iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1106  }
1107  y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1108  y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1109 
1110  x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1111  x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1112 
1113  if ( x1 == x2 ) randAngle = x2;
1114  else
1115  {
1116  if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1117  else
1118  {
1119  randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1120  }
1121  }
1122  }
1123  return randAngle;
1124 }
Double_t y2[nxs]
Double_t y1[nxs]
Double_t x2[nxs]
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
Double_t x1[nxs]
double G4double
Definition: G4Types.hh:76
G4PhysicsTable * fAngleTable
Here is the caller graph for this function:

◆ Initialise()

void G4DiffuseElastic::Initialise ( )

Definition at line 138 of file G4DiffuseElastic.cc.

139 {
140 
141  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
142 
143  const G4ElementTable* theElementTable = G4Element::GetElementTable();
144 
145  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
146 
147  for( jEl = 0; jEl < numOfEl; ++jEl) // application element loop
148  {
149  fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
152 
153  if( verboseLevel > 0 )
154  {
155  G4cout<<"G4DiffuseElastic::Initialise() the element: "
156  <<(*theElementTable)[jEl]->GetName()<<G4endl;
157  }
159  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
160 
161  BuildAngleTable();
162  fAngleBank.push_back(fAngleTable);
163  }
164  return;
165 }
G4double CalculateNuclearRad(G4double A)
std::vector< G4String > fElementNameVector
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
std::vector< G4PhysicsTable * > fAngleBank
#define G4endl
Definition: G4ios.hh:61
std::vector< G4Element * > G4ElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
std::vector< G4double > fElementNumberVector
G4double GetAtomicMassAmu(const G4String &symb) const
G4PhysicsTable * fAngleTable
Here is the call graph for this function:

◆ InitialiseOnFly()

void G4DiffuseElastic::InitialiseOnFly ( G4double  Z,
G4double  A 
)

Definition at line 973 of file G4DiffuseElastic.cc.

974 {
975  fAtomicNumber = Z; // atomic number
976  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
977 
979 
980  if( verboseLevel > 0 )
981  {
982  G4cout<<"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
983  <<Z<<"; and A = "<<A<<G4endl;
984  }
986 
987  BuildAngleTable();
988 
989  fAngleBank.push_back(fAngleTable);
990 
991  return;
992 }
G4double CalculateNuclearRad(G4double A)
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
Float_t Z
std::vector< G4PhysicsTable * > fAngleBank
#define G4endl
Definition: G4ios.hh:61
std::vector< G4double > fElementNumberVector
G4double GetAtomicMassAmu(const G4String &symb) const
G4PhysicsTable * fAngleTable
Here is the call graph for this function:
Here is the caller graph for this function:

◆ IntegralElasticProb()

G4double G4DiffuseElastic::IntegralElasticProb ( const G4ParticleDefinition particle,
G4double  theta,
G4double  momentum,
G4double  A 
)

Definition at line 695 of file G4DiffuseElastic.cc.

699 {
700  G4double result;
701  fParticle = particle;
702  fWaveVector = momentum/hbarc;
703  fAtomicWeight = A;
704 
706 
707 
709 
710  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
711  result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
712 
713  return result;
714 }
G4double CalculateNuclearRad(G4double A)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
double A(double temperature)
float hbarc
Definition: hepunit.py:265
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
G4double GetIntegrandFunction(G4double theta)
Here is the call graph for this function:

◆ IsApplicable()

G4bool G4DiffuseElastic::IsApplicable ( const G4HadProjectile projectile,
G4Nucleus nucleus 
)
inlinevirtual

Reimplemented from G4HadronicInteraction.

Definition at line 233 of file G4DiffuseElastic.hh.

235 {
236  if( ( projectile.GetDefinition() == G4Proton::Proton() ||
237  projectile.GetDefinition() == G4Neutron::Neutron() ||
238  projectile.GetDefinition() == G4PionPlus::PionPlus() ||
239  projectile.GetDefinition() == G4PionMinus::PionMinus() ||
240  projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
241  projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
242 
243  nucleus.GetZ_asInt() >= 2 ) return true;
244  else return false;
245 }
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
const G4ParticleDefinition * GetDefinition() const
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
Here is the call graph for this function:

◆ NeutronTuniform()

G4double G4DiffuseElastic::NeutronTuniform ( G4int  Z)

Definition at line 824 of file G4DiffuseElastic.cc.

825 {
826  G4double elZ = G4double(Z);
827  elZ -= 1.;
828  // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.;
829  G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.;
830  return Tkin;
831 }
Float_t Z
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleInvariantT()

G4double G4DiffuseElastic::SampleInvariantT ( const G4ParticleDefinition p,
G4double  plab,
G4int  Z,
G4int  A 
)
virtual

Reimplemented from G4HadronElastic.

Definition at line 786 of file G4DiffuseElastic.cc.

788 {
789  fParticle = aParticle;
790  G4double m1 = fParticle->GetPDGMass(), t;
791  G4double totElab = std::sqrt(m1*m1+p*p);
793  G4LorentzVector lv1(p,0.0,0.0,totElab);
794  G4LorentzVector lv(0.0,0.0,0.0,mass2);
795  lv += lv1;
796 
797  G4ThreeVector bst = lv.boostVector();
798  lv1.boost(-bst);
799 
800  G4ThreeVector p1 = lv1.vect();
801  G4double momentumCMS = p1.mag();
802 
803  if( aParticle == theNeutron)
804  {
805  G4double Tmax = NeutronTuniform( Z );
806  G4double pCMS2 = momentumCMS*momentumCMS;
807  G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1;
808 
809  if( Tkin <= Tmax )
810  {
811  t = 4.*pCMS2*G4UniformRand();
812  // G4cout<<Tkin<<", "<<Tmax<<", "<<std::sqrt(t)<<"; ";
813  return t;
814  }
815  }
816 
817  t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
818 
819  return t;
820 }
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double NeutronTuniform(G4int Z)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
Float_t Z
double mag() const
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
G4ParticleDefinition * theNeutron
Here is the call graph for this function:

◆ SampleT()

G4double G4DiffuseElastic::SampleT ( const G4ParticleDefinition aParticle,
G4double  p,
G4double  A 
)

Definition at line 720 of file G4DiffuseElastic.cc.

721 {
722  G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
723  G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
724  return t;
725 }
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
double A(double temperature)
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleTableT()

G4double G4DiffuseElastic::SampleTableT ( const G4ParticleDefinition aParticle,
G4double  p,
G4double  Z,
G4double  A 
)

Definition at line 838 of file G4DiffuseElastic.cc.

840 {
841  G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
842  G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
843  // G4double t = p*p*alpha; // -t !!!
844  return t;
845 }
double A(double temperature)
Float_t Z
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
double G4double
Definition: G4Types.hh:76
static const G4double alpha
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleTableThetaCMS()

G4double G4DiffuseElastic::SampleTableThetaCMS ( const G4ParticleDefinition aParticle,
G4double  p,
G4double  Z,
G4double  A 
)

Definition at line 853 of file G4DiffuseElastic.cc.

855 {
856  size_t iElement;
857  G4int iMomentum, iAngle;
858  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
859  G4double m1 = particle->GetPDGMass();
860 
861  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
862  {
863  if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
864  }
865  if ( iElement == fElementNumberVector.size() )
866  {
867  InitialiseOnFly(Z,A); // table preparation, if needed
868 
869  // iElement--;
870 
871  // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
872  // << " is not found, return zero angle" << G4endl;
873  // return 0.; // no table for this element
874  }
875  // G4cout<<"iElement = "<<iElement<<G4endl;
876 
877  fAngleTable = fAngleBank[iElement];
878 
879  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
880 
881  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
882  {
883  if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
884  }
885  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
886  if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
887 
888  // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
889 
890  if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
891  {
892  position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
893 
894  // G4cout<<"position = "<<position<<G4endl;
895 
896  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
897  {
898  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
899  }
900  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
901 
902  // G4cout<<"iAngle = "<<iAngle<<G4endl;
903 
904  randAngle = GetScatteringAngle(iMomentum, iAngle, position);
905 
906  // G4cout<<"randAngle = "<<randAngle<<G4endl;
907  }
908  else // kinE inside between energy table edges
909  {
910  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
911  position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
912 
913  // G4cout<<"position = "<<position<<G4endl;
914 
915  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
916  {
917  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
918  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
919  }
920  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
921 
922  // G4cout<<"iAngle = "<<iAngle<<G4endl;
923 
924  theta2 = GetScatteringAngle(iMomentum, iAngle, position);
925 
926  // G4cout<<"theta2 = "<<theta2<<G4endl;
927  E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
928 
929  // G4cout<<"E2 = "<<E2<<G4endl;
930 
931  iMomentum--;
932 
933  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
934 
935  // G4cout<<"position = "<<position<<G4endl;
936 
937  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
938  {
939  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
940  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
941  }
942  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
943 
944  theta1 = GetScatteringAngle(iMomentum, iAngle, position);
945 
946  // G4cout<<"theta1 = "<<theta1<<G4endl;
947 
948  E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
949 
950  // G4cout<<"E1 = "<<E1<<G4endl;
951 
952  W = 1.0/(E2 - E1);
953  W1 = (E2 - kinE)*W;
954  W2 = (kinE - E1)*W;
955 
956  randAngle = W1*theta1 + W2*theta2;
957 
958  // randAngle = theta2;
959  // G4cout<<"randAngle = "<<randAngle<<G4endl;
960  }
961  // G4double angle = randAngle;
962  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
963 
964  if(randAngle < 0.) randAngle = 0.;
965 
966  return randAngle;
967 }
void InitialiseOnFly(G4double Z, G4double A)
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
int G4int
Definition: G4Types.hh:78
G4double GetLowEdgeEnergy(size_t binNumber) const
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
Float_t Z
std::vector< G4PhysicsTable * > fAngleBank
G4PhysicsLogVector * fEnergyVector
double G4double
Definition: G4Types.hh:76
std::vector< G4double > fElementNumberVector
G4PhysicsTable * fAngleTable
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleThetaCMS()

G4double G4DiffuseElastic::SampleThetaCMS ( const G4ParticleDefinition aParticle,
G4double  p,
G4double  A 
)

Definition at line 733 of file G4DiffuseElastic.cc.

735 {
736  G4int i, iMax = 100;
737  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
738 
739  fParticle = particle;
740  fWaveVector = momentum/hbarc;
741  fAtomicWeight = A;
742 
744 
745  thetaMax = 10.174/fWaveVector/fNuclearRadius;
746 
747  if (thetaMax > pi) thetaMax = pi;
748 
750 
751  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
752  norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
753 
754  norm *= G4UniformRand();
755 
756  for(i = 1; i <= iMax; i++)
757  {
758  theta1 = (i-1)*thetaMax/iMax;
759  theta2 = i*thetaMax/iMax;
760  sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
761 
762  if ( sum >= norm )
763  {
764  result = 0.5*(theta1 + theta2);
765  break;
766  }
767  }
768  if (i > iMax ) result = 0.5*(theta1 + theta2);
769 
770  G4double sigma = pi*thetaMax/iMax;
771 
772  result += G4RandGauss::shoot(0.,sigma);
773 
774  if(result < 0.) result = 0.;
775  if(result > thetaMax) result = thetaMax;
776 
777  return result;
778 }
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double CalculateNuclearRad(G4double A)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
Float_t norm
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
float hbarc
Definition: hepunit.py:265
static const double pi
Definition: G4SIunits.hh:74
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
G4double GetIntegrandFunction(G4double theta)
Here is the call graph for this function:
Here is the caller graph for this function:

◆ SampleThetaLab()

G4double G4DiffuseElastic::SampleThetaLab ( const G4HadProjectile aParticle,
G4double  tmass,
G4double  A 
)

Definition at line 1135 of file G4DiffuseElastic.cc.

1137 {
1138  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1139  G4double m1 = theParticle->GetPDGMass();
1140  G4double plab = aParticle->GetTotalMomentum();
1141  G4LorentzVector lv1 = aParticle->Get4Momentum();
1142  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1143  lv += lv1;
1144 
1145  G4ThreeVector bst = lv.boostVector();
1146  lv1.boost(-bst);
1147 
1148  G4ThreeVector p1 = lv1.vect();
1149  G4double ptot = p1.mag();
1150  G4double tmax = 4.0*ptot*ptot;
1151  G4double t = 0.0;
1152 
1153 
1154  //
1155  // Sample t
1156  //
1157 
1158  t = SampleT( theParticle, ptot, A);
1159 
1160  // NaN finder
1161  if(!(t < 0.0 || t >= 0.0))
1162  {
1163  if (verboseLevel > 0)
1164  {
1165  G4cout << "G4DiffuseElastic:WARNING: A = " << A
1166  << " mom(GeV)= " << plab/GeV
1167  << " S-wave will be sampled"
1168  << G4endl;
1169  }
1170  t = G4UniformRand()*tmax;
1171  }
1172  if(verboseLevel>1)
1173  {
1174  G4cout <<" t= " << t << " tmax= " << tmax
1175  << " ptot= " << ptot << G4endl;
1176  }
1177  // Sampling of angles in CM system
1178 
1179  G4double phi = G4UniformRand()*twopi;
1180  G4double cost = 1. - 2.0*t/tmax;
1181  G4double sint;
1182 
1183  if( cost >= 1.0 )
1184  {
1185  cost = 1.0;
1186  sint = 0.0;
1187  }
1188  else if( cost <= -1.0)
1189  {
1190  cost = -1.0;
1191  sint = 0.0;
1192  }
1193  else
1194  {
1195  sint = std::sqrt((1.0-cost)*(1.0+cost));
1196  }
1197  if (verboseLevel>1)
1198  {
1199  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1200  }
1201  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1202  v1 *= ptot;
1203  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1204 
1205  nlv1.boost(bst);
1206 
1207  G4ThreeVector np1 = nlv1.vect();
1208 
1209  // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1210 
1211  G4double theta = np1.theta();
1212 
1213  return theta;
1214 }
const G4LorentzVector & Get4Momentum() const
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
double theta() const
double mag() const
HepLorentzVector & boost(double, double, double)
static const double twopi
Definition: G4SIunits.hh:75
static const double GeV
Definition: G4SIunits.hh:214
G4double GetTotalMomentum() const
const G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ SetHEModelLowLimit()

void G4DiffuseElastic::SetHEModelLowLimit ( G4double  value)
inline

Definition at line 257 of file G4DiffuseElastic.hh.

258 {
259  lowEnergyLimitHE = value;
260 }

◆ SetLowestEnergyLimit()

void G4DiffuseElastic::SetLowestEnergyLimit ( G4double  value)
inline

Definition at line 267 of file G4DiffuseElastic.hh.

268 {
269  lowestEnergyLimit = value;
270 }

◆ SetPlabLowLimit()

void G4DiffuseElastic::SetPlabLowLimit ( G4double  value)
inline

Definition at line 252 of file G4DiffuseElastic.hh.

253 {
254  plabLowLimit = value;
255 }

◆ SetQModelLowLimit()

void G4DiffuseElastic::SetQModelLowLimit ( G4double  value)
inline

Definition at line 262 of file G4DiffuseElastic.hh.

263 {
264  lowEnergyLimitQ = value;
265 }

◆ SetRecoilKinEnergyLimit()

void G4DiffuseElastic::SetRecoilKinEnergyLimit ( G4double  value)
inline

Definition at line 247 of file G4DiffuseElastic.hh.

248 {
249  lowEnergyRecoilLimit = value;
250 }
G4double lowEnergyRecoilLimit

◆ TestAngleTable()

void G4DiffuseElastic::TestAngleTable ( const G4ParticleDefinition theParticle,
G4double  partMom,
G4double  Z,
G4double  A 
)

Definition at line 1342 of file G4DiffuseElastic.cc.

1344 {
1345  fAtomicNumber = Z; // atomic number
1346  fAtomicWeight = A; // number of nucleons
1348 
1349 
1350 
1351  G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1352  <<Z<<"; and A = "<<A<<G4endl;
1353 
1355 
1356 
1357 
1358 
1359  G4int i=0, j;
1360  G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
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.;
1364  G4double epsilon = 0.001;
1365 
1367 
1369 
1370  fWaveVector = partMom/hbarc;
1371 
1373  G4double kR2 = kR*kR;
1374  G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1375  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1376 
1377  alphaMax = kRmax*kRmax/kR2;
1378 
1379  if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1380 
1381  alphaCoulomb = kRcoul*kRcoul/kR2;
1382 
1383  if( z )
1384  {
1385  a = partMom/m1; // beta*gamma for m1
1386  fBeta = a/std::sqrt(1+a*a);
1388  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1389  }
1390  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1391 
1392  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1393 
1394 
1395  fAddCoulomb = false;
1396 
1397  for(j = 1; j < fAngleBin; j++)
1398  {
1399  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1400  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1401 
1402  alpha1 = alphaMax*(j-1)/fAngleBin;
1403  alpha2 = alphaMax*( j )/fAngleBin;
1404 
1405  if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1406 
1407  deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1408  deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1409  deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1410  alpha1, alpha2,epsilon);
1411 
1412  // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1413  // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1414 
1415  sumL10 += deltaL10;
1416  sumL96 += deltaL96;
1417  sumAG += deltaAG;
1418 
1419  G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1420  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1421 
1422  angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1423  }
1424  fAngleTable->insertAt(i,angleVector);
1425  fAngleBank.push_back(fAngleTable);
1426 
1427  /*
1428  // Integral over all angle range - Bad accuracy !!!
1429 
1430  sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1431  sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1432  sumAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1433  0., alpha2,epsilon);
1434  G4cout<<G4endl;
1435  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1436  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1437  */
1438  return;
1439 }
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
G4double CalculateNuclearRad(G4double A)
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)
int G4int
Definition: G4Types.hh:78
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
float hbarc
Definition: hepunit.py:265
Float_t Z
std::vector< G4PhysicsTable * > fAngleBank
static const double degree
Definition: G4SIunits.hh:143
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
const G4ParticleDefinition * fParticle
double G4double
Definition: G4Types.hh:76
std::vector< G4double > fElementNumberVector
G4PhysicsTable * fAngleTable
G4double GetPDGCharge() const
double epsilon(double density, double temperature)
G4double GetIntegrandFunction(G4double theta)
Here is the call graph for this function:

◆ ThetaCMStoThetaLab()

G4double G4DiffuseElastic::ThetaCMStoThetaLab ( const G4DynamicParticle aParticle,
G4double  tmass,
G4double  thetaCMS 
)

Definition at line 1223 of file G4DiffuseElastic.cc.

1225 {
1226  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1227  G4double m1 = theParticle->GetPDGMass();
1228  // G4double plab = aParticle->GetTotalMomentum();
1229  G4LorentzVector lv1 = aParticle->Get4Momentum();
1230  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1231 
1232  lv += lv1;
1233 
1234  G4ThreeVector bst = lv.boostVector();
1235 
1236  lv1.boost(-bst);
1237 
1238  G4ThreeVector p1 = lv1.vect();
1239  G4double ptot = p1.mag();
1240 
1241  G4double phi = G4UniformRand()*twopi;
1242  G4double cost = std::cos(thetaCMS);
1243  G4double sint;
1244 
1245  if( cost >= 1.0 )
1246  {
1247  cost = 1.0;
1248  sint = 0.0;
1249  }
1250  else if( cost <= -1.0)
1251  {
1252  cost = -1.0;
1253  sint = 0.0;
1254  }
1255  else
1256  {
1257  sint = std::sqrt((1.0-cost)*(1.0+cost));
1258  }
1259  if (verboseLevel>1)
1260  {
1261  G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1262  }
1263  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1264  v1 *= ptot;
1265  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1266 
1267  nlv1.boost(bst);
1268 
1269  G4ThreeVector np1 = nlv1.vect();
1270 
1271 
1272  G4double thetaLab = np1.theta();
1273 
1274  return thetaLab;
1275 }
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double theta() const
double mag() const
HepLorentzVector & boost(double, double, double)
static const double twopi
Definition: G4SIunits.hh:75
G4LorentzVector Get4Momentum() const
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

◆ ThetaLabToThetaCMS()

G4double G4DiffuseElastic::ThetaLabToThetaCMS ( const G4DynamicParticle aParticle,
G4double  tmass,
G4double  thetaLab 
)

Definition at line 1283 of file G4DiffuseElastic.cc.

1285 {
1286  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1287  G4double m1 = theParticle->GetPDGMass();
1288  G4double plab = aParticle->GetTotalMomentum();
1289  G4LorentzVector lv1 = aParticle->Get4Momentum();
1290  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1291 
1292  lv += lv1;
1293 
1294  G4ThreeVector bst = lv.boostVector();
1295 
1296  // lv1.boost(-bst);
1297 
1298  // G4ThreeVector p1 = lv1.vect();
1299  // G4double ptot = p1.mag();
1300 
1301  G4double phi = G4UniformRand()*twopi;
1302  G4double cost = std::cos(thetaLab);
1303  G4double sint;
1304 
1305  if( cost >= 1.0 )
1306  {
1307  cost = 1.0;
1308  sint = 0.0;
1309  }
1310  else if( cost <= -1.0)
1311  {
1312  cost = -1.0;
1313  sint = 0.0;
1314  }
1315  else
1316  {
1317  sint = std::sqrt((1.0-cost)*(1.0+cost));
1318  }
1319  if (verboseLevel>1)
1320  {
1321  G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1322  }
1323  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1324  v1 *= plab;
1325  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1326 
1327  nlv1.boost(-bst);
1328 
1329  G4ThreeVector np1 = nlv1.vect();
1330 
1331 
1332  G4double thetaCMS = np1.theta();
1333 
1334  return thetaCMS;
1335 }
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double theta() const
HepLorentzVector & boost(double, double, double)
static const double twopi
Definition: G4SIunits.hh:75
G4LorentzVector Get4Momentum() const
G4ParticleDefinition * GetDefinition() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
Here is the call graph for this function:

Member Data Documentation

◆ fAddCoulomb

G4bool G4DiffuseElastic::fAddCoulomb
private

Definition at line 229 of file G4DiffuseElastic.hh.

◆ fAm

G4double G4DiffuseElastic::fAm
private

Definition at line 228 of file G4DiffuseElastic.hh.

◆ fAngleBank

std::vector<G4PhysicsTable*> G4DiffuseElastic::fAngleBank
private

Definition at line 216 of file G4DiffuseElastic.hh.

◆ fAngleBin

G4int G4DiffuseElastic::fAngleBin
private

Definition at line 212 of file G4DiffuseElastic.hh.

◆ fAngleTable

G4PhysicsTable* G4DiffuseElastic::fAngleTable
private

Definition at line 215 of file G4DiffuseElastic.hh.

◆ fAtomicNumber

G4double G4DiffuseElastic::fAtomicNumber
private

Definition at line 224 of file G4DiffuseElastic.hh.

◆ fAtomicWeight

G4double G4DiffuseElastic::fAtomicWeight
private

Definition at line 223 of file G4DiffuseElastic.hh.

◆ fBeta

G4double G4DiffuseElastic::fBeta
private

Definition at line 226 of file G4DiffuseElastic.hh.

◆ fElementNameVector

std::vector<G4String> G4DiffuseElastic::fElementNameVector
private

Definition at line 219 of file G4DiffuseElastic.hh.

◆ fElementNumberVector

std::vector<G4double> G4DiffuseElastic::fElementNumberVector
private

Definition at line 218 of file G4DiffuseElastic.hh.

◆ fEnergyBin

G4int G4DiffuseElastic::fEnergyBin
private

Definition at line 211 of file G4DiffuseElastic.hh.

◆ fEnergyVector

G4PhysicsLogVector* G4DiffuseElastic::fEnergyVector
private

Definition at line 214 of file G4DiffuseElastic.hh.

◆ fNuclearRadius

G4double G4DiffuseElastic::fNuclearRadius
private

Definition at line 225 of file G4DiffuseElastic.hh.

◆ fParticle

const G4ParticleDefinition* G4DiffuseElastic::fParticle
private

Definition at line 221 of file G4DiffuseElastic.hh.

◆ fWaveVector

G4double G4DiffuseElastic::fWaveVector
private

Definition at line 222 of file G4DiffuseElastic.hh.

◆ fZommerfeld

G4double G4DiffuseElastic::fZommerfeld
private

Definition at line 227 of file G4DiffuseElastic.hh.

◆ lowEnergyLimitHE

G4double G4DiffuseElastic::lowEnergyLimitHE
private

Definition at line 206 of file G4DiffuseElastic.hh.

◆ lowEnergyLimitQ

G4double G4DiffuseElastic::lowEnergyLimitQ
private

Definition at line 207 of file G4DiffuseElastic.hh.

◆ lowEnergyRecoilLimit

G4double G4DiffuseElastic::lowEnergyRecoilLimit
private

Definition at line 205 of file G4DiffuseElastic.hh.

◆ lowestEnergyLimit

G4double G4DiffuseElastic::lowestEnergyLimit
private

Definition at line 208 of file G4DiffuseElastic.hh.

◆ plabLowLimit

G4double G4DiffuseElastic::plabLowLimit
private

Definition at line 209 of file G4DiffuseElastic.hh.

◆ theAlpha

G4ParticleDefinition* G4DiffuseElastic::theAlpha
private

Definition at line 200 of file G4DiffuseElastic.hh.

◆ theDeuteron

G4ParticleDefinition* G4DiffuseElastic::theDeuteron
private

Definition at line 199 of file G4DiffuseElastic.hh.

◆ theNeutron

G4ParticleDefinition* G4DiffuseElastic::theNeutron
private

Definition at line 198 of file G4DiffuseElastic.hh.

◆ thePionMinus

const G4ParticleDefinition* G4DiffuseElastic::thePionMinus
private

Definition at line 203 of file G4DiffuseElastic.hh.

◆ thePionPlus

const G4ParticleDefinition* G4DiffuseElastic::thePionPlus
private

Definition at line 202 of file G4DiffuseElastic.hh.

◆ theProton

G4ParticleDefinition* G4DiffuseElastic::theProton
private

Definition at line 192 of file G4DiffuseElastic.hh.


The documentation for this class was generated from the following files: