Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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
 
G4PhysicsLogVector * GetPAIxscVector () const
 
G4PhysicsLogVector * GetPAIdEdxVector () const
 
G4PhysicsLogVector * GetPAIphotonVector () const
 
G4PhysicsLogVector * GetPAIelectronVector () const
 
G4PhysicsLogVector * GetChCosSqVector () const
 
G4PhysicsLogVector * GetChWidthVector () const
 

Detailed Description

Definition at line 47 of file G4InitXscPAI.hh.

Constructor & Destructor Documentation

◆ G4InitXscPAI()

G4InitXscPAI::G4InitXscPAI ( const G4MaterialCutsCouple *  matCC)
explicit

Definition at line 69 of file G4InitXscPAI.cc.

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

◆ ~G4InitXscPAI()

G4InitXscPAI::~G4InitXscPAI ( )
virtual

Definition at line 114 of file G4InitXscPAI.cc.

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

Member Function Documentation

◆ DifPAIdEdx()

G4double G4InitXscPAI::DifPAIdEdx ( G4double  omega)

Definition at line 476 of file G4InitXscPAI.cc.

477{
478 G4double dEdx = omega*DifPAIxSection(omega);
479 return dEdx;
480}
G4double
double G4double
Definition: G4Types.hh:83
G4InitXscPAI::DifPAIxSection
G4double DifPAIxSection(G4double omega)
Definition: G4InitXscPAI.cc:413

Referenced by IntegralPAIdEdx().

◆ DifPAIxSection()

G4double G4InitXscPAI::DifPAIxSection ( G4double  omega)

Definition at line 413 of file G4InitXscPAI.cc.

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

Referenced by DifPAIdEdx(), and IntegralPAIxSection().

◆ GetBinPAI()

G4int G4InitXscPAI::GetBinPAI ( ) const
inline

Definition at line 105 of file G4InitXscPAI.hh.

105{ return fPAIbin ; }

◆ GetChCosSqVector()

G4PhysicsLogVector * G4InitXscPAI::GetChCosSqVector ( ) const
inline

Definition at line 116 of file G4InitXscPAI.hh.

116{ return fChCosSqVector;}

◆ GetChWidthVector()

G4PhysicsLogVector * G4InitXscPAI::GetChWidthVector ( ) const
inline

Definition at line 117 of file G4InitXscPAI.hh.

117{ return fChWidthVector;}

◆ GetIntervalNumber()

G4int G4InitXscPAI::GetIntervalNumber ( ) const
inline

Definition at line 104 of file G4InitXscPAI.hh.

104{ return fIntervalNumber ; }

◆ GetMatSandiaMatrix()

G4double G4InitXscPAI::GetMatSandiaMatrix ( G4int  i,
G4int  j 
) const
inline

Definition at line 109 of file G4InitXscPAI.hh.

110 { return (*(*fMatSandiaMatrix)[i])[j]; }

◆ GetNormalizationCof()

G4double G4InitXscPAI::GetNormalizationCof ( ) const
inline

Definition at line 107 of file G4InitXscPAI.hh.

107{ return fNormalizationCof ; }

◆ GetPAIdEdxVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIdEdxVector ( ) const
inline

Definition at line 113 of file G4InitXscPAI.hh.

113{ return fPAIdEdxVector;}

◆ GetPAIelectronVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIelectronVector ( ) const
inline

Definition at line 115 of file G4InitXscPAI.hh.

115{ return fPAIelectronVector;}

◆ GetPAIphotonVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIphotonVector ( ) const
inline

Definition at line 114 of file G4InitXscPAI.hh.

114{ return fPAIphotonVector;}

◆ GetPAIxscVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIxscVector ( ) const
inline

Definition at line 112 of file G4InitXscPAI.hh.

112{ return fPAIxscVector;}

◆ GetPhotonLambda()

G4double G4InitXscPAI::GetPhotonLambda ( G4double  omega)

Definition at line 932 of file G4InitXscPAI.cc.

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

◆ GetStepCerenkovLoss()

G4double G4InitXscPAI::GetStepCerenkovLoss ( G4double  step)

Definition at line 981 of file G4InitXscPAI.cc.

982{
983 G4double loss = 0.0 ;
984 loss *= step;
985
986 return loss ;
987}

◆ GetStepEnergyLoss()

