40 throw G4HadronicException(__FILE__, __LINE__,
"G4StatMFMicroPartition::copy_constructor meant to not be accessable");
48 throw G4HadronicException(__FILE__, __LINE__,
"G4StatMFMicroPartition::operator= meant to not be accessable");
66 void G4StatMFMicroPartition::CoulombFreeEnergy(
G4int anA)
78 if (anA == 0 || anA == 1)
80 _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA);
82 else if (anA == 2 || anA == 3 || anA == 4)
85 _theCoulombFreeEnergy.push_back(CoulombConstFactor*0.5*std::pow(anA,5./3.));
89 _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA*std::pow(anA,5./3.));
93 G4double G4StatMFMicroPartition::GetCoulombEnergy(
void)
102 for (
unsigned int i = 0; i < _thePartition.size(); i++)
103 CoulombEnergy += _theCoulombFreeEnergy[i] -
elm_coupling*(3./5.)*
104 ZA*ZA*std::pow(static_cast<G4double>(_thePartition[i]),5./3.)/
107 return CoulombEnergy;
118 for (
unsigned int i = 0; i < _thePartition.size(); i++)
120 if (_thePartition[i] == 0 || _thePartition[i] == 1)
122 PartitionEnergy += _theCoulombFreeEnergy[i];
124 else if (_thePartition[i] == 2)
128 + _theCoulombFreeEnergy[i];
130 else if (_thePartition[i] == 3)
134 + _theCoulombFreeEnergy[i];
136 else if (_thePartition[i] == 4)
140 + _theCoulombFreeEnergy[i]
141 + 4.*T*T/InvLevelDensity(4.);
148 T*T/InvLevelDensity(_thePartition[i]))
153 (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/theA)*_thePartition[i] +
157 std::pow(static_cast<G4double>(_thePartition[i]),2./3.) +
160 _theCoulombFreeEnergy[i];
164 PartitionEnergy +=
elm_coupling*(3./5.)*theZ*theZ*CoulombFactor/
166 + (3./2.)*T*(_thePartition.size()-1);
168 return PartitionEnergy;
174 G4double PartitionEnergy = GetPartitionEnergy(0.0);
178 if (std::fabs(U + FreeInternalE0 - PartitionEnergy) < 0.003)
return -1.0;
187 G4double Da = (U + FreeInternalE0 - GetPartitionEnergy(Ta))/U;
188 G4double Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
191 while (Da*Db > 0.0 && maxit < 1000)
195 Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
198 G4double eps = 1.0e-14*std::abs(Ta-Tb);
200 for (
G4int i = 0; i < 1000; i++)
203 if (std::fabs(Ta-Tb) <= eps)
return Tmid;
204 G4double Dmid = (U + FreeInternalE0 - GetPartitionEnergy(Tmid))/U;
205 if (std::fabs(Dmid) < 0.003)
return Tmid;
218 G4cerr <<
"G4StatMFMicroPartition::CalcPartitionTemperature: I can't calculate the temperature"
229 G4double T = CalcPartitionTemperature(U,FreeInternalE0);
230 if ( T <= 0.0)
return _Probability = 0.0;
237 for (i = 0; i < _thePartition.size() - 1; i++)
240 for (
unsigned int ii = i+1; i< _thePartition.size(); i++)
242 if (_thePartition[i] == _thePartition[ii]) f++;
250 for (i = 0; i < _thePartition.size(); i++)
252 ProbDegeneracy *= GetDegeneracyFactor(static_cast<G4int>(_thePartition[i]));
253 ProbA32 *=
static_cast<G4double>(_thePartition[i])*
254 std::sqrt(static_cast<G4double>(_thePartition[i]));
259 for (i = 0; i < _thePartition.size(); i++)
262 if (_thePartition[i] == 4)
265 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i]);
268 else if (_thePartition[i] > 4)
271 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i])
273 * std::pow(static_cast<G4double>(_thePartition[i]),2.0/3.0);
279 ThermalWaveLenght3 = ThermalWaveLenght3*ThermalWaveLenght3*ThermalWaveLenght3;
282 G4double kappa = (1. +
elm_coupling*(std::pow(static_cast<G4double>(_thePartition.size()),1./3.)-1.0)
284 kappa = kappa*kappa*kappa;
290 (_thePartition.size()-1.0)*std::log(FreeVolume/ThermalWaveLenght3) +
291 1.5*(_thePartition.size()-1.0) - (3./2.)*std::log(
G4double(theA)));
293 PartitionEntropy += std::log(ProbDegeneracy) + TranslationalS;
294 _Entropy = PartitionEntropy;
297 G4double exponent = PartitionEntropy-SCompound;
298 if (exponent > 700.0) exponent = 700.0;
299 return _Probability = std::exp(exponent);
302 G4double G4StatMFMicroPartition::GetDegeneracyFactor(
G4int A)
307 if (A > 4) DegFactor = 1.0;
308 else if (A == 1) DegFactor = 4.0;
309 else if (A == 2) DegFactor = 3.0;
310 else if (A == 3) DegFactor = 4.0;
311 else if (A == 4) DegFactor = 1.0;
318 std::vector<G4int> FragmentsZ;
325 for (
unsigned int i = 0; i < _thePartition.size(); i++)
329 if (Af > 1.5 && Af < 4.5) ZMean = 0.5*Af;
330 else ZMean = Af*Z0/A0;
331 G4double ZDispersion = std::sqrt(Af * MeanT/CC);
337 while (Zf < 0 || Zf > Af);
338 FragmentsZ.push_back(Zf);
341 ZBalance = Z0 - SumZ;
343 while (std::abs(ZBalance) > 1);
344 FragmentsZ[0] += ZBalance;
347 for (
unsigned int i = 0; i < _thePartition.size(); i++)
static G4double GetGamma0()
G4bool operator!=(const G4StatMFMicroPartition &right) const
ThreeVector shoot(const G4int Ap, const G4int Af)
static G4double GetKappaCoulomb()
G4bool operator==(const G4StatMFMicroPartition &right) const
void CreateFragment(G4int A, G4int Z)
G4StatMFChannel * ChooseZ(G4int A0, G4int Z0, G4double MeanT)
G4StatMFMicroPartition(G4int A, G4int Z)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double CalcPartitionProbability(G4double U, G4double FreeInternalE0, G4double SCompound)
static G4double DBetaDT(G4double T)
static G4double Beta(G4double T)
G4GLOB_DLL std::ostream G4cerr