Geant4 9.6.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 48 of file G4InitXscPAI.hh.

Constructor & Destructor Documentation

◆ G4InitXscPAI()

G4InitXscPAI::G4InitXscPAI ( const G4MaterialCutsCouple *  matCC)

Definition at line 70 of file G4InitXscPAI.cc.

71 : fPAIxscVector(NULL),
72 fPAIdEdxVector(NULL),
73 fPAIphotonVector(NULL),
74 fPAIelectronVector(NULL),
75 fChCosSqVector(NULL),
76 fChWidthVector(NULL)
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 }
102 KillCloseIntervals();
103 Normalisation();
104 fBetaGammaSq = fTmax = 0.0;
105 fIntervalTmax = fCurrentInterval = 0;
106}
G4int
int G4int
Definition: G4Types.hh:66
G4InitXscPAI::KillCloseIntervals
void KillCloseIntervals()
Definition: G4InitXscPAI.cc:131
G4InitXscPAI::Normalisation
void Normalisation()
Definition: G4InitXscPAI.cc:162
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:150
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:179
G4Material::GetElectronDensity
G4double GetElectronDensity() const
Definition: G4Material.hh:216
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:261
G4SandiaTable::GetMaxInterval
G4int GetMaxInterval() const
Definition: G4SandiaTable.hh:189
G4SandiaTable::GetSandiaMatTable
G4double GetSandiaMatTable(G4int, G4int)
Definition: G4SandiaTable.cc:802

◆ ~G4InitXscPAI()

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

◆ DifPAIdEdx()

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
double G4double
Definition: G4Types.hh:64
G4InitXscPAI::DifPAIxSection
G4double DifPAIxSection(G4double omega)
Definition: G4InitXscPAI.cc:414

Referenced by IntegralPAIdEdx().

◆ DifPAIxSection()

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 G4double betaBohr2 = fine_structure_const*fine_structure_const ;
424 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
G4InitXscPAI::IntegralTerm
G4double IntegralTerm(G4double omega)
Definition: G4InitXscPAI.cc:256
G4InitXscPAI::ImPartDielectricConst
G4double ImPartDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4InitXscPAI.cc:294
G4InitXscPAI::RePartDielectricConst
G4double RePartDielectricConst(G4double energy)
Definition: G4InitXscPAI.cc:343
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:70

Referenced by DifPAIdEdx(), and IntegralPAIxSection().

◆ GetBinPAI()

G4int G4InitXscPAI::GetBinPAI ( ) const
inline

Definition at line 106 of file G4InitXscPAI.hh.

106{ return fPAIbin ; }

◆ GetChCosSqVector()

G4PhysicsLogVector * G4InitXscPAI::GetChCosSqVector ( ) const
inline

Definition at line 117 of file G4InitXscPAI.hh.

117{ return fChCosSqVector;}

◆ GetChWidthVector()

G4PhysicsLogVector * G4InitXscPAI::GetChWidthVector ( ) const
inline

Definition at line 118 of file G4InitXscPAI.hh.

118{ return fChWidthVector;}

◆ GetIntervalNumber()

G4int G4InitXscPAI::GetIntervalNumber ( ) const
inline

Definition at line 105 of file G4InitXscPAI.hh.

105{ return fIntervalNumber ; }

◆ GetMatSandiaMatrix()

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

Definition at line 110 of file G4InitXscPAI.hh.

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

◆ GetNormalizationCof()

G4double G4InitXscPAI::GetNormalizationCof ( ) const
inline

Definition at line 108 of file G4InitXscPAI.hh.

108{ return fNormalizationCof ; }

◆ GetPAIdEdxVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIdEdxVector ( ) const
inline

Definition at line 114 of file G4InitXscPAI.hh.

114{ return fPAIdEdxVector;}

◆ GetPAIelectronVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIelectronVector ( ) const
inline

Definition at line 116 of file G4InitXscPAI.hh.

116{ return fPAIelectronVector;}

◆ GetPAIphotonVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIphotonVector ( ) const
inline

Definition at line 115 of file G4InitXscPAI.hh.

115{ return fPAIphotonVector;}

◆ GetPAIxscVector()

G4PhysicsLogVector * G4InitXscPAI::GetPAIxscVector ( ) const
inline

