Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4InitXscPAI Class Reference

#include <G4InitXscPAI.hh>

Public Member Functions

 G4InitXscPAI (const G4MaterialCutsCouple *matCC)
 
virtual ~G4InitXscPAI ()
 
void KillCloseIntervals ()
 
void Normalisation ()
 
G4double RutherfordIntegral (G4int intervalNumber, G4double limitLow, G4double limitHigh)
 
G4double IntegralTerm (G4double omega)
 
G4double ImPartDielectricConst (G4int intervalNumber, G4double energy)
 
G4double RePartDielectricConst (G4double energy)
 
G4double ModuleSqDielectricConst (G4int intervalNumber, G4double energy)
 
G4double DifPAIxSection (G4double omega)
 
G4double DifPAIdEdx (G4double omega)
 
G4double PAIdNdxCherenkov (G4double omega)
 
G4double PAIdNdxPlasmon (G4double omega)
 
void IntegralPAIxSection (G4double bg2, G4double Tmax)
 
void IntegralCherenkov (G4double bg2, G4double Tmax)
 
void IntegralPlasmon (G4double bg2, G4double Tmax)
 
void IntegralPAIdEdx (G4double bg2, G4double Tmax)
 
G4double GetPhotonLambda (G4double omega)
 
G4double GetStepEnergyLoss (G4double step)
 
G4double GetStepCerenkovLoss (G4double step)
 
G4double GetStepPlasmonLoss (G4double step)
 
G4int GetIntervalNumber () const
 
G4int GetBinPAI () const
 
G4double GetNormalizationCof () const
 
G4double GetMatSandiaMatrix (G4int i, G4int j) const
 
G4PhysicsLogVectorGetPAIxscVector () const
 
G4PhysicsLogVectorGetPAIdEdxVector () const
 
G4PhysicsLogVectorGetPAIphotonVector () const
 
G4PhysicsLogVectorGetPAIelectronVector () const
 
G4PhysicsLogVectorGetChCosSqVector () const
 
G4PhysicsLogVectorGetChWidthVector () const
 

Detailed Description

Definition at line 48 of file G4InitXscPAI.hh.

Constructor & Destructor Documentation

G4InitXscPAI::G4InitXscPAI ( const G4MaterialCutsCouple matCC)
explicit

Definition at line 70 of file G4InitXscPAI.cc.

71  : fPAIxscVector(nullptr),
72  fPAIdEdxVector(nullptr),
73  fPAIphotonVector(nullptr),
74  fPAIelectronVector(nullptr),
75  fChCosSqVector(nullptr),
76  fChWidthVector(nullptr)
77 {
78  G4int i, j, matIndex;
79 
80  fDensity = matCC->GetMaterial()->GetDensity();
81  fElectronDensity = matCC->GetMaterial()->GetElectronDensity();
82  matIndex = matCC->GetMaterial()->GetIndex();
83 
84  fSandia = new G4SandiaTable(matIndex);
85  fIntervalNumber = fSandia->GetMaxInterval()-1;
86 
87  fMatSandiaMatrix = new G4OrderedTable();
88 
89  for (i = 0; i < fIntervalNumber; i++)
90  {
91  fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
92  }
93  for (i = 0; i < fIntervalNumber; i++)
94  {
95  (*(*fMatSandiaMatrix)[i])[0] = fSandia->GetSandiaMatTable(i,0);
96 
97  for(j = 1; j < 5 ; j++)
98  {
99  (*(*fMatSandiaMatrix)[i])[j] = fSandia->GetSandiaMatTable(i,j)*fDensity;
100  }
101  }
103  Normalisation();
104  fBetaGammaSq = fTmax = 0.0;
105  fIntervalTmax = fCurrentInterval = 0;
106 }
size_t GetIndex() const
Definition: G4Material.hh:262
G4double GetDensity() const
Definition: G4Material.hh:180
G4int GetMaxInterval() const
void KillCloseIntervals()
int G4int
Definition: G4Types.hh:78
void Normalisation()
G4double GetElectronDensity() const
Definition: G4Material.hh:217
G4double GetSandiaMatTable(G4int, G4int) const
const G4Material * GetMaterial() const

Here is the call graph for this function:

G4InitXscPAI::~G4InitXscPAI ( )
virtual

Definition at line 115 of file G4InitXscPAI.cc.

116 {
117  if(fPAIxscVector) delete fPAIxscVector;
118  if(fPAIdEdxVector) delete fPAIdEdxVector;
119  if(fPAIphotonVector) delete fPAIphotonVector;
120  if(fPAIelectronVector) delete fPAIelectronVector;
121  if(fChCosSqVector) delete fChCosSqVector;
122  if(fChWidthVector) delete fChWidthVector;
123  delete fSandia;
124  delete fMatSandiaMatrix;
125 }

