Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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)
 
G4int GetVerboseLevel () const
 
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
 
virtual const std::pair
< G4double, G4double
GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double
GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 

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 ( )

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 
87  theProton = G4Proton::Proton();
88  theNeutron = G4Neutron::Neutron();
89  theDeuteron = G4Deuteron::Deuteron();
90  theAlpha = G4Alpha::Alpha();
91  thePionPlus = G4PionPlus::PionPlus();
92  thePionMinus = G4PionMinus::PionMinus();
93 
94  fEnergyBin = 200;
95  fAngleBin = 200;
96 
97  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
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 }
G4HadronElastic(const G4String &name="hElasticLHEP")
static constexpr double TeV
Definition: G4SIunits.hh:218
void SetMinEnergy(G4double anEnergy)
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static constexpr double GeV
Definition: G4SIunits.hh:217
void SetMaxEnergy(const G4double anEnergy)
static constexpr double MeV
Definition: G4SIunits.hh:214
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

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 }

Member Function Documentation

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 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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 }
G4double G4ParticleHPJENDLHEData::G4double result
tuple x
Definition: test.py:50
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:

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 
1007  fAngleTable = new G4PhysicsTable( fEnergyBin );
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 
1016  G4double kR = fWaveVector*fNuclearRadius;
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);
1041  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
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)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4double GetPDGMass() const
tuple z
Definition: test.py:28
void insertAt(size_t, G4PhysicsVector *)
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
static constexpr double pi
Definition: SystemOfUnits.h:54
G4double GetIntegrandFunction(G4double theta)

Here is the call graph for this function:

Here is the caller graph for this function:

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 constexpr double Bohr_radius
static constexpr double hbarc
const G4int n
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:

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 powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4double A23(G4double A) const
Definition: G4Pow.hh:160
double A(double temperature)
G4double A13(G4double A) const
Definition: G4Pow.hh:132
static constexpr double fermi
Definition: SystemOfUnits.h:83
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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 }
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4double GetPDGMass() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 440 of file G4DiffuseElastic.hh.

441 {
442  fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
443 
444  return fZommerfeld;
445 }
static constexpr double fine_structure_const

Here is the caller graph for this function:

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 }
tuple x
Definition: test.py:50
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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 constexpr double hbarc
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
const G4int n
tuple z
Definition: test.py:28
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:

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 }
static c2_factory< G4double > c2
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
static constexpr double hbarc
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4GLOB_DLL std::ostream G4cout
const G4int n
tuple z
Definition: test.py:28
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
static constexpr double twopi
Definition: SystemOfUnits.h:55
tuple c1
Definition: plottest35.py:14

Here is the call graph for this function:

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 constexpr double Bohr_radius
static constexpr double hbarc
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4GLOB_DLL std::ostream G4cout
const G4int n
static constexpr double fermi
Definition: SystemOfUnits.h:83
tuple z
Definition: test.py:28
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
static constexpr double pi
Definition: SystemOfUnits.h:54

Here is the call graph for this function:

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)
G4double BesselJone(G4double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double BesselOneByArg(G4double z)
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103

Here is the call graph for this function:

Here is the caller graph for this function:

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)
G4double BesselJone(G4double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double BesselOneByArg(G4double z)
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103

Here is the call graph for this function:

Here is the caller graph for this function:

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)
G4double BesselJone(G4double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double DampFactor(G4double z)
G4double BesselOneByArg(G4double z)
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double fermi
Definition: G4SIunits.hh:103
static const G4double alpha

Here is the call graph for this function:

Here is the caller graph for this function:

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;
252  fNuclearRadius = CalculateNuclearRad(A);
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);
264  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
265  fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
266  }
267  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
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)
tuple z
Definition: test.py:28
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:

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;
181  fNuclearRadius = CalculateNuclearRad(A);
182 
183  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
184 
185  return sigma;
186 }
G4double CalculateNuclearRad(G4double A)
double A(double temperature)
G4double GetDiffElasticProb(G4double theta)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4DiffuseElastic::GetIntegrandFunction ( G4double  theta)

Definition at line 679 of file G4DiffuseElastic.cc.

680 {
682 
683  result = GetDiffElasticSumProbA(alpha);
684 
685  // result *= 2*pi*std::sin(theta);
686 
687  return result;
688 }
G4double G4ParticleHPJENDLHEData::G4double result
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:

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:503
int G4int
Definition: G4Types.hh:78
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
G4IonTable * GetIonTable() const
double A(double temperature)
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const

Here is the call graph for this function:

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:503
int G4int
Definition: G4Types.hh:78
G4IonTable * GetIonTable() const
double A(double temperature)
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const

Here is the call graph for this function:

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:503
int G4int
Definition: G4Types.hh:78
G4IonTable * GetIonTable() const
double A(double temperature)
static G4Triton * Triton()
Definition: G4Triton.cc:95
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const

Here is the call graph for this function:

G4double G4DiffuseElastic::GetNuclearRadius ( )
inline

Definition at line 192 of file G4DiffuseElastic.hh.