G4double G4InitXscPAI::GetStepEnergyLoss ( G4double  step)

Definition at line 969 of file G4InitXscPAI.cc.

970{
971 G4double loss = 0.0 ;
972 loss *= step;
973
974 return loss ;
975}

◆ GetStepPlasmonLoss()

G4double G4InitXscPAI::GetStepPlasmonLoss ( G4double  step)

Definition at line 993 of file G4InitXscPAI.cc.

994{
995
996
997 G4double loss = 0.0 ;
998 loss *= step;
999 return loss ;
1000}

◆ ImPartDielectricConst()

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

Definition at line 293 of file G4InitXscPAI.cc.

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

Referenced by DifPAIxSection(), IntegralCherenkov(), ModuleSqDielectricConst(), PAIdNdxCherenkov(), and PAIdNdxPlasmon().

◆ IntegralCherenkov()

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

Definition at line 757 of file G4InitXscPAI.cc.

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

◆ IntegralPAIdEdx()

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

Definition at line 677 of file G4InitXscPAI.cc.

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

◆ IntegralPAIxSection()

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

Definition at line 596 of file G4InitXscPAI.cc.

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

◆ IntegralPlasmon()

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

Definition at line 853 of file G4InitXscPAI.cc.

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

◆ IntegralTerm()

G4double G4InitXscPAI::IntegralTerm ( G4double  omega)

Definition at line 255 of file G4InitXscPAI.cc.

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

Referenced by DifPAIxSection(), and PAIdNdxPlasmon().

◆ KillCloseIntervals()

void G4InitXscPAI::KillCloseIntervals ( )

Definition at line 130 of file G4InitXscPAI.cc.

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

Referenced by G4InitXscPAI().

◆ ModuleSqDielectricConst()

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

Definition at line 319 of file G4InitXscPAI.cc.

321{
322 G4double eIm2, eRe2, result;
323
324 result = ImPartDielectricConst(k,omega);
325 eIm2 = result*result;
326
327 result = RePartDielectricConst(omega);
328 eRe2 = result*result;
329
330 result = eIm2 + eRe2;
331
332 return result ;
333}

Referenced by IntegralCherenkov().

◆ Normalisation()

void G4InitXscPAI::Normalisation ( )

Definition at line 161 of file G4InitXscPAI.cc.

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

Referenced by G4InitXscPAI().

◆ PAIdNdxCherenkov()

G4double G4InitXscPAI::PAIdNdxCherenkov ( G4double  omega)

Definition at line 486 of file G4InitXscPAI.cc.

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

Referenced by IntegralCherenkov().

◆ PAIdNdxPlasmon()

G4double G4InitXscPAI::PAIdNdxPlasmon ( G4double  omega)

Definition at line 550 of file G4InitXscPAI.cc.

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

Referenced by IntegralPlasmon().

◆ RePartDielectricConst()

G4double G4InitXscPAI::RePartDielectricConst ( G4double  energy)

Definition at line 342 of file G4InitXscPAI.cc.

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

Referenced by DifPAIxSection(), IntegralCherenkov(), ModuleSqDielectricConst(), PAIdNdxCherenkov(), and PAIdNdxPlasmon().

◆ RutherfordIntegral()

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

Definition at line 231 of file G4InitXscPAI.cc.

234{
235 G4double c1, c2, c3, a1, a2, a3, a4 ;
236
237 a1 = (*(*fMatSandiaMatrix)[k])[1];
238 a2 = (*(*fMatSandiaMatrix)[k])[2];
239 a3 = (*(*fMatSandiaMatrix)[k])[3];
240 a4 = (*(*fMatSandiaMatrix)[k])[4];
241 // G4cout<<"RI: x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
242 c1 = (x2 - x1)/x1/x2 ;
243 c2 = (x2 - x1)*(x2 + x1)/x1/x1/x2/x2 ;
244 c3 = (x2 - x1)*(x1*x1 + x1*x2 + x2*x2)/x1/x1/x1/x2/x2/x2 ;
245 // G4cout<<" RI: c1 = "<<c1<<"; "<<"c2 = "<<c2<<"; "<<"c3 = "<<c3<<G4endl;
246
247 return a1*log(x2/x1) + a2*c1 + a3*c2/2 + a4*c3/3 ;
248
249} // end of RutherfordIntegral

Referenced by IntegralTerm(), and Normalisation().

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