source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroCanonical.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: 9.6 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: G4StatMFMicroCanonical.cc,v 1.7 2008/07/25 11:20:47 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara
32
33
34#include "G4StatMFMicroCanonical.hh"
35#include "G4HadronicException.hh"
36#include <numeric>
37
38
39// Copy constructor
40G4StatMFMicroCanonical::G4StatMFMicroCanonical(const G4StatMFMicroCanonical &
41 ) : G4VStatMFEnsemble()
42{
43 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::copy_constructor meant to not be accessable");
44}
45
46// Operators
47
48G4StatMFMicroCanonical & G4StatMFMicroCanonical::
49operator=(const G4StatMFMicroCanonical & )
50{
51 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::operator= meant to not be accessable");
52 return *this;
53}
54
55
56G4bool G4StatMFMicroCanonical::operator==(const G4StatMFMicroCanonical & ) const
57{
58 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::operator== meant to not be accessable");
59 return false;
60}
61
62
63G4bool G4StatMFMicroCanonical::operator!=(const G4StatMFMicroCanonical & ) const
64{
65 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::operator!= meant to not be accessable");
66 return true;
67}
68
69
70
71// constructor
72G4StatMFMicroCanonical::G4StatMFMicroCanonical(G4Fragment const & theFragment)
73{
74 // Perform class initialization
75 Initialize(theFragment);
76
77}
78
79
80// destructor
81G4StatMFMicroCanonical::~G4StatMFMicroCanonical()
82{
83 // garbage collection
84 if (!_ThePartitionManagerVector.empty()) {
85 std::for_each(_ThePartitionManagerVector.begin(),
86 _ThePartitionManagerVector.end(),
87 DeleteFragment());
88 }
89}
90
91
92
93// Initialization method
94
95void G4StatMFMicroCanonical::Initialize(const G4Fragment & theFragment)
96{
97
98 std::vector<G4StatMFMicroManager*>::iterator it;
99
100 // Excitation Energy
101 G4double U = theFragment.GetExcitationEnergy();
102
103 G4double A = theFragment.GetA();
104 G4double Z = theFragment.GetZ();
105
106 // Configuration temperature
107 G4double TConfiguration = std::sqrt(8.0*U/A);
108
109 // Free internal energy at Temperature T = 0
110 __FreeInternalE0 = A*(
111 // Volume term (for T = 0)
112 -G4StatMFParameters::GetE0() +
113 // Symmetry term
114 G4StatMFParameters::GetGamma0()*(1.0-2.0*Z/A)*(1.0-2.0*Z/A)
115 ) +
116 // Surface term (for T = 0)
117 G4StatMFParameters::GetBeta0()*std::pow(A,2.0/3.0) +
118 // Coulomb term
119 elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*std::pow(A,1.0/3.0));
120
121 // Statistical weight
122 G4double W = 0.0;
123
124
125 // Mean breakup multiplicity
126 __MeanMultiplicity = 0.0;
127
128 // Mean channel temperature
129 __MeanTemperature = 0.0;
130
131 // Mean channel entropy
132 __MeanEntropy = 0.0;
133
134 // Calculate entropy of compound nucleus
135 G4double SCompoundNucleus = CalcEntropyOfCompoundNucleus(theFragment,TConfiguration);
136
137 // Statistical weight of compound nucleus
138 _WCompoundNucleus = 1.0; // std::exp(SCompoundNucleus - SCompoundNucleus);
139
140 W += _WCompoundNucleus;
141
142
143
144 // Maximal fragment multiplicity allowed in direct simulation
145 G4int MaxMult = G4StatMFMicroCanonical::MaxAllowedMultiplicity;
146 if (A > 110) MaxMult -= 1;
147
148
149
150 for (G4int m = 2; m <= MaxMult; m++) {
151 G4StatMFMicroManager * aMicroManager =
152 new G4StatMFMicroManager(theFragment,m,__FreeInternalE0,SCompoundNucleus);
153 _ThePartitionManagerVector.push_back(aMicroManager);
154 }
155
156 // W is the total probability
157 W = std::accumulate(_ThePartitionManagerVector.begin(),
158 _ThePartitionManagerVector.end(),
159 W,SumProbabilities());
160
161 // Normalization of statistical weights
162 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
163 {
164 (*it)->Normalize(W);
165 }
166
167 _WCompoundNucleus /= W;
168
169 __MeanMultiplicity += 1.0 * _WCompoundNucleus;
170 __MeanTemperature += TConfiguration * _WCompoundNucleus;
171 __MeanEntropy += SCompoundNucleus * _WCompoundNucleus;
172
173
174 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
175 {
176 __MeanMultiplicity += (*it)->GetMeanMultiplicity();
177 __MeanTemperature += (*it)->GetMeanTemperature();
178 __MeanEntropy += (*it)->GetMeanEntropy();
179 }
180
181 return;
182}
183
184
185
186
187
188G4double G4StatMFMicroCanonical::CalcFreeInternalEnergy(const G4Fragment & theFragment,
189 const G4double T)
190{
191 G4double A = theFragment.GetA();
192 G4double Z = theFragment.GetZ();
193 G4double A13 = std::pow(A,1.0/3.0);
194
195 G4double InvLevelDensityPar = G4StatMFParameters::GetEpsilon0()*(1.0 + 3.0/(A-1.0));
196
197 G4double VolumeTerm = (-G4StatMFParameters::GetE0()+T*T/InvLevelDensityPar)*A;
198
199 G4double SymmetryTerm = G4StatMFParameters::GetGamma0()*(A - 2.0*Z)*(A - 2.0*Z)/A;
200
201 G4double SurfaceTerm = (G4StatMFParameters::Beta(T)-T*G4StatMFParameters::DBetaDT(T))*A13*A13;
202
203 G4double CoulombTerm = elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*A13);
204
205 return VolumeTerm + SymmetryTerm + SurfaceTerm + CoulombTerm;
206}
207
208
209
210G4double G4StatMFMicroCanonical::CalcEntropyOfCompoundNucleus(const G4Fragment & theFragment,G4double & TConf)
211 // Calculates Temperature and Entropy of compound nucleus
212{
213 const G4double A = theFragment.GetA();
214 // const G4double Z = theFragment.GetZ();
215 const G4double U = theFragment.GetExcitationEnergy();
216 const G4double A13 = std::pow(A,1.0/3.0);
217
218 G4double Ta = std::max(std::sqrt(U/(0.125*A)),0.0012*MeV);
219 G4double Tb = Ta;
220
221 G4double ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Ta);
222 G4double Da = (U+__FreeInternalE0-ECompoundNucleus)/U;
223 G4double Db = 0.0;
224
225
226 G4double InvLevelDensity = CalcInvLevelDensity(static_cast<G4int>(A));
227
228
229 // bracketing the solution
230 if (Da == 0.0) {
231 TConf = Ta;
232 return 2*Ta*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Ta)*A13*A13;
233 } else if (Da < 0.0) {
234 do {
235 Tb -= 0.5*Tb;
236 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
237 Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
238 } while (Db < 0.0);
239 } else {
240 do {
241 Tb += 0.5*Tb;
242 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
243 Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
244 } while (Db > 0.0);
245 }
246
247
248 G4double eps = 1.0e-14 * std::abs(Tb-Ta);
249
250 for (G4int i = 0; i < 1000; i++) {
251 G4double Tc = (Ta+Tb)/2.0;
252 if (std::abs(Ta-Tb) <= eps) {
253 TConf = Tc;
254 return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
255 }
256 ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tc);
257 G4double Dc = (U+__FreeInternalE0-ECompoundNucleus)/U;
258
259 if (Dc == 0.0) {
260 TConf = Tc;
261 return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
262 }
263
264 if (Da*Dc < 0.0) {
265 Tb = Tc;
266 Db = Dc;
267 } else {
268 Ta = Tc;
269 Da = Dc;
270 }
271 }
272
273 G4cerr <<
274 "G4StatMFMicrocanoncal::CalcEntropyOfCompoundNucleus: I can't calculate the temperature"
275 << G4endl;
276
277 return 0.0;
278}
279
280
281
282
283G4StatMFChannel * G4StatMFMicroCanonical::ChooseAandZ(const G4Fragment & theFragment)
284 // Choice of fragment atomic numbers and charges
285{
286 // We choose a multiplicity (1,2,3,...) and then a channel
287 G4double RandNumber = G4UniformRand();
288
289 if (RandNumber < _WCompoundNucleus) {
290
291 G4StatMFChannel * aChannel = new G4StatMFChannel;
292 aChannel->CreateFragment(theFragment.GetA(),theFragment.GetZ());
293 return aChannel;
294
295 } else {
296
297 G4double AccumWeight = _WCompoundNucleus;
298 std::vector<G4StatMFMicroManager*>::iterator it;
299 for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it) {
300 AccumWeight += (*it)->GetProbability();
301 if (RandNumber < AccumWeight) {
302 return (*it)->ChooseChannel(theFragment.GetA(),theFragment.GetZ(),__MeanTemperature);
303 }
304 }
305 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::ChooseAandZ: wrong normalization!");
306 }
307
308 return 0;
309}
310
311
312
313
314G4double G4StatMFMicroCanonical::CalcInvLevelDensity(const G4int anA)
315{
316 // Calculate Inverse Density Level
317 // Epsilon0*(1 + 3 /(Af - 1))
318 if (anA == 1) return 0.0;
319 else return
320 G4StatMFParameters::GetEpsilon0()*(1.0+3.0/(anA - 1.0));
321}
Note: See TracBrowser for help on using the repository browser.