192 {return fNuclearRadius;};
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 }
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

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
150  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( fAtomicNumber ) );
151  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
152 
153  if( verboseLevel > 0 )
154  {
155  G4cout<<"G4DiffuseElastic::Initialise() the element: "
156  <<(*theElementTable)[jEl]->GetName()<<G4endl;
157  }
158  fElementNumberVector.push_back(fAtomicNumber);
159  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
160 
161  BuildAngleTable();
162  fAngleBank.push_back(fAngleTable);
163  }
164  return;
165 }
G4double CalculateNuclearRad(G4double A)
static G4NistManager * Instance()
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61
std::vector< G4Element * > G4ElementTable
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398

Here is the call graph for this function:

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 
978  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
979 
980  if( verboseLevel > 0 )
981  {
982  G4cout<<"G4DiffuseElastic::InitialiseOnFly() the element with Z = "
983  <<Z<<"; and A = "<<A<<G4endl;
984  }
985  fElementNumberVector.push_back(fAtomicNumber);
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)
G4double GetAtomicMassAmu(const G4String &symb) const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 695 of file G4DiffuseElastic.cc.

699 {
701  fParticle = particle;
702  fWaveVector = momentum/hbarc;
703  fAtomicWeight = A;
704 
705  fNuclearRadius = CalculateNuclearRad(A);
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 G4ParticleHPJENDLHEData::G4double result
G4double CalculateNuclearRad(G4double A)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
double A(double temperature)
double G4double
Definition: G4Types.hh:76
G4double GetIntegrandFunction(G4double theta)

Here is the call graph for this function:

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
const G4ParticleDefinition * GetDefinition() const
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113

Here is the call graph for this function:

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 }
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:

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)
G4double GetPDGMass() const
double G4double
Definition: G4Types.hh:76
double mag() const

Here is the call graph for this function:

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 }
const char * p
Definition: xmltok.h:285
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:

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 }
const char * p
Definition: xmltok.h:285
double A(double temperature)
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:

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 GetLowEdgeEnergy(size_t binNumber) const
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
int G4int
Definition: G4Types.hh:78
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
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:

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 
743  fNuclearRadius = CalculateNuclearRad(A);
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 G4ParticleHPJENDLHEData::G4double result
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)
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double A(double temperature)
static constexpr double pi
Definition: G4SIunits.hh:75
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:

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 }
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
static constexpr double twopi
Definition: G4SIunits.hh:76
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & boost(double, double, double)
const G4LorentzVector & Get4Momentum() const
double theta() const
G4double GetPDGMass() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
double mag() const
G4double GetTotalMomentum() const

Here is the call graph for this function:

void G4DiffuseElastic::SetHEModelLowLimit ( G4double  value)
inline

Definition at line 257 of file G4DiffuseElastic.hh.

258 {
259  lowEnergyLimitHE = value;
260 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
void G4DiffuseElastic::SetLowestEnergyLimit ( G4double  value)
inline

Definition at line 267 of file G4DiffuseElastic.hh.

268 {
269  lowestEnergyLimit = value;
270 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
void G4DiffuseElastic::SetPlabLowLimit ( G4double  value)
inline

Definition at line 252 of file G4DiffuseElastic.hh.

253 {
254  plabLowLimit = value;
255 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
void G4DiffuseElastic::SetQModelLowLimit ( G4double  value)
inline

Definition at line 262 of file G4DiffuseElastic.hh.

263 {
264  lowEnergyLimitQ = value;
265 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
void G4DiffuseElastic::SetRecoilKinEnergyLimit ( G4double  value)
inline

Definition at line 247 of file G4DiffuseElastic.hh.

248 {
249  lowEnergyRecoilLimit = value;
250 }
const XML_Char int const XML_Char * value
Definition: expat.h:331
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
1347  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1348 
1349 
1350 
1351  G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1352  <<Z<<"; and A = "<<A<<G4endl;
1353 
1354  fElementNumberVector.push_back(fAtomicNumber);
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 
1368  fAngleTable = new G4PhysicsTable(fEnergyBin);
1369 
1370  fWaveVector = partMom/hbarc;
1371 
1372  G4double kR = fWaveVector*fNuclearRadius;
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);
1387  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
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)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void PutValue(size_t index, G4double energy, 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)
static constexpr double degree
Definition: G4SIunits.hh:144
G4double GetPDGMass() const
tuple z
Definition: test.py:28
void insertAt(size_t, G4PhysicsVector *)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
double epsilon(double density, double temperature)
G4double GetIntegrandFunction(G4double theta)

Here is the call graph for this function:

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 }
G4ParticleDefinition * GetDefinition() const
static constexpr double twopi
Definition: G4SIunits.hh:76
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
HepLorentzVector & boost(double, double, double)
G4LorentzVector Get4Momentum() const
double theta() const
G4double GetPDGMass() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
double mag() const

Here is the call graph for this function:

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 }
G4ParticleDefinition * GetDefinition() const
static constexpr double twopi
Definition: G4SIunits.hh:76
G4double GetTotalMomentum() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
HepLorentzVector & boost(double, double, double)
G4LorentzVector Get4Momentum() const
double theta() const
G4double GetPDGMass() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:


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