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 "G4StatMFMicroPartition.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4HadronicException.hh"
00036
00037
00038 G4StatMFMicroPartition::G4StatMFMicroPartition(const G4StatMFMicroPartition & )
00039 {
00040 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::copy_constructor meant to not be accessable");
00041 }
00042
00043
00044
00045 G4StatMFMicroPartition & G4StatMFMicroPartition::
00046 operator=(const G4StatMFMicroPartition & )
00047 {
00048 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator= meant to not be accessable");
00049 return *this;
00050 }
00051
00052
00053 G4bool G4StatMFMicroPartition::operator==(const G4StatMFMicroPartition & ) const
00054 {
00055
00056 return false;
00057 }
00058
00059
00060 G4bool G4StatMFMicroPartition::operator!=(const G4StatMFMicroPartition & ) const
00061 {
00062
00063 return true;
00064 }
00065
00066
00067
00068 void G4StatMFMicroPartition::CoulombFreeEnergy(const G4double anA)
00069 {
00070
00071 G4double CoulombConstFactor = 1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
00072
00073 CoulombConstFactor = elm_coupling * (3./5.) *
00074 (1. - CoulombConstFactor)/G4StatMFParameters::Getr0();
00075
00076
00077
00078 if (anA == 0 || anA == 1)
00079 {
00080 _theCoulombFreeEnergy.push_back(CoulombConstFactor*(theZ/theA)*(theZ/theA));
00081 }
00082 else if (anA == 2 || anA == 3 || anA == 4)
00083 {
00084
00085 _theCoulombFreeEnergy.push_back(CoulombConstFactor*0.5*std::pow(anA,5./3.));
00086 }
00087 else
00088 {
00089 _theCoulombFreeEnergy.push_back(CoulombConstFactor*(theZ/theA)*(theZ/theA)*
00090 std::pow(anA,5./3.));
00091
00092 }
00093 }
00094
00095
00096 G4double G4StatMFMicroPartition::GetCoulombEnergy(void)
00097 {
00098 G4double CoulombFactor = 1.0/
00099 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
00100
00101 G4double CoulombEnergy = elm_coupling*(3./5.)*theZ*theZ*CoulombFactor/
00102 (G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.));
00103
00104 for (unsigned int i = 0; i < _thePartition.size(); i++)
00105 CoulombEnergy += _theCoulombFreeEnergy[i] - elm_coupling*(3./5.)*
00106 (theZ/theA)*(theZ/theA)*std::pow(static_cast<G4double>(_thePartition[i]),5./3.)/
00107 G4StatMFParameters::Getr0();
00108
00109 return CoulombEnergy;
00110 }
00111
00112 G4double G4StatMFMicroPartition::GetPartitionEnergy(const G4double T)
00113 {
00114 G4double CoulombFactor = 1.0/
00115 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
00116
00117 G4double PartitionEnergy = 0.0;
00118
00119
00120
00121 for (unsigned int i = 0; i < _thePartition.size(); i++)
00122 {
00123 if (_thePartition[i] == 0 || _thePartition[i] == 1)
00124 {
00125 PartitionEnergy += _theCoulombFreeEnergy[i];
00126 }
00127 else if (_thePartition[i] == 2)
00128 {
00129 PartitionEnergy +=
00130 -2.796
00131 + _theCoulombFreeEnergy[i];
00132 }
00133 else if (_thePartition[i] == 3)
00134 {
00135 PartitionEnergy +=
00136 -9.224
00137 + _theCoulombFreeEnergy[i];
00138 }
00139 else if (_thePartition[i] == 4)
00140 {
00141 PartitionEnergy +=
00142 -30.11
00143 + _theCoulombFreeEnergy[i]
00144 + 4.*T*T/InvLevelDensity(4.);
00145 }
00146 else
00147 {
00148 PartitionEnergy +=
00149
00150 (- G4StatMFParameters::GetE0() +
00151 T*T/InvLevelDensity(_thePartition[i]))
00152 *_thePartition[i] +
00153
00154
00155 G4StatMFParameters::GetGamma0()*
00156 (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/theA)*_thePartition[i] +
00157
00158
00159 (G4StatMFParameters::Beta(T) - T*G4StatMFParameters::DBetaDT(T))*
00160 std::pow(static_cast<G4double>(_thePartition[i]),2./3.) +
00161
00162
00163 _theCoulombFreeEnergy[i];
00164 }
00165 }
00166
00167 PartitionEnergy += elm_coupling*(3./5.)*theZ*theZ*CoulombFactor/
00168 (G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.))
00169 + (3./2.)*T*(_thePartition.size()-1);
00170
00171 return PartitionEnergy;
00172 }
00173
00174
00175 G4double G4StatMFMicroPartition::CalcPartitionTemperature(const G4double U,
00176 const G4double FreeInternalE0)
00177 {
00178 G4double PartitionEnergy = GetPartitionEnergy(0.0);
00179
00180
00181
00182 if (std::abs(U + FreeInternalE0 - PartitionEnergy) < 0.003) return -1.0;
00183
00184
00185
00186
00187 G4double Ta = 0.001;
00188 G4double Tb = std::max(std::sqrt(8.0*U/theA),0.0012*MeV);
00189 G4double Tmid = 0.0;
00190
00191 G4double Da = (U + FreeInternalE0 - GetPartitionEnergy(Ta))/U;
00192 G4double Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
00193
00194 G4int maxit = 0;
00195 while (Da*Db > 0.0 && maxit < 1000)
00196 {
00197 ++maxit;
00198 Tb += 0.5*Tb;
00199 Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
00200 }
00201
00202 G4double eps = 1.0e-14*std::abs(Ta-Tb);
00203
00204 for (G4int i = 0; i < 1000; i++)
00205 {
00206 Tmid = (Ta+Tb)/2.0;
00207 if (std::abs(Ta-Tb) <= eps) return Tmid;
00208 G4double Dmid = (U + FreeInternalE0 - GetPartitionEnergy(Tmid))/U;
00209 if (std::abs(Dmid) < 0.003) return Tmid;
00210 if (Da*Dmid < 0.0)
00211 {
00212 Tb = Tmid;
00213 Db = Dmid;
00214 }
00215 else
00216 {
00217 Ta = Tmid;
00218 Da = Dmid;
00219 }
00220 }
00221
00222 G4cerr << "G4StatMFMicroPartition::CalcPartitionTemperature: I can't calculate the temperature"
00223 << G4endl;
00224
00225 return -1.0;
00226
00227 }
00228
00229
00230 G4double G4StatMFMicroPartition::CalcPartitionProbability(const G4double U,
00231 const G4double FreeInternalE0,
00232 const G4double SCompound)
00233 {
00234 G4double T = CalcPartitionTemperature(U,FreeInternalE0);
00235 if ( T <= 0.0) return _Probability = 0.0;
00236 _Temperature = T;
00237
00238
00239
00240 G4double Fact = 1.0;
00241 unsigned int i;
00242 for (i = 0; i < _thePartition.size() - 1; i++)
00243 {
00244 G4double f = 1.0;
00245 for (unsigned int ii = i+1; i< _thePartition.size(); i++)
00246 {
00247 if (_thePartition[i] == _thePartition[ii]) f++;
00248 }
00249 Fact *= f;
00250 }
00251
00252 G4double ProbDegeneracy = 1.0;
00253 G4double ProbA32 = 1.0;
00254
00255 for (i = 0; i < _thePartition.size(); i++)
00256 {
00257 ProbDegeneracy *= GetDegeneracyFactor(static_cast<G4int>(_thePartition[i]));
00258 ProbA32 *= static_cast<G4double>(_thePartition[i])*
00259 std::sqrt(static_cast<G4double>(_thePartition[i]));
00260 }
00261
00262
00263 G4double PartitionEntropy = 0.0;
00264 for (i = 0; i < _thePartition.size(); i++)
00265 {
00266
00267 if (_thePartition[i] == 4)
00268 {
00269 PartitionEntropy +=
00270 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i]);
00271 }
00272
00273 else if (_thePartition[i] > 4)
00274 {
00275 PartitionEntropy +=
00276 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i])
00277 - G4StatMFParameters::DBetaDT(T)
00278 * std::pow(static_cast<G4double>(_thePartition[i]),2.0/3.0);
00279 }
00280 }
00281
00282
00283 G4double ThermalWaveLenght3 = 16.15*fermi/std::sqrt(T);
00284 ThermalWaveLenght3 = ThermalWaveLenght3*ThermalWaveLenght3*ThermalWaveLenght3;
00285
00286
00287 G4double kappa = (1. + elm_coupling*(std::pow(static_cast<G4double>(_thePartition.size()),1./3.)-1.0)
00288 /(G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.)));
00289 kappa = kappa*kappa*kappa;
00290 kappa -= 1.;
00291 G4double V0 = (4./3.)*pi*theA*G4StatMFParameters::Getr0()*G4StatMFParameters::Getr0()*
00292 G4StatMFParameters::Getr0();
00293 G4double FreeVolume = kappa*V0;
00294 G4double TranslationalS = std::max(0.0, std::log(ProbA32/Fact) +
00295 (_thePartition.size()-1.0)*std::log(FreeVolume/ThermalWaveLenght3) +
00296 1.5*(_thePartition.size()-1.0) - (3./2.)*std::log(theA));
00297
00298 PartitionEntropy += std::log(ProbDegeneracy) + TranslationalS;
00299 _Entropy = PartitionEntropy;
00300
00301
00302 G4double exponent = PartitionEntropy-SCompound;
00303 if (exponent > 700.0) exponent = 700.0;
00304 return _Probability = std::exp(exponent);
00305 }
00306
00307
00308
00309 G4double G4StatMFMicroPartition::GetDegeneracyFactor(const G4int A)
00310 {
00311
00312
00313 G4double DegFactor = 0;
00314 if (A > 4) DegFactor = 1.0;
00315 else if (A == 1) DegFactor = 4.0;
00316 else if (A == 2) DegFactor = 3.0;
00317 else if (A == 3) DegFactor = 2.0+2.0;
00318 else if (A == 4) DegFactor = 1.0;
00319 return DegFactor;
00320 }
00321
00322
00323 G4StatMFChannel * G4StatMFMicroPartition::ChooseZ(const G4double A0, const G4double Z0, const G4double MeanT)
00324
00325 {
00326 std::vector<G4int> FragmentsZ;
00327
00328 G4int ZBalance = 0;
00329 do
00330 {
00331 G4double CC = G4StatMFParameters::GetGamma0()*8.0;
00332 G4int SumZ = 0;
00333 for (unsigned int i = 0; i < _thePartition.size(); i++)
00334 {
00335 G4double ZMean;
00336 G4double Af = _thePartition[i];
00337 if (Af > 1.5 && Af < 4.5) ZMean = 0.5*Af;
00338 else ZMean = Af*Z0/A0;
00339 G4double ZDispersion = std::sqrt(Af * MeanT/CC);
00340 G4int Zf;
00341 do
00342 {
00343 Zf = static_cast<G4int>(G4RandGauss::shoot(ZMean,ZDispersion));
00344 }
00345 while (Zf < 0 || Zf > Af);
00346 FragmentsZ.push_back(Zf);
00347 SumZ += Zf;
00348 }
00349 ZBalance = static_cast<G4int>(Z0) - SumZ;
00350 }
00351 while (std::abs(ZBalance) > 1.1);
00352 FragmentsZ[0] += ZBalance;
00353
00354 G4StatMFChannel * theChannel = new G4StatMFChannel;
00355 for (unsigned int i = 0; i < _thePartition.size(); i++)
00356 {
00357 theChannel->CreateFragment(_thePartition[i],FragmentsZ[i]);
00358 }
00359
00360 return theChannel;
00361 }