source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroCanonical.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.1 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: G4StatMFMacroCanonical.cc,v 1.8 2008/11/19 14:33:31 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// by V. Lara
31// --------------------------------------------------------------------
32//
33// Modified:
34// 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
35// Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
36// Moscow, pshenich@fias.uni-frankfurt.de) fixed infinite loop for
37// a fagment with Z=A; fixed memory leak
38
39#include "G4StatMFMacroCanonical.hh"
40
41
42// constructor
43G4StatMFMacroCanonical::G4StatMFMacroCanonical(const G4Fragment & theFragment)
44{
45
46 // Get memory for clusters
47 _theClusters.push_back(new G4StatMFMacroNucleon); // Size 1
48 _theClusters.push_back(new G4StatMFMacroBiNucleon); // Size 2
49 _theClusters.push_back(new G4StatMFMacroTriNucleon); // Size 3
50 _theClusters.push_back(new G4StatMFMacroTetraNucleon); // Size 4
51 for (G4int i = 4; i < theFragment.GetA(); i++)
52 _theClusters.push_back(new G4StatMFMacroMultiNucleon(i+1)); // Size 5 ... A
53
54 // Perform class initialization
55 Initialize(theFragment);
56
57}
58
59
60// destructor
61G4StatMFMacroCanonical::~G4StatMFMacroCanonical()
62{
63 // garbage collection
64 if (!_theClusters.empty())
65 {
66 std::for_each(_theClusters.begin(),_theClusters.end(),DeleteFragment());
67 }
68}
69
70// operators definitions
71G4StatMFMacroCanonical &
72G4StatMFMacroCanonical::operator=(const G4StatMFMacroCanonical & )
73{
74 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroCanonical::operator= meant to not be accessable");
75 return *this;
76}
77
78G4bool G4StatMFMacroCanonical::operator==(const G4StatMFMacroCanonical & ) const
79{
80 return false;
81}
82
83
84G4bool G4StatMFMacroCanonical::operator!=(const G4StatMFMacroCanonical & ) const
85{
86 return true;
87}
88
89
90// Initialization method
91
92
93void G4StatMFMacroCanonical::Initialize(const G4Fragment & theFragment)
94{
95
96 G4double A = theFragment.GetA();
97 G4double Z = theFragment.GetZ();
98
99 // Free Internal energy at T = 0
100 __FreeInternalE0 = A*( -G4StatMFParameters::GetE0() + // Volume term (for T = 0)
101 G4StatMFParameters::GetGamma0()* // Symmetry term
102 (1.0-2.0*Z/A)*(1.0-2.0*Z/A) ) +
103 G4StatMFParameters::GetBeta0()*std::pow(A,2.0/3.0) + // Surface term (for T = 0)
104 (3.0/5.0)*elm_coupling*Z*Z/(G4StatMFParameters::Getr0()* // Coulomb term
105 std::pow(A,1.0/3.0));
106
107
108
109 CalculateTemperature(theFragment);
110
111 return;
112}
113
114
115
116
117
118void G4StatMFMacroCanonical::CalculateTemperature(const G4Fragment & theFragment)
119{
120 // Excitation Energy
121 G4double U = theFragment.GetExcitationEnergy();
122
123 G4double A = theFragment.GetA();
124 G4double Z = theFragment.GetZ();
125
126 // Fragment Multiplicity
127 G4double FragMult = std::max((1.0+(2.31/MeV)*(U/A - 3.5*MeV))*A/100.0, 2.0);
128
129
130 // Parameter Kappa
131 _Kappa = (1.0+elm_coupling*(std::pow(FragMult,1./3.)-1)/
132 (G4StatMFParameters::Getr0()*std::pow(A,1./3.)));
133 _Kappa = _Kappa*_Kappa*_Kappa - 1.0;
134
135
136 G4StatMFMacroTemperature * theTemp = new
137 G4StatMFMacroTemperature(A,Z,U,__FreeInternalE0,_Kappa,&_theClusters);
138
139 __MeanTemperature = theTemp->CalcTemperature();
140 _ChemPotentialNu = theTemp->GetChemicalPotentialNu();
141 _ChemPotentialMu = theTemp->GetChemicalPotentialMu();
142 __MeanMultiplicity = theTemp->GetMeanMultiplicity();
143 __MeanEntropy = theTemp->GetEntropy();
144
145 delete theTemp;
146
147 return;
148}
149
150
151// --------------------------------------------------------------------------
152
153G4StatMFChannel * G4StatMFMacroCanonical::ChooseAandZ(const G4Fragment &theFragment)
154 // Calculate total fragments multiplicity, fragment atomic numbers and charges
155{
156 G4double A = theFragment.GetA();
157 G4double Z = theFragment.GetZ();
158
159 std::vector<G4double> ANumbers(static_cast<G4int>(A));
160
161 G4double Multiplicity = ChooseA(A,ANumbers);
162
163
164 std::vector<G4double> FragmentsA;
165
166 G4int i = 0;
167 for (i = 0; i < A; i++)
168 {
169 for (G4int j = 0; j < ANumbers[i]; j++) FragmentsA.push_back(i+1);
170 }
171
172
173 // Sort fragments in decreasing order
174 G4int im = 0;
175 for (G4int j = 0; j < Multiplicity; j++)
176 {
177 G4double FragmentsAMax = 0.0;
178 im = j;
179 for (i = j; i < Multiplicity; i++)
180 {
181 if (FragmentsA[i] <= FragmentsAMax) continue;
182 else
183 {
184 im = i;
185 FragmentsAMax = FragmentsA[im];
186 }
187 }
188
189 if (im != j)
190 {
191 FragmentsA[im] = FragmentsA[j];
192 FragmentsA[j] = FragmentsAMax;
193 }
194 }
195
196 return ChooseZ(static_cast<G4int>(Z),FragmentsA);
197}
198
199
200
201G4double G4StatMFMacroCanonical::ChooseA(const G4double A, std::vector<G4double> & ANumbers)
202 // Determines fragments multiplicities and compute total fragment multiplicity
203{
204 G4double multiplicity = 0.0;
205 G4int i;
206
207
208 std::vector<G4double> AcumMultiplicity;
209 AcumMultiplicity.reserve(static_cast<G4int>(A));
210
211 AcumMultiplicity.push_back((*(_theClusters.begin()))->GetMeanMultiplicity());
212 for (std::vector<G4VStatMFMacroCluster*>::iterator it = _theClusters.begin()+1;
213 it != _theClusters.end(); ++it)
214 {
215 AcumMultiplicity.push_back((*it)->GetMeanMultiplicity()+AcumMultiplicity.back());
216 }
217
218 G4int CheckA;
219 do {
220 CheckA = -1;
221 G4int SumA = 0;
222 G4int ThisOne = 0;
223 multiplicity = 0.0;
224 for (i = 0; i < A; i++) ANumbers[i] = 0.0;
225 do {
226 G4double RandNumber = G4UniformRand()*__MeanMultiplicity;
227 for (i = 0; i < A; i++) {
228 if (RandNumber < AcumMultiplicity[i]) {
229 ThisOne = i;
230 break;
231 }
232 }
233 multiplicity++;
234 ANumbers[ThisOne] = ANumbers[ThisOne]+1;
235 SumA += ThisOne+1;
236 CheckA = static_cast<G4int>(A) - SumA;
237
238 } while (CheckA > 0);
239
240 } while (CheckA < 0 || std::abs(__MeanMultiplicity - multiplicity) > std::sqrt(__MeanMultiplicity) + 1./2.);
241
242 return multiplicity;
243}
244
245
246G4StatMFChannel * G4StatMFMacroCanonical::ChooseZ(const G4int & Z,
247 std::vector<G4double> & FragmentsA)
248 //
249{
250 std::vector<G4double> FragmentsZ;
251
252 G4double DeltaZ = 0.0;
253 G4double CP = (3./5.)*(elm_coupling/G4StatMFParameters::Getr0())*
254 (1.0 - 1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.));
255
256 G4int multiplicity = FragmentsA.size();
257
258 do
259 {
260 FragmentsZ.clear();
261 G4int SumZ = 0;
262 for (G4int i = 0; i < multiplicity; i++)
263 {
264 G4double A = FragmentsA[i];
265 if (A <= 1.0)
266 {
267 G4double RandNumber = G4UniformRand();
268 if (RandNumber < (*_theClusters.begin())->GetZARatio())
269 {
270 FragmentsZ.push_back(1.0);
271 SumZ += static_cast<G4int>(FragmentsZ[i]);
272 }
273 else FragmentsZ.push_back(0.0);
274 }
275 else
276 {
277 G4double RandZ;
278 G4double CC = 8.0*G4StatMFParameters::GetGamma0()+2.0*CP*std::pow(FragmentsA[i],2./3.);
279 G4double ZMean;
280 if (FragmentsA[i] > 1.5 && FragmentsA[i] < 4.5) ZMean = 0.5*FragmentsA[i];
281 else ZMean = FragmentsA[i]*(4.0*G4StatMFParameters::GetGamma0()+_ChemPotentialNu)/CC;
282 G4double ZDispersion = std::sqrt(FragmentsA[i]*__MeanTemperature/CC);
283 G4int z;
284 do
285 {
286 RandZ = G4RandGauss::shoot(ZMean,ZDispersion);
287 z = static_cast<G4int>(RandZ+0.5);
288 } while (z < 0 || z > A);
289 FragmentsZ.push_back(z);
290 SumZ += z;
291 }
292 }
293 DeltaZ = Z - SumZ;
294 }
295 while (std::abs(DeltaZ) > 1.1);
296
297 // DeltaZ can be 0, 1 or -1
298 G4int idx = 0;
299 if (DeltaZ < 0.0)
300 {
301 while (FragmentsZ[idx] < 0.5) ++idx;
302 }
303 FragmentsZ[idx] += DeltaZ;
304
305 G4StatMFChannel * theChannel = new G4StatMFChannel;
306 for (G4int i = multiplicity-1; i >= 0; i--)
307 {
308 theChannel->CreateFragment(FragmentsA[i],FragmentsZ[i]);
309 }
310
311
312 return theChannel;
313}
Note: See TracBrowser for help on using the repository browser.