00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4ProtonField.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4VNuclearDensity.hh"
00043 #include "G4FermiMomentum.hh"
00044 #include "G4V3DNucleus.hh"
00045
00046 G4ProtonField::G4ProtonField(G4V3DNucleus * aNucleus) :
00047 G4VNuclearField(aNucleus), theDensity(theNucleus->GetNuclearDensity())
00048 {
00049 theA = theNucleus->GetMassNumber();
00050 theZ = theNucleus->GetCharge();
00051 theBarrier = GetBarrier();
00052 theRadius = 2.*theNucleus->GetOuterRadius();
00053 theFermi.Init(theA, theZ);
00054 G4double aR=0;
00055 while(aR<theRadius)
00056 {
00057 G4ThreeVector aPosition(0,0,aR);
00058 G4double density = GetDensity(aPosition);
00059 G4double fermiMom = GetFermiMomentum(density);
00060 theFermiMomBuffer.push_back(fermiMom);
00061 aR+=0.3*fermi;
00062 }
00063 {
00064 G4ThreeVector aPosition(0,0,theRadius);
00065 G4double density = GetDensity(aPosition);
00066 G4double fermiMom = GetFermiMomentum(density);
00067 theFermiMomBuffer.push_back(fermiMom);
00068 }
00069 {
00070 G4ThreeVector aPosition(0,0,theRadius+0.001*fermi);
00071 theFermiMomBuffer.push_back(0);
00072 }
00073 {
00074 G4ThreeVector aPosition(0,0,1.*m);
00075 theFermiMomBuffer.push_back(0);
00076 }
00077 }
00078
00079
00080 G4ProtonField::~G4ProtonField()
00081 { }
00082
00083 G4double G4ProtonField::GetField(const G4ThreeVector & aPosition)
00084 {
00085
00086 G4double x = aPosition.mag();
00087 G4int index = static_cast<G4int>(x/(0.3*fermi) );
00088 if(index+2>static_cast<G4int>(theFermiMomBuffer.size())) return theFermiMomBuffer.back();
00089 G4double y1 = theFermiMomBuffer[index];
00090 G4double y2 = theFermiMomBuffer[index+1];
00091 G4double x1 = (0.3*fermi)*index;
00092 G4double x2 = (0.3*fermi)*(index+1);
00093 G4double fermiMom = y1 + (x-x1)*(y2-y1)/(x2-x1);
00094 G4double y = -1*(fermiMom*fermiMom)/(2*proton_mass_c2)+theBarrier;
00095
00096 return y;
00097 }
00098
00099 G4double G4ProtonField::GetBarrier()
00100 {
00101 G4double coulombBarrier = (1.44/1.14) * MeV * theZ / (1.0 + std::pow(theA,1./3.));
00102
00103 G4double bindingEnergy =0;
00104
00105
00106
00107
00108 return bindingEnergy/theA+coulombBarrier;
00109 }