Member Function Documentation

G4double G4InitXscPAI::DifPAIdEdx ( G4double  omega)

Definition at line 477 of file G4InitXscPAI.cc.

478 {
479  G4double dEdx = omega*DifPAIxSection(omega);
480  return dEdx;
481 }
G4double DifPAIxSection(G4double omega)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4InitXscPAI::DifPAIxSection ( G4double  omega)

Definition at line 414 of file G4InitXscPAI.cc.

415 {
416  G4int i = fCurrentInterval;
417  G4double betaGammaSq = fBetaGammaSq;
418  G4double integralTerm = IntegralTerm(omega);
419  G4double be2,cof,x1,x2,x3,x4,x5,x6,x7,x8,result ;
420  G4double epsilonRe = RePartDielectricConst(omega);
421  G4double epsilonIm = ImPartDielectricConst(i,omega);
422  G4double be4 ;
423  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
424  static const G4double betaBohr4 = betaBohr2*betaBohr2*4.0 ;
425  be2 = betaGammaSq/(1 + betaGammaSq) ;
426  be4 = be2*be2 ;
427 
428  cof = 1 ;
429  x1 = log(2*electron_mass_c2/omega) ;
430 
431  if( betaGammaSq < 0.01 ) x2 = log(be2) ;
432  else
433  {
434  x2 = -log( (1/betaGammaSq - epsilonRe)*
435  (1/betaGammaSq - epsilonRe) +
436  epsilonIm*epsilonIm )/2 ;
437  }
438  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
439  {
440  x6=0 ;
441  }
442  else
443  {
444  x3 = -epsilonRe + 1/betaGammaSq ;
445  x5 = -1 - epsilonRe + be2*((1 +epsilonRe)*(1 + epsilonRe) +
446  epsilonIm*epsilonIm) ;
447 
448  x7 = atan2(epsilonIm,x3) ;
449  x6 = x5 * x7 ;
450  }
451  // if(fImPartDielectricConst[i] == 0) x6 = 0 ;
452 
453  x4 = ((x1 + x2)*epsilonIm + x6)/hbarc ;
454  // if( x4 < 0.0 ) x4 = 0.0 ;
455  x8 = (1 + epsilonRe)*(1 + epsilonRe) +
456  epsilonIm*epsilonIm;
457 
458  result = (x4 + cof*integralTerm/omega/omega) ;
459  if(result < 1.0e-8) result = 1.0e-8 ;
460  result *= fine_structure_const/be2/pi ;
461  // result *= (1-exp(-beta/betaBohr))*(1-exp(-beta/betaBohr)) ;
462  // result *= (1-exp(-be2/betaBohr2)) ;
463  result *= (1-exp(-be4/betaBohr4)) ;
464  if(fDensity >= fSolidDensity)
465  {
466  result /= x8 ;
467  }
468  return result ;
469 
470 } // end of DifPAIxSection
G4double G4ParticleHPJENDLHEData::G4double result
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double IntegralTerm(G4double omega)

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4InitXscPAI::GetBinPAI ( ) const
inline

Definition at line 106 of file G4InitXscPAI.hh.

106 { return fPAIbin ; }
G4PhysicsLogVector* G4InitXscPAI::GetChCosSqVector ( ) const
inline

Definition at line 117 of file G4InitXscPAI.hh.

117 { return fChCosSqVector;}
G4PhysicsLogVector* G4InitXscPAI::GetChWidthVector ( ) const
inline

Definition at line 118 of file G4InitXscPAI.hh.

118 { return fChWidthVector;}
G4int G4InitXscPAI::GetIntervalNumber ( ) const
inline

Definition at line 105 of file G4InitXscPAI.hh.

105 { return fIntervalNumber ; }
G4double G4InitXscPAI::GetMatSandiaMatrix ( G4int  i,
G4int  j 
) const
inline

Definition at line 110 of file G4InitXscPAI.hh.

111  { return (*(*fMatSandiaMatrix)[i])[j]; }
G4double G4InitXscPAI::GetNormalizationCof ( ) const
inline

Definition at line 108 of file G4InitXscPAI.hh.

108 { return fNormalizationCof ; }
G4PhysicsLogVector* G4InitXscPAI::GetPAIdEdxVector ( ) const
inline

Definition at line 114 of file G4InitXscPAI.hh.

114 { return fPAIdEdxVector;}
G4PhysicsLogVector* G4InitXscPAI::GetPAIelectronVector ( ) const
inline

Definition at line 116 of file G4InitXscPAI.hh.

116 { return fPAIelectronVector;}
G4PhysicsLogVector* G4InitXscPAI::GetPAIphotonVector ( ) const
inline

