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 #include "G4StatMFFragment.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4HadronicException.hh"
00035
00036
00037 G4StatMFFragment::G4StatMFFragment(const G4StatMFFragment & )
00038 {
00039 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFFragment::copy_constructor meant to not be accessable");
00040 }
00041
00042
00043
00044 G4StatMFFragment & G4StatMFFragment::
00045 operator=(const G4StatMFFragment & )
00046 {
00047 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFFragment::operator= meant to not be accessable");
00048 return *this;
00049 }
00050
00051
00052 G4bool G4StatMFFragment::operator==(const G4StatMFFragment & ) const
00053 {
00054
00055 return false;
00056 }
00057
00058
00059 G4bool G4StatMFFragment::operator!=(const G4StatMFFragment & ) const
00060 {
00061
00062 return true;
00063 }
00064
00065
00066
00067 G4double G4StatMFFragment::GetCoulombEnergy(void) const
00068 {
00069 if (theZ <= 0.1) return 0.0;
00070 G4double Coulomb = (3./5.)*(elm_coupling*theZ*theZ)*
00071 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
00072 (G4StatMFParameters::Getr0()*std::pow(theA,1./3.));
00073
00074 return Coulomb;
00075 }
00076
00077
00078 G4double G4StatMFFragment::GetEnergy(const G4double T) const
00079 {
00080 if (theA < 1 || theZ < 0 || theZ > theA) {
00081 G4cerr << "G4StatMFFragment::GetEnergy: A = " << theA
00082 << ", Z = " << theZ << G4endl;
00083 throw G4HadronicException(__FILE__, __LINE__,
00084 "G4StatMFFragment::GetEnergy: Wrong values for A and Z!");
00085 }
00086 G4double BulkEnergy = G4NucleiProperties::GetMassExcess(static_cast<G4int>(theA),
00087 static_cast<G4int>(theZ));
00088
00089 if (theA < 4) return BulkEnergy - GetCoulombEnergy();
00090
00091 G4double SurfaceEnergy;
00092 if (G4StatMFParameters::DBetaDT(T) == 0.0) SurfaceEnergy = 0.0;
00093 else SurfaceEnergy = (5./2.)*std::pow(theA,2.0/3.0)*T*T*
00094 G4StatMFParameters::GetBeta0()/
00095 (G4StatMFParameters::GetCriticalTemp()*
00096 G4StatMFParameters::GetCriticalTemp());
00097
00098
00099 G4double ExchangeEnergy = theA*T*T/GetInvLevelDensity();
00100 if (theA != 4) ExchangeEnergy += SurfaceEnergy;
00101
00102 return BulkEnergy + ExchangeEnergy - GetCoulombEnergy();
00103
00104 }
00105
00106
00107 G4double G4StatMFFragment::GetInvLevelDensity(void) const
00108 {
00109
00110
00111 if (theA == 1) return 0.0;
00112 else return
00113 G4StatMFParameters::GetEpsilon0()*(1.0+3.0/(theA - 1.0));
00114 }
00115
00116
00117
00118 G4Fragment * G4StatMFFragment::GetFragment(const G4double T)
00119 {
00120 G4double U = CalcExcitationEnergy(T);
00121
00122 G4double M = GetNuclearMass();
00123
00124 G4LorentzVector FourMomentum(_momentum,std::sqrt(_momentum.mag2()+(M+U)*(M+U)));
00125
00126 G4Fragment * theFragment = new G4Fragment(static_cast<G4int>(theA),static_cast<G4int>(theZ),FourMomentum);
00127
00128 return theFragment;
00129 }
00130
00131
00132 G4double G4StatMFFragment::CalcExcitationEnergy(const G4double T)
00133 {
00134 if (theA <= 3) return 0.0;
00135
00136 G4double BulkEnergy = theA*T*T/GetInvLevelDensity();
00137
00138
00139 if (theA == 4) return BulkEnergy;
00140
00141
00142 G4double SurfaceEnergy = 0.0;
00143 if (std::abs(G4StatMFParameters::DBetaDT(T)) > 1.0e-20)
00144
00145
00146 SurfaceEnergy = (5./2.)*std::pow(theA,2.0/3.0)*(G4StatMFParameters::Beta(T) -
00147 T*G4StatMFParameters::DBetaDT(T) - G4StatMFParameters::GetBeta0());
00148
00149 return BulkEnergy + SurfaceEnergy;
00150 }
00151
00152