source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroPartition.cc@ 1199

Last change on this file since 1199 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 11.9 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4StatMFMicroPartition.cc,v 1.8 2008/07/25 11:20:47 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// by V. Lara
31// --------------------------------------------------------------------
32
33#include "G4StatMFMicroPartition.hh"
34#include "G4HadronicException.hh"
35
36
37// Copy constructor
38G4StatMFMicroPartition::G4StatMFMicroPartition(const G4StatMFMicroPartition & )
39{
40 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::copy_constructor meant to not be accessable");
41}
42
43// Operators
44
45G4StatMFMicroPartition & G4StatMFMicroPartition::
46operator=(const G4StatMFMicroPartition & )
47{
48 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator= meant to not be accessable");
49 return *this;
50}
51
52
53G4bool G4StatMFMicroPartition::operator==(const G4StatMFMicroPartition & ) const
54{
55 //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator== meant to not be accessable");
56 return false;
57}
58
59
60G4bool G4StatMFMicroPartition::operator!=(const G4StatMFMicroPartition & ) const
61{
62 //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator!= meant to not be accessable");
63 return true;
64}
65
66
67
68void G4StatMFMicroPartition::CoulombFreeEnergy(const G4double anA)
69{
70 // This Z independent factor in the Coulomb free energy
71 G4double CoulombConstFactor = 1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
72
73 CoulombConstFactor = elm_coupling * (3./5.) *
74 (1. - CoulombConstFactor)/G4StatMFParameters::Getr0();
75
76 // We use the aproximation Z_f ~ Z/A * A_f
77
78 if (anA == 0 || anA == 1)
79 {
80 _theCoulombFreeEnergy.push_back(CoulombConstFactor*(theZ/theA)*(theZ/theA));
81 }
82 else if (anA == 2 || anA == 3 || anA == 4)
83 {
84 // Z/A ~ 1/2
85 _theCoulombFreeEnergy.push_back(CoulombConstFactor*0.5*std::pow(anA,5./3.));
86 }
87 else // anA > 4
88 {
89 _theCoulombFreeEnergy.push_back(CoulombConstFactor*(theZ/theA)*(theZ/theA)*
90 std::pow(anA,5./3.));
91
92 }
93}
94
95
96G4double G4StatMFMicroPartition::GetCoulombEnergy(void)
97{
98 G4double CoulombFactor = 1.0/
99 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
100
101 G4double CoulombEnergy = elm_coupling*(3./5.)*theZ*theZ*CoulombFactor/
102 (G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.));
103
104 for (unsigned int i = 0; i < _thePartition.size(); i++)
105 CoulombEnergy += _theCoulombFreeEnergy[i] - elm_coupling*(3./5.)*
106 (theZ/theA)*(theZ/theA)*std::pow(static_cast<G4double>(_thePartition[i]),5./3.)/
107 G4StatMFParameters::Getr0();
108
109 return CoulombEnergy;
110}
111
112G4double G4StatMFMicroPartition::GetPartitionEnergy(const G4double T)
113{
114 G4double CoulombFactor = 1.0/
115 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0);
116
117 G4double PartitionEnergy = 0.0;
118
119
120 // We use the aprox that Z_f ~ Z/A * A_f
121 for (unsigned int i = 0; i < _thePartition.size(); i++)
122 {
123 if (_thePartition[i] == 0 || _thePartition[i] == 1)
124 {
125 PartitionEnergy += _theCoulombFreeEnergy[i];
126 }
127 else if (_thePartition[i] == 2)
128 {
129 PartitionEnergy +=
130 -2.796 // Binding Energy of deuteron ??????
131 + _theCoulombFreeEnergy[i];
132 }
133 else if (_thePartition[i] == 3)
134 {
135 PartitionEnergy +=
136 -9.224 // Binding Energy of trtion/He3 ??????
137 + _theCoulombFreeEnergy[i];
138 }
139 else if (_thePartition[i] == 4)
140 {
141 PartitionEnergy +=
142 -30.11 // Binding Energy of ALPHA ??????
143 + _theCoulombFreeEnergy[i]
144 + 4.*T*T/InvLevelDensity(4.);
145 }
146 else
147 {
148 PartitionEnergy +=
149 //Volume term
150 (- G4StatMFParameters::GetE0() +
151 T*T/InvLevelDensity(_thePartition[i]))
152 *_thePartition[i] +
153
154 // Symmetry term
155 G4StatMFParameters::GetGamma0()*
156 (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/theA)*_thePartition[i] +
157
158 // Surface term
159 (G4StatMFParameters::Beta(T) - T*G4StatMFParameters::DBetaDT(T))*
160 std::pow(static_cast<G4double>(_thePartition[i]),2./3.) +
161
162 // Coulomb term
163 _theCoulombFreeEnergy[i];
164 }
165 }
166
167 PartitionEnergy += elm_coupling*(3./5.)*theZ*theZ*CoulombFactor/
168 (G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.))
169 + (3./2.)*T*(_thePartition.size()-1);
170
171 return PartitionEnergy;
172}
173
174
175G4double G4StatMFMicroPartition::CalcPartitionTemperature(const G4double U,
176 const G4double FreeInternalE0)
177{
178 G4double PartitionEnergy = GetPartitionEnergy(0.0);
179
180 // If this happens, T = 0 MeV, which means that probability for this
181 // partition will be 0
182 if (std::abs(U + FreeInternalE0 - PartitionEnergy) < 0.003) return -1.0;
183
184 // Calculate temperature by midpoint method
185
186 // Bracketing the solution
187 G4double Ta = 0.001;
188 G4double Tb = std::max(std::sqrt(8.0*U/theA),0.0012*MeV);
189 G4double Tmid = 0.0;
190
191 G4double Da = (U + FreeInternalE0 - GetPartitionEnergy(Ta))/U;
192 G4double Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
193
194 G4int maxit = 0;
195 while (Da*Db > 0.0 && maxit < 1000)
196 {
197 ++maxit;
198 Tb += 0.5*Tb;
199 Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
200 }
201
202 G4double eps = 1.0e-14*std::abs(Ta-Tb);
203
204 for (G4int i = 0; i < 1000; i++)
205 {
206 Tmid = (Ta+Tb)/2.0;
207 if (std::abs(Ta-Tb) <= eps) return Tmid;
208 G4double Dmid = (U + FreeInternalE0 - GetPartitionEnergy(Tmid))/U;
209 if (std::abs(Dmid) < 0.003) return Tmid;
210 if (Da*Dmid < 0.0)
211 {
212 Tb = Tmid;
213 Db = Dmid;
214 }
215 else
216 {
217 Ta = Tmid;
218 Da = Dmid;
219 }
220 }
221 // if we arrive here the temperature could not be calculated
222 G4cerr << "G4StatMFMicroPartition::CalcPartitionTemperature: I can't calculate the temperature"
223 << G4endl;
224 // and set probability to 0 returning T < 0
225 return -1.0;
226
227}
228
229
230G4double G4StatMFMicroPartition::CalcPartitionProbability(const G4double U,
231 const G4double FreeInternalE0,
232 const G4double SCompound)
233{
234 G4double T = CalcPartitionTemperature(U,FreeInternalE0);
235 if ( T <= 0.0) return _Probability = 0.0;
236 _Temperature = T;
237
238
239 // Factorial of fragment multiplicity
240 G4double Fact = 1.0;
241 unsigned int i;
242 for (i = 0; i < _thePartition.size() - 1; i++)
243 {
244 G4double f = 1.0;
245 for (unsigned int ii = i+1; i< _thePartition.size(); i++)
246 {
247 if (_thePartition[i] == _thePartition[ii]) f++;
248 }
249 Fact *= f;
250 }
251
252 G4double ProbDegeneracy = 1.0;
253 G4double ProbA32 = 1.0;
254
255 for (i = 0; i < _thePartition.size(); i++)
256 {
257 ProbDegeneracy *= GetDegeneracyFactor(static_cast<G4int>(_thePartition[i]));
258 ProbA32 *= static_cast<G4double>(_thePartition[i])*
259 std::sqrt(static_cast<G4double>(_thePartition[i]));
260 }
261
262 // Compute entropy
263 G4double PartitionEntropy = 0.0;
264 for (i = 0; i < _thePartition.size(); i++)
265 {
266 // interaction entropy for alpha
267 if (_thePartition[i] == 4)
268 {
269 PartitionEntropy +=
270 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i]);
271 }
272 // interaction entropy for Af > 4
273 else if (_thePartition[i] > 4)
274 {
275 PartitionEntropy +=
276 2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i])
277 - G4StatMFParameters::DBetaDT(T)
278 * std::pow(static_cast<G4double>(_thePartition[i]),2.0/3.0);
279 }
280 }
281
282 // Thermal Wave Lenght = std::sqrt(2 pi hbar^2 / nucleon_mass T)
283 G4double ThermalWaveLenght3 = 16.15*fermi/std::sqrt(T);
284 ThermalWaveLenght3 = ThermalWaveLenght3*ThermalWaveLenght3*ThermalWaveLenght3;
285
286 // Translational Entropy
287 G4double kappa = (1. + elm_coupling*(std::pow(static_cast<G4double>(_thePartition.size()),1./3.)-1.0)
288 /(G4StatMFParameters::Getr0()*std::pow(static_cast<G4double>(theA),1./3.)));
289 kappa = kappa*kappa*kappa;
290 kappa -= 1.;
291 G4double V0 = (4./3.)*pi*theA*G4StatMFParameters::Getr0()*G4StatMFParameters::Getr0()*
292 G4StatMFParameters::Getr0();
293 G4double FreeVolume = kappa*V0;
294 G4double TranslationalS = std::max(0.0, std::log(ProbA32/Fact) +
295 (_thePartition.size()-1.0)*std::log(FreeVolume/ThermalWaveLenght3) +
296 1.5*(_thePartition.size()-1.0) - (3./2.)*std::log(theA));
297
298 PartitionEntropy += std::log(ProbDegeneracy) + TranslationalS;
299 _Entropy = PartitionEntropy;
300
301 // And finally compute probability of fragment configuration
302 G4double exponent = PartitionEntropy-SCompound;
303 if (exponent > 700.0) exponent = 700.0;
304 return _Probability = std::exp(exponent);
305}
306
307
308
309G4double G4StatMFMicroPartition::GetDegeneracyFactor(const G4int A)
310{
311 // Degeneracy factors are statistical factors
312 // DegeneracyFactor for nucleon is (2S_n + 1)(2I_n + 1) = 4
313 G4double DegFactor = 0;
314 if (A > 4) DegFactor = 1.0;
315 else if (A == 1) DegFactor = 4.0; // nucleon
316 else if (A == 2) DegFactor = 3.0; // Deuteron
317 else if (A == 3) DegFactor = 2.0+2.0; // Triton + He3
318 else if (A == 4) DegFactor = 1.0; // alpha
319 return DegFactor;
320}
321
322
323G4StatMFChannel * G4StatMFMicroPartition::ChooseZ(const G4double A0, const G4double Z0, const G4double MeanT)
324// Gives fragments charges
325{
326 std::vector<G4int> FragmentsZ;
327
328 G4int ZBalance = 0;
329 do
330 {
331 G4double CC = G4StatMFParameters::GetGamma0()*8.0;
332 G4int SumZ = 0;
333 for (unsigned int i = 0; i < _thePartition.size(); i++)
334 {
335 G4double ZMean;
336 G4double Af = _thePartition[i];
337 if (Af > 1.5 && Af < 4.5) ZMean = 0.5*Af;
338 else ZMean = Af*Z0/A0;
339 G4double ZDispersion = std::sqrt(Af * MeanT/CC);
340 G4int Zf;
341 do
342 {
343 Zf = static_cast<G4int>(G4RandGauss::shoot(ZMean,ZDispersion));
344 }
345 while (Zf < 0 || Zf > Af);
346 FragmentsZ.push_back(Zf);
347 SumZ += Zf;
348 }
349 ZBalance = static_cast<G4int>(Z0) - SumZ;
350 }
351 while (std::abs(ZBalance) > 1.1);
352 FragmentsZ[0] += ZBalance;
353
354 G4StatMFChannel * theChannel = new G4StatMFChannel;
355 for (unsigned int i = 0; i < _thePartition.size(); i++)
356 {
357 theChannel->CreateFragment(_thePartition[i],FragmentsZ[i]);
358 }
359
360 return theChannel;
361}
Note: See TracBrowser for help on using the repository browser.