Definition at line 115 of file G4InitXscPAI.hh.

115 { return fPAIphotonVector;}
G4PhysicsLogVector* G4InitXscPAI::GetPAIxscVector ( ) const
inline

Definition at line 113 of file G4InitXscPAI.hh.

113 { return fPAIxscVector;}
G4double G4InitXscPAI::GetPhotonLambda ( G4double  omega)

Definition at line 933 of file G4InitXscPAI.cc.

934 {
935  G4int i ;
936  G4double omega2, omega3, omega4, a1, a2, a3, a4, lambda ;
937 
938  omega2 = omega*omega ;
939  omega3 = omega2*omega ;
940  omega4 = omega2*omega2 ;
941 
942  for(i = 0; i < fIntervalNumber;i++)
943  {
944  if( omega < (*(*fMatSandiaMatrix)[i])[0] ) break ;
945  }
946  if( i == 0 )
947  {
948  G4cout<<"Warning: energy in G4InitXscPAI::GetPhotonLambda < I1"<<G4endl;
949  }
950  else i-- ;
951 
952  a1 = (*(*fMatSandiaMatrix)[i])[1];
953  a2 = (*(*fMatSandiaMatrix)[i])[2];
954  a3 = (*(*fMatSandiaMatrix)[i])[3];
955  a4 = (*(*fMatSandiaMatrix)[i])[4];
956 
957  lambda = 1./(a1/omega + a2/omega2 + a3/omega3 + a4/omega4);
958 
959  return lambda ;
960 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double G4InitXscPAI::GetStepCerenkovLoss ( G4double  step)

Definition at line 982 of file G4InitXscPAI.cc.

983 {
984  G4double loss = 0.0 ;
985  loss *= step;
986 
987  return loss ;
988 }
double G4double
Definition: G4Types.hh:76
G4double G4InitXscPAI::GetStepEnergyLoss ( G4double  step)

Definition at line 970 of file G4InitXscPAI.cc.

971 {
972  G4double loss = 0.0 ;
973  loss *= step;
974 
975  return loss ;
976 }
double G4double
Definition: G4Types.hh:76
G4double G4InitXscPAI::GetStepPlasmonLoss ( G4double  step)

Definition at line 994 of file G4InitXscPAI.cc.

995 {
996 
997 
998  G4double loss = 0.0 ;
999  loss *= step;
1000  return loss ;
1001 }
double G4double
Definition: G4Types.hh:76
G4double G4InitXscPAI::ImPartDielectricConst ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 294 of file G4InitXscPAI.cc.

296 {
297  G4double energy2,energy3,energy4,a1,a2,a3,a4,result;
298 
299  a1 = (*(*fMatSandiaMatrix)[k])[1];
300  a2 = (*(*fMatSandiaMatrix)[k])[2];
301  a3 = (*(*fMatSandiaMatrix)[k])[3];
302  a4 = (*(*fMatSandiaMatrix)[k])[4];
303 
304  energy2 = energy1*energy1;
305  energy3 = energy2*energy1;
306  energy4 = energy3*energy1;
307 
308  result = a1/energy1+a2/energy2+a3/energy3+a4/energy4 ;
309  result *= hbarc/energy1 ;
310 
311  return result ;
312 
313 } // end of ImPartDielectricConst
G4double G4ParticleHPJENDLHEData::G4double result
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

void G4InitXscPAI::IntegralCherenkov ( G4double  bg2,
G4double  Tmax 
)

Definition at line 758 of file G4InitXscPAI.cc.

759 {
760  G4int i,k,i1,i2;
761  G4double energy1, energy2, beta2, module2, cos2, width, result = 0.;
762 
763  fBetaGammaSq = bg2;
764  fTmax = Tmax;
765  beta2 = bg2/(1+bg2);
766 
767  if(fPAIphotonVector) delete fPAIphotonVector;
768  if(fChCosSqVector) delete fChCosSqVector;
769  if(fChWidthVector) delete fChWidthVector;
770 
771  fPAIphotonVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
772  fChCosSqVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
773  fChWidthVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
774 
775  fPAIphotonVector->PutValue(fPAIbin-1,result);
776  fChCosSqVector->PutValue(fPAIbin-1,1.);
777  fChWidthVector->PutValue(fPAIbin-1,1e-7);
778 
779  for( i = fIntervalNumber - 1; i >= 0; i-- )
780  {
781  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
782  }
783  if (i < 0) i = 0; // Tmax should be more than
784  // first ionisation potential
785  fIntervalTmax = i;
786 
788 
789  for( k = fPAIbin - 2; k >= 0; k-- )
790  {
791  energy1 = fPAIphotonVector->GetLowEdgeEnergy(k);
792  energy2 = fPAIphotonVector->GetLowEdgeEnergy(k+1);
793 
794  for( i = fIntervalTmax; i >= 0; i-- )
795  {
796  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
797  }
798  if(i < 0) i = 0;
799  i2 = i;
800 
801  for( i = fIntervalTmax; i >= 0; i-- )
802  {
803  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
804  }
805  if(i < 0) i = 0;
806  i1 = i;
807 
808  module2 = ModuleSqDielectricConst(i1,energy1);
809  cos2 = RePartDielectricConst(energy1)/module2/beta2;
810  width = ImPartDielectricConst(i1,energy1)/module2/beta2;
811 
812  fChCosSqVector->PutValue(k,cos2);
813  fChWidthVector->PutValue(k,width);
814 
815  if( i1 == i2 )
816  {
817  fCurrentInterval = i1;
818  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxCherenkov,
819  energy1,energy2);
820  fPAIphotonVector->PutValue(k,result);
821 
822  }
823  else
824  {
825  for( i = i2; i >= i1; i-- )
826  {
827  fCurrentInterval = i;
828 
829  if( i==i2 ) result += integral.Legendre10(this,
831  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
832 
833  else if( i == i1 ) result += integral.Legendre10(this,
835  (*(*fMatSandiaMatrix)[i+1])[0]);
836 
837  else result += integral.Legendre10(this,
839  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
840  }
841  fPAIphotonVector->PutValue(k,result);
842  }
843  // G4cout<<k<<"\t"<<result<<G4endl;
844  }
845  return;
846 } // end of IntegralCerenkov
G4double G4ParticleHPJENDLHEData::G4double result
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
G4double PAIdNdxCherenkov(G4double omega)
G4double ModuleSqDielectricConst(G4int intervalNumber, G4double energy)
#define width
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double RePartDielectricConst(G4double energy)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
void PutValue(size_t index, G4double theValue)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4InitXscPAI::IntegralPAIdEdx ( G4double  bg2,
G4double  Tmax 
)

Definition at line 678 of file G4InitXscPAI.cc.

679 {
680  G4int i,k,i1,i2;
681  G4double energy1, energy2, result = 0.;
682 
683  fBetaGammaSq = bg2;
684  fTmax = Tmax;
685 
686  if(fPAIdEdxVector) delete fPAIdEdxVector;
687 
688  fPAIdEdxVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
689  fPAIdEdxVector->PutValue(fPAIbin-1,result);
690 
691  for( i = fIntervalNumber - 1; i >= 0; i-- )
692  {
693  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
694  }
695  if (i < 0) i = 0; // Tmax should be more than
696  // first ionisation potential
697  fIntervalTmax = i;
698 
700 
701  for( k = fPAIbin - 2; k >= 0; k-- )
702  {
703  energy1 = fPAIdEdxVector->GetLowEdgeEnergy(k);
704  energy2 = fPAIdEdxVector->GetLowEdgeEnergy(k+1);
705 
706  for( i = fIntervalTmax; i >= 0; i-- )
707  {
708  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
709  }
710  if(i < 0) i = 0;
711  i2 = i;
712 
713  for( i = fIntervalTmax; i >= 0; i-- )
714  {
715  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
716  }
717  if(i < 0) i = 0;
718  i1 = i;
719 
720  if( i1 == i2 )
721  {
722  fCurrentInterval = i1;
723  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIdEdx,
724  energy1,energy2);
725  fPAIdEdxVector->PutValue(k,result);
726  }
727  else
728  {
729  for( i = i2; i >= i1; i-- )
730  {
731  fCurrentInterval = i;
732 
733  if( i==i2 ) result += integral.Legendre10(this,
735  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
736 
737  else if( i == i1 ) result += integral.Legendre10(this,
738  &G4InitXscPAI::DifPAIdEdx,energy1,
739  (*(*fMatSandiaMatrix)[i+1])[0]);
740 
741  else result += integral.Legendre10(this,
743  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
744  }
745  fPAIdEdxVector->PutValue(k,result);
746  }
747  // G4cout<<k<<"\t"<<result<<G4endl;
748  }
749  return ;
750 }
G4double G4ParticleHPJENDLHEData::G4double result
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
void PutValue(size_t index, G4double theValue)
G4double DifPAIdEdx(G4double omega)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4InitXscPAI::IntegralPAIxSection ( G4double  bg2,
G4double  Tmax 
)

Definition at line 597 of file G4InitXscPAI.cc.

598 {
599  G4int i,k,i1,i2;
600  G4double energy1, energy2, result = 0.;
601 
602  fBetaGammaSq = bg2;
603  fTmax = Tmax;
604 
605  if(fPAIxscVector) delete fPAIxscVector;
606 
607  fPAIxscVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
608  fPAIxscVector->PutValue(fPAIbin-1,result);
609 
610  for( i = fIntervalNumber - 1; i >= 0; i-- )
611  {
612  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
613  }
614  if (i < 0) i = 0; // Tmax should be more than
615  // first ionisation potential
616  fIntervalTmax = i;
617 
619 
620  for( k = fPAIbin - 2; k >= 0; k-- )
621  {
622  energy1 = fPAIxscVector->GetLowEdgeEnergy(k);
623  energy2 = fPAIxscVector->GetLowEdgeEnergy(k+1);
624 
625  for( i = fIntervalTmax; i >= 0; i-- )
626  {
627  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
628  }
629  if(i < 0) i = 0;
630  i2 = i;
631 
632  for( i = fIntervalTmax; i >= 0; i-- )
633  {
634  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
635  }
636  if(i < 0) i = 0;
637  i1 = i;
638 
639  if( i1 == i2 )
640  {
641  fCurrentInterval = i1;
642  result += integral.Legendre10(this,&G4InitXscPAI::DifPAIxSection,
643  energy1,energy2);
644  fPAIxscVector->PutValue(k,result);
645  }
646  else
647  {
648  for( i = i2; i >= i1; i-- )
649  {
650  fCurrentInterval = i;
651 
652  if( i==i2 ) result += integral.Legendre10(this,
654  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
655 
656  else if( i == i1 ) result += integral.Legendre10(this,
658  (*(*fMatSandiaMatrix)[i+1])[0]);
659 
660  else result += integral.Legendre10(this,
662  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
663  }
664  fPAIxscVector->PutValue(k,result);
665  }
666  // G4cout<<k<<"\t"<<result<<G4endl;
667  }
668  return ;
669 }
G4double G4ParticleHPJENDLHEData::G4double result
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
G4double DifPAIxSection(G4double omega)
void PutValue(size_t index, G4double theValue)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4InitXscPAI::IntegralPlasmon ( G4double  bg2,
G4double  Tmax 
)

Definition at line 854 of file G4InitXscPAI.cc.

855 {
856  G4int i,k,i1,i2;
857  G4double energy1, energy2, result = 0.;
858 
859  fBetaGammaSq = bg2;
860  fTmax = Tmax;
861 
862  if(fPAIelectronVector) delete fPAIelectronVector;
863 
864  fPAIelectronVector = new G4PhysicsLogVector( (*(*fMatSandiaMatrix)[0])[0], fTmax, fPAIbin);
865  fPAIelectronVector->PutValue(fPAIbin-1,result);
866 
867  for( i = fIntervalNumber - 1; i >= 0; i-- )
868  {
869  if( Tmax >= (*(*fMatSandiaMatrix)[i])[0] ) break;
870  }
871  if (i < 0) i = 0; // Tmax should be more than
872  // first ionisation potential
873  fIntervalTmax = i;
874 
876 
877  for( k = fPAIbin - 2; k >= 0; k-- )
878  {
879  energy1 = fPAIelectronVector->GetLowEdgeEnergy(k);
880  energy2 = fPAIelectronVector->GetLowEdgeEnergy(k+1);
881 
882  for( i = fIntervalTmax; i >= 0; i-- )
883  {
884  if( energy2 > (*(*fMatSandiaMatrix)[i])[0] ) break;
885  }
886  if(i < 0) i = 0;
887  i2 = i;
888 
889  for( i = fIntervalTmax; i >= 0; i-- )
890  {
891  if( energy1 > (*(*fMatSandiaMatrix)[i])[0] ) break;
892  }
893  if(i < 0) i = 0;
894  i1 = i;
895 
896  if( i1 == i2 )
897  {
898  fCurrentInterval = i1;
899  result += integral.Legendre10(this,&G4InitXscPAI::PAIdNdxPlasmon,
900  energy1,energy2);
901  fPAIelectronVector->PutValue(k,result);
902  }
903  else
904  {
905  for( i = i2; i >= i1; i-- )
906  {
907  fCurrentInterval = i;
908 
909  if( i==i2 ) result += integral.Legendre10(this,
911  (*(*fMatSandiaMatrix)[i])[0] ,energy2);
912 
913  else if( i == i1 ) result += integral.Legendre10(this,
915  (*(*fMatSandiaMatrix)[i+1])[0]);
916 
917  else result += integral.Legendre10(this,
919  (*(*fMatSandiaMatrix)[i])[0] ,(*(*fMatSandiaMatrix)[i+1])[0]);
920  }
921  fPAIelectronVector->PutValue(k,result);
922  }
923  // G4cout<<k<<"\t"<<result<<G4endl;
924  }
925  return;
926 } // end of IntegralPlasmon
G4double G4ParticleHPJENDLHEData::G4double result
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
void PutValue(size_t index, G4double theValue)
double G4double
Definition: G4Types.hh:76
G4double PAIdNdxPlasmon(G4double omega)

Here is the call graph for this function:

G4double G4InitXscPAI::IntegralTerm ( G4double  omega)

Definition at line 256 of file G4InitXscPAI.cc.

257 {
258  G4int i;
259  G4double energy1, energy2, result = 0.;
260 
261  for( i = 0; i <= fIntervalTmax; i++ )
262  {
263  if(i == fIntervalTmax)
264  {
265  energy1 = (*(*fMatSandiaMatrix)[i])[0];
266  result += RutherfordIntegral(i,energy1,omega);
267  }
268  else
269  {
270  if( omega <= (*(*fMatSandiaMatrix)[i+1])[0])
271  {
272  energy1 = (*(*fMatSandiaMatrix)[i])[0];
273  result += RutherfordIntegral(i,energy1,omega);
274  break;
275  }
276  else
277  {
278  energy1 = (*(*fMatSandiaMatrix)[i])[0];
279  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
280  result += RutherfordIntegral(i,energy1,energy2);
281  }
282  }
283  // G4cout<<"IntegralTerm<<"("<<omega<<")"<<" = "<<result<<G4endl;
284  }
285  return result;
286 }
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4InitXscPAI::KillCloseIntervals ( )

Definition at line 131 of file G4InitXscPAI.cc.

132 {
133  G4int i, j, k;
134  G4double energy1, energy2;
135 
136  for( i = 0 ; i < fIntervalNumber - 1 ; i++ )
137  {
138  energy1 = (*(*fMatSandiaMatrix)[i])[0];
139  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
140 
141  if( energy2 - energy1 > 1.5*fDelta*(energy1 + energy2) ) continue ;
142  else
143  {
144  for(j = i; j < fIntervalNumber-1; j++)
145  {
146  for( k = 0; k < 5; k++ )
147  {
148  (*(*fMatSandiaMatrix)[j])[k] = (*(*fMatSandiaMatrix)[j+1])[k];
149  }
150  }
151  fIntervalNumber-- ;
152  i-- ;
153  }
154  }
155 
156 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76

Here is the caller graph for this function:

G4double G4InitXscPAI::ModuleSqDielectricConst ( G4int  intervalNumber,
G4double  energy 
)

Definition at line 320 of file G4InitXscPAI.cc.

322 {
323  G4double eIm2, eRe2, result;
324 
325  result = ImPartDielectricConst(k,omega);
326  eIm2 = result*result;
327 
328  result = RePartDielectricConst(omega);
329  eRe2 = result*result;
330 
331  result = eIm2 + eRe2;
332 
333  return result ;
334 }
G4double G4ParticleHPJENDLHEData::G4double result
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double RePartDielectricConst(G4double energy)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4InitXscPAI::Normalisation ( )

Definition at line 162 of file G4InitXscPAI.cc.

163 {
164  G4int i, j;
165  G4double energy1, energy2, /*delta,*/ cof; // , shift;
166 
167  energy1 = (*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
168  energy2 = 2.*(*(*fMatSandiaMatrix)[fIntervalNumber-1])[0];
169 
170 
171  cof = RutherfordIntegral(fIntervalNumber-1,energy1,energy2);
172 
173  for( i = fIntervalNumber-2; i >= 0; i-- )
174  {
175  energy1 = (*(*fMatSandiaMatrix)[i])[0];
176  energy2 = (*(*fMatSandiaMatrix)[i+1])[0];
177 
178  cof += RutherfordIntegral(i,energy1,energy2);
179  // G4cout<<"norm. cof = "<<cof<<G4endl;
180  }
181  fNormalizationCof = 2*pi*pi*hbarc*hbarc*fine_structure_const/electron_mass_c2 ;
182  fNormalizationCof *= fElectronDensity;
183  //delta = fNormalizationCof - cof;
184  fNormalizationCof /= cof;
185  // G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof
186  // <<"; at delta ="<<delta<<G4endl ;
187 
188  for (i = 0; i < fIntervalNumber; i++) // renormalisation on QM sum rule
189  {
190  for(j = 1; j < 5 ; j++)
191  {
192  (*(*fMatSandiaMatrix)[i])[j] *= fNormalizationCof;
193  }
194  }
195  /*
196  if(delta > 0) // shift the first energy interval
197  {
198  for(i=1;i<100;i++)
199  {
200  energy1 = (1.-i/100.)*(*(*fMatSandiaMatrix)[0])[0];
201  energy2 = (*(*fMatSandiaMatrix)[0])[0];
202  shift = RutherfordIntegral(0,energy1,energy2);
203  G4cout<<shift<<"\t";
204  if(shift >= delta) break;
205  }
206  (*(*fMatSandiaMatrix)[0])[0] = energy1;
207  cof += shift;
208  }
209  else if(delta < 0)
210  {
211  for(i=1;i<100;i++)
212  {
213  energy1 = (*(*fMatSandiaMatrix)[0])[0];
214  energy2 = (*(*fMatSandiaMatrix)[0])[0] +
215  ( (*(*fMatSandiaMatrix)[0])[0] - (*(*fMatSandiaMatrix)[0])[0] )*i/100.;
216  shift = RutherfordIntegral(0,energy1,energy2);
217  if( shift >= std::abs(delta) ) break;
218  }
219  (*(*fMatSandiaMatrix)[0])[0] = energy2;
220  cof -= shift;
221  }
222  G4cout<<G4cout<<"G4InitXscPAI::fNormalizationCof/cof = "<<fNormalizationCof/cof
223  <<"; at delta ="<<delta<<" and i = "<<i<<G4endl ;
224  */
225 }
int G4int
Definition: G4Types.hh:78
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4InitXscPAI::PAIdNdxCherenkov ( G4double  omega)

Definition at line 487 of file G4InitXscPAI.cc.

488 {
489  G4int i = fCurrentInterval;
490  G4double betaGammaSq = fBetaGammaSq;
491  G4double epsilonRe = RePartDielectricConst(omega);
492  G4double epsilonIm = ImPartDielectricConst(i,omega);
493 
494  G4double /*cof,*/ logarithm, x3, x5, argument, modul2, dNdxC ;
495  G4double be2, be4;
496 
497  //cof = 1.0 ;
498  static const G4double cofBetaBohr = 4.0 ;
499  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
500  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
501 
502  be2 = betaGammaSq/(1 + betaGammaSq) ;
503  be4 = be2*be2 ;
504 
505  if( betaGammaSq < 0.01 ) logarithm = log(1.0+betaGammaSq) ; // 0.0 ;
506  else
507  {
508  logarithm = -log( (1/betaGammaSq - epsilonRe)*
509  (1/betaGammaSq - epsilonRe) +
510  epsilonIm*epsilonIm )*0.5 ;
511  logarithm += log(1+1.0/betaGammaSq) ;
512  }
513 
514  if( epsilonIm == 0.0 || betaGammaSq < 0.01 )
515  {
516  argument = 0.0 ;
517  }
518  else
519  {
520  x3 = -epsilonRe + 1.0/betaGammaSq ;
521  x5 = -1.0 - epsilonRe +
522  be2*((1.0 +epsilonRe)*(1.0 + epsilonRe) +
523  epsilonIm*epsilonIm) ;
524  if( x3 == 0.0 ) argument = 0.5*pi;
525  else argument = atan2(epsilonIm,x3) ;
526  argument *= x5 ;
527  }
528  dNdxC = ( logarithm*epsilonIm + argument )/hbarc ;
529 
530  if(dNdxC < 1.0e-8) dNdxC = 1.0e-8 ;
531 
532  dNdxC *= fine_structure_const/be2/pi ;
533 
534  dNdxC *= (1-exp(-be4/betaBohr4)) ;
535 
536  if(fDensity >= fSolidDensity)
537  {
538  modul2 = (1.0 + epsilonRe)*(1.0 + epsilonRe) +
539  epsilonIm*epsilonIm;
540  dNdxC /= modul2 ;
541  }
542  return dNdxC ;
543 
544 } // end of PAIdNdxCerenkov
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4InitXscPAI::PAIdNdxPlasmon ( G4double  omega)

Definition at line 551 of file G4InitXscPAI.cc.

552 {
553  G4int i = fCurrentInterval;
554  G4double betaGammaSq = fBetaGammaSq;
555  G4double integralTerm = IntegralTerm(omega);
556  G4double epsilonRe = RePartDielectricConst(omega);
557  G4double epsilonIm = ImPartDielectricConst(i,omega);
558 
559  G4double cof, resonance, modul2, dNdxP ;
560  G4double be2, be4;
561 
562  cof = 1 ;
563  static const G4double cofBetaBohr = 4.0 ;
564  static const G4double betaBohr2 = fine_structure_const*fine_structure_const ;
565  static const G4double betaBohr4 = betaBohr2*betaBohr2*cofBetaBohr ;
566 
567  be2 = betaGammaSq/(1 + betaGammaSq) ;
568  be4 = be2*be2 ;
569 
570  resonance = log(2*electron_mass_c2*be2/omega) ;
571  resonance *= epsilonIm/hbarc ;
572 
573 
574  dNdxP = ( resonance + cof*integralTerm/omega/omega ) ;
575 
576  if( dNdxP < 1.0e-8 ) dNdxP = 1.0e-8 ;
577 
578  dNdxP *= fine_structure_const/be2/pi ;
579  dNdxP *= (1-exp(-be4/betaBohr4)) ;
580 
581  if( fDensity >= fSolidDensity )
582  {
583  modul2 = (1 + epsilonRe)*(1 + epsilonRe) +
584  epsilonIm*epsilonIm;
585  dNdxP /= modul2 ;
586  }
587  return dNdxP ;
588 
589 } // end of PAIdNdxPlasmon
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
G4double RePartDielectricConst(G4double energy)
int G4int
Definition: G4Types.hh:78
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
G4double IntegralTerm(G4double omega)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4InitXscPAI::RePartDielectricConst ( G4double  energy)

Definition at line 343 of file G4InitXscPAI.cc.

344 {
345  G4int i;
346  G4double x0, x02, x03, x04, x05, x1, x2, a1,a2,a3,a4,xx1 ,xx2 , xx12,
347  c1, c2, c3, cof1, cof2, xln1, xln2, xln3, result ;
348 
349  x0 = enb ;
350  result = 0 ;
351 
352  for( i = 0; i < fIntervalNumber-1; i++)
353  {
354  x1 = (*(*fMatSandiaMatrix)[i])[0];
355  x2 = (*(*fMatSandiaMatrix)[i+1])[0] ;
356 
357  a1 = (*(*fMatSandiaMatrix)[i])[1];
358  a2 = (*(*fMatSandiaMatrix)[i])[2];
359  a3 = (*(*fMatSandiaMatrix)[i])[3];
360  a4 = (*(*fMatSandiaMatrix)[i])[4];
361 
362  if( std::abs(x0-x1) < 0.5*(x0+x1)*fDelta )
363  {
364  if(x0 >= x1) x0 = x1*(1+fDelta);
365  else x0 = x1*(1-fDelta);
366  }
367  if( std::abs(x0-x2) < 0.5*(x0+x2)*fDelta )
368  {
369  if(x0 >= x2) x0 = x2*(1+fDelta);
370  else x0 = x2*(1-fDelta);
371  }
372  xx1 = x1 - x0 ;
373  xx2 = x2 - x0 ;
374  xx12 = xx2/xx1 ;
375 
376  if( xx12 < 0 ) xx12 = -xx12;
377 
378  xln1 = log(x2/x1) ;
379  xln2 = log(xx12) ;
380  xln3 = log((x2 + x0)/(x1 + x0)) ;
381 
382  x02 = x0*x0 ;
383  x03 = x02*x0 ;
384  x04 = x03*x0 ;
385  x05 = x04*x0;
386 
387  c1 = (x2 - x1)/x1/x2 ;
388  c2 = (x2 - x1)*(x2 +x1)/x1/x1/x2/x2 ;
389  c3 = (x2 -x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
390 
391  result -= (a1/x02 + a3/x04)*xln1 ;
392  result -= (a2/x02 + a4/x04)*c1 ;
393  result -= a3*c2/2/x02 ;
394  result -= a4*c3/3/x02 ;
395 
396  cof1 = a1/x02 + a3/x04 ;
397  cof2 = a2/x03 + a4/x05 ;
398 
399  result += 0.5*(cof1 +cof2)*xln2 ;
400  result += 0.5*(cof1 - cof2)*xln3 ;
401  }
402  result *= 2*hbarc/pi ;
403 
404  return result ;
405 
406 } // end of RePartDielectricConst
G4double G4ParticleHPJENDLHEData::G4double result
static c2_factory< G4double > c2
int G4int
Definition: G4Types.hh:78
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
tuple c1
Definition: plottest35.py:14

Here is the caller graph for this function:

G4double G4InitXscPAI::RutherfordIntegral ( G4int  intervalNumber,
G4double  limitLow,
G4double  limitHigh 
)

Definition at line 232 of file G4InitXscPAI.cc.

235 {
236  G4double c1, c2, c3, a1, a2, a3, a4 ;
237 
238  a1 = (*(*fMatSandiaMatrix)[k])[1];
239  a2 = (*(*fMatSandiaMatrix)[k])[2];
240  a3 = (*(*fMatSandiaMatrix)[k])[3];
241  a4 = (*(*fMatSandiaMatrix)[k])[4];
242  // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
243  c1 = (x2 - x1)/x1/x2 ;
244  c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
245  c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
246  // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
247 
248  return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
249 
250 } // end of RutherfordIntegral
static c2_factory< G4double > c2
double G4double
Definition: G4Types.hh:76
tuple c1
Definition: plottest35.py:14

Here is the caller graph for this function:


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