Definition at line 113 of file G4InitXscPAI.hh.

113{ return fPAIxscVector;}

◆ GetPhotonLambda()

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}
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout

◆ GetStepCerenkovLoss()

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}

◆ GetStepEnergyLoss()

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}

◆ GetStepPlasmonLoss()

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}

◆ ImPartDielectricConst()

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

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

◆ IntegralCherenkov()

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
787 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
830 &G4InitXscPAI::PAIdNdxCherenkov,
831 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
832
833 else if( i == i1 ) result += integral.Legendre10(this,
834 &G4InitXscPAI::PAIdNdxCherenkov,energy1,
835 (*(*fMatSandiaMatrix)[i+1])[0]);
836
837 else result += integral.Legendre10(this,
838 &G4InitXscPAI::PAIdNdxCherenkov,
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
G4InitXscPAI::ModuleSqDielectricConst
G4double ModuleSqDielectricConst(G4int intervalNumber, G4double energy)
Definition: G4InitXscPAI.cc:320
G4InitXscPAI::PAIdNdxCherenkov
G4double PAIdNdxCherenkov(G4double omega)
Definition: G4InitXscPAI.cc:487
G4PhysicsVector::GetLowEdgeEnergy
virtual G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:158
G4PhysicsVector::PutValue
void PutValue(size_t index, G4double theValue)

◆ IntegralPAIdEdx()

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
699 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
734 &G4InitXscPAI::DifPAIdEdx,
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,
742 &G4InitXscPAI::DifPAIdEdx,
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}
G4InitXscPAI::DifPAIdEdx
G4double DifPAIdEdx(G4double omega)
Definition: G4InitXscPAI.cc:477

◆ IntegralPAIxSection()

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
618 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
653 &G4InitXscPAI::DifPAIxSection,
654 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
655
656 else if( i == i1 ) result += integral.Legendre10(this,
657 &G4InitXscPAI::DifPAIxSection,energy1,
658 (*(*fMatSandiaMatrix)[i+1])[0]);
659
660 else result += integral.Legendre10(this,
661 &G4InitXscPAI::DifPAIxSection,
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}

◆ IntegralPlasmon()

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
875 G4Integrator<G4InitXscPAI,G4double(G4InitXscPAI::*)(G4double)> integral;
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,
910 &G4InitXscPAI::PAIdNdxPlasmon,
911 (*(*fMatSandiaMatrix)[i])[0] ,energy2);
912
913 else if( i == i1 ) result += integral.Legendre10(this,
914 &G4InitXscPAI::PAIdNdxPlasmon,energy1,
915 (*(*fMatSandiaMatrix)[i+1])[0]);
916
917 else result += integral.Legendre10(this,
918 &G4InitXscPAI::PAIdNdxPlasmon,
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
G4InitXscPAI::PAIdNdxPlasmon
G4double PAIdNdxPlasmon(G4double omega)
Definition: G4InitXscPAI.cc:551

◆ IntegralTerm()

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}
G4InitXscPAI::RutherfordIntegral
G4double RutherfordIntegral(G4int intervalNumber, G4double limitLow, G4double limitHigh)
Definition: G4InitXscPAI.cc:232

Referenced by DifPAIxSection(), and PAIdNdxPlasmon().

◆ KillCloseIntervals()

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}

Referenced by G4InitXscPAI().

◆ ModuleSqDielectricConst()

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}

Referenced by IntegralCherenkov().

◆ Normalisation()

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}

Referenced by G4InitXscPAI().

◆ PAIdNdxCherenkov()

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, betaBohr2,betaBohr4,cofBetaBohr ;
496
497 //cof = 1.0 ;
498 cofBetaBohr = 4.0 ;
499 betaBohr2 = fine_structure_const*fine_structure_const ;
500 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

Referenced by IntegralCherenkov().

◆ PAIdNdxPlasmon()

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, betaBohr2, betaBohr4, cofBetaBohr ;
561
562 cof = 1 ;
563 cofBetaBohr = 4.0 ;
564 betaBohr2 = fine_structure_const*fine_structure_const ;
565 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

Referenced by IntegralPlasmon().

◆ RePartDielectricConst()

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

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

◆ RutherfordIntegral()

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

Referenced by IntegralTerm(), and Normalisation().

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