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 <numeric>
00033
00034 #include "G4StatMFMicroCanonical.hh"
00035 #include "G4PhysicalConstants.hh"
00036 #include "G4SystemOfUnits.hh"
00037 #include "G4HadronicException.hh"
00038 #include "G4Pow.hh"
00039
00040
00041 G4StatMFMicroCanonical::G4StatMFMicroCanonical(G4Fragment const & theFragment)
00042 {
00043
00044 Initialize(theFragment);
00045
00046 }
00047
00048
00049
00050 G4StatMFMicroCanonical::~G4StatMFMicroCanonical()
00051 {
00052
00053 if (!_ThePartitionManagerVector.empty()) {
00054 std::for_each(_ThePartitionManagerVector.begin(),
00055 _ThePartitionManagerVector.end(),
00056 DeleteFragment());
00057 }
00058 }
00059
00060
00061
00062
00063
00064 void G4StatMFMicroCanonical::Initialize(const G4Fragment & theFragment)
00065 {
00066
00067 std::vector<G4StatMFMicroManager*>::iterator it;
00068
00069
00070 G4double U = theFragment.GetExcitationEnergy();
00071
00072 G4int A = theFragment.GetA_asInt();
00073 G4int Z = theFragment.GetZ_asInt();
00074 G4double x = 1.0 - 2.0*Z/G4double(A);
00075 G4Pow* g4pow = G4Pow::GetInstance();
00076
00077
00078 G4double TConfiguration = std::sqrt(8.0*U/G4double(A));
00079
00080
00081 __FreeInternalE0 = A*(
00082
00083 -G4StatMFParameters::GetE0() +
00084
00085 G4StatMFParameters::GetGamma0()*x*x
00086 ) +
00087
00088 G4StatMFParameters::GetBeta0()*g4pow->Z23(A) +
00089
00090 elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*g4pow->Z13(A));
00091
00092
00093 G4double W = 0.0;
00094
00095
00096
00097 __MeanMultiplicity = 0.0;
00098
00099
00100 __MeanTemperature = 0.0;
00101
00102
00103 __MeanEntropy = 0.0;
00104
00105
00106 G4double SCompoundNucleus = CalcEntropyOfCompoundNucleus(theFragment,TConfiguration);
00107
00108
00109 _WCompoundNucleus = 1.0;
00110
00111 W += _WCompoundNucleus;
00112
00113
00114
00115
00116 G4int MaxMult = G4StatMFMicroCanonical::MaxAllowedMultiplicity;
00117 if (A > 110) MaxMult -= 1;
00118
00119
00120
00121 for (G4int im = 2; im <= MaxMult; im++) {
00122 G4StatMFMicroManager * aMicroManager =
00123 new G4StatMFMicroManager(theFragment,im,__FreeInternalE0,SCompoundNucleus);
00124 _ThePartitionManagerVector.push_back(aMicroManager);
00125 }
00126
00127
00128 W = std::accumulate(_ThePartitionManagerVector.begin(),
00129 _ThePartitionManagerVector.end(),
00130 W,SumProbabilities());
00131
00132
00133 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
00134 {
00135 (*it)->Normalize(W);
00136 }
00137
00138 _WCompoundNucleus /= W;
00139
00140 __MeanMultiplicity += 1.0 * _WCompoundNucleus;
00141 __MeanTemperature += TConfiguration * _WCompoundNucleus;
00142 __MeanEntropy += SCompoundNucleus * _WCompoundNucleus;
00143
00144
00145 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
00146 {
00147 __MeanMultiplicity += (*it)->GetMeanMultiplicity();
00148 __MeanTemperature += (*it)->GetMeanTemperature();
00149 __MeanEntropy += (*it)->GetMeanEntropy();
00150 }
00151
00152 return;
00153 }
00154
00155
00156 G4double G4StatMFMicroCanonical::CalcFreeInternalEnergy(const G4Fragment & theFragment,
00157 G4double T)
00158 {
00159 G4int A = theFragment.GetA_asInt();
00160 G4int Z = theFragment.GetZ_asInt();
00161 G4double A13 = G4Pow::GetInstance()->Z13(A);
00162
00163 G4double InvLevelDensityPar = G4StatMFParameters::GetEpsilon0()*(1.0 + 3.0/G4double(A-1));
00164
00165 G4double VolumeTerm = (-G4StatMFParameters::GetE0()+T*T/InvLevelDensityPar)*A;
00166
00167 G4double SymmetryTerm = G4StatMFParameters::GetGamma0()*(A - 2*Z)*(A - 2*Z)/G4double(A);
00168
00169 G4double SurfaceTerm = (G4StatMFParameters::Beta(T)-T*G4StatMFParameters::DBetaDT(T))*A13*A13;
00170
00171 G4double CoulombTerm = elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*A13);
00172
00173 return VolumeTerm + SymmetryTerm + SurfaceTerm + CoulombTerm;
00174 }
00175
00176 G4double
00177 G4StatMFMicroCanonical::CalcEntropyOfCompoundNucleus(const G4Fragment & theFragment,
00178 G4double & TConf)
00179
00180 {
00181 G4int A = theFragment.GetA_asInt();
00182
00183 G4double U = theFragment.GetExcitationEnergy();
00184 G4double A13 = G4Pow::GetInstance()->Z13(A);
00185
00186 G4double Ta = std::max(std::sqrt(U/(0.125*A)),0.0012*MeV);
00187 G4double Tb = Ta;
00188
00189 G4double ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Ta);
00190 G4double Da = (U+__FreeInternalE0-ECompoundNucleus)/U;
00191 G4double Db = 0.0;
00192
00193 G4double InvLevelDensity = CalcInvLevelDensity(static_cast<G4int>(A));
00194
00195
00196 if (Da == 0.0) {
00197 TConf = Ta;
00198 return 2*Ta*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Ta)*A13*A13;
00199 } else if (Da < 0.0) {
00200 do {
00201 Tb -= 0.5*Tb;
00202 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
00203 Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
00204 } while (Db < 0.0);
00205 } else {
00206 do {
00207 Tb += 0.5*Tb;
00208 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
00209 Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
00210 } while (Db > 0.0);
00211 }
00212
00213 G4double eps = 1.0e-14 * std::fabs(Tb-Ta);
00214
00215 for (G4int i = 0; i < 1000; i++) {
00216 G4double Tc = (Ta+Tb)/2.0;
00217 if (std::abs(Ta-Tb) <= eps) {
00218 TConf = Tc;
00219 return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
00220 }
00221 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tc);
00222 G4double Dc = (U+__FreeInternalE0-ECompoundNucleus)/U;
00223
00224 if (Dc == 0.0) {
00225 TConf = Tc;
00226 return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
00227 }
00228
00229 if (Da*Dc < 0.0) {
00230 Tb = Tc;
00231 Db = Dc;
00232 } else {
00233 Ta = Tc;
00234 Da = Dc;
00235 }
00236 }
00237
00238 G4cerr <<
00239 "G4StatMFMicrocanoncal::CalcEntropyOfCompoundNucleus: I can't calculate the temperature"
00240 << G4endl;
00241
00242 return 0.0;
00243 }
00244
00245 G4StatMFChannel * G4StatMFMicroCanonical::ChooseAandZ(const G4Fragment & theFragment)
00246
00247 {
00248
00249 G4double RandNumber = G4UniformRand();
00250
00251 if (RandNumber < _WCompoundNucleus) {
00252
00253 G4StatMFChannel * aChannel = new G4StatMFChannel;
00254 aChannel->CreateFragment(theFragment.GetA_asInt(),theFragment.GetZ_asInt());
00255 return aChannel;
00256
00257 } else {
00258
00259 G4double AccumWeight = _WCompoundNucleus;
00260 std::vector<G4StatMFMicroManager*>::iterator it;
00261 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it) {
00262 AccumWeight += (*it)->GetProbability();
00263 if (RandNumber < AccumWeight) {
00264 return (*it)->ChooseChannel(theFragment.GetA(),theFragment.GetZ(),__MeanTemperature);
00265 }
00266 }
00267 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::ChooseAandZ: wrong normalization!");
00268 }
00269
00270 return 0;
00271 }
00272
00273 G4double G4StatMFMicroCanonical::CalcInvLevelDensity(G4int anA)
00274 {
00275
00276
00277 if (anA == 1) return 0.0;
00278 else return
00279 G4StatMFParameters::GetEpsilon0()*(1.0+3.0/(anA - 1.0));
00280 }