source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMF.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: 10.3 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: G4StatMF.cc,v 1.6 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#include "G4StatMF.hh"
34
35
36
37// Default constructor
38G4StatMF::G4StatMF() : _theEnsemble(0) {}
39
40
41// Destructor
42G4StatMF::~G4StatMF() {} //{if (_theEnsemble != 0) delete _theEnsemble;}
43
44
45// Copy constructor
46G4StatMF::G4StatMF(const G4StatMF & ) : G4VMultiFragmentation()
47{
48 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::copy_constructor meant to not be accessable");
49}
50
51
52// Operators
53
54G4StatMF & G4StatMF::operator=(const G4StatMF & )
55{
56 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator= meant to not be accessable");
57 return *this;
58}
59
60
61G4bool G4StatMF::operator==(const G4StatMF & )
62{
63 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator== meant to not be accessable");
64 return false;
65}
66
67
68G4bool G4StatMF::operator!=(const G4StatMF & )
69{
70 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator!= meant to not be accessable");
71 return true;
72}
73
74
75
76
77
78G4FragmentVector * G4StatMF::BreakItUp(const G4Fragment &theFragment)
79{
80 // G4FragmentVector * theResult = new G4FragmentVector;
81
82 if (theFragment.GetExcitationEnergy() <= 0.0) {
83 G4FragmentVector * theResult = new G4FragmentVector;
84 theResult->push_back(new G4Fragment(theFragment));
85 return 0;
86 }
87
88
89 // Maximun average multiplicity: M_0 = 2.6 for A ~ 200
90 // and M_0 = 3.3 for A <= 110
91 G4double MaxAverageMultiplicity =
92 G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA()));
93
94
95 // We'll use two kinds of ensembles
96 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
97 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
98
99
100 //-------------------------------------------------------
101 // Direct simulation part (Microcanonical ensemble)
102 //-------------------------------------------------------
103
104 // Microcanonical ensemble initialization
105 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
106
107 G4int Iterations = 0;
108 G4int IterationsLimit = 100000;
109 G4double Temperature = 0.0;
110
111 G4bool FirstTime = true;
112 G4StatMFChannel * theChannel = 0;
113
114 G4bool ChannelOk;
115 do { // Try to de-excite as much as IterationLimit permits
116 do {
117
118 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
119 if (theMeanMult <= MaxAverageMultiplicity) {
120 // G4cout << "MICROCANONICAL" << G4endl;
121 // Choose fragments atomic numbers and charges from direct simulation
122 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
123 _theEnsemble = theMicrocanonicalEnsemble;
124 } else {
125 //-----------------------------------------------------
126 // Non direct simulation part (Macrocanonical Ensemble)
127 //-----------------------------------------------------
128 if (FirstTime) {
129 // Macrocanonical ensemble initialization
130 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
131 _theEnsemble = theMacrocanonicalEnsemble;
132 FirstTime = false;
133 }
134 // G4cout << "MACROCANONICAL" << G4endl;
135 // Select calculated fragment total multiplicity,
136 // fragment atomic numbers and fragment charges.
137 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
138 }
139
140 ChannelOk = theChannel->CheckFragments();
141 if (!ChannelOk) delete theChannel;
142
143 } while (!ChannelOk);
144
145
146 if (theChannel->GetMultiplicity() <= 1) {
147 G4FragmentVector * theResult = new G4FragmentVector;
148 theResult->push_back(new G4Fragment(theFragment));
149 delete theMicrocanonicalEnsemble;
150 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
151 delete theChannel;
152 return theResult;
153 }
154
155 //--------------------------------------
156 // Second part of simulation procedure.
157 //--------------------------------------
158
159 // Find temperature of breaking channel.
160 Temperature = _theEnsemble->GetMeanTemperature(); // Initial guess for Temperature
161
162 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
163
164 // Do not forget to delete this unusable channel, for which we failed to find the temperature,
165 // otherwise for very proton-reach nuclei it would lead to memory leak due to large
166 // number of iterations. N.B. "theChannel" is created in G4StatMFMacroCanonical::ChooseZ()
167
168 // G4cout << " Iteration # " << Iterations << " Mean Temperature = " << Temperature << G4endl;
169
170 delete theChannel;
171
172 } while (Iterations++ < IterationsLimit );
173
174
175
176 // If Iterations >= IterationsLimit means that we couldn't solve for temperature
177 if (Iterations >= IterationsLimit)
178 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
179
180
181 G4FragmentVector * theResult = theChannel->
182 GetFragments(theFragment.GetA(),theFragment.GetZ(),Temperature);
183
184
185
186 // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
187 // Original nucleus 4-momentum in CM system
188 G4LorentzVector InitialMomentum(theFragment.GetMomentum());
189 InitialMomentum.boost(-InitialMomentum.boostVector());
190 G4double ScaleFactor = 0.0;
191 G4double SavedScaleFactor = 0.0;
192 do {
193 G4double FragmentsEnergy = 0.0;
194 G4FragmentVector::iterator j;
195 for (j = theResult->begin(); j != theResult->end(); j++)
196 FragmentsEnergy += (*j)->GetMomentum().e();
197 SavedScaleFactor = ScaleFactor;
198 ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
199 G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
200 for (j = theResult->begin(); j != theResult->end(); j++) {
201 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
202 G4double Mass = (*j)->GetMomentum().m();
203 G4LorentzVector NewMomentum;
204 NewMomentum.setVect(ScaledMomentum);
205 NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
206 (*j)->SetMomentum(NewMomentum);
207 }
208 } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
209 // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
210
211 // Perform Lorentz boost
212 G4FragmentVector::iterator i;
213 for (i = theResult->begin(); i != theResult->end(); i++) {
214 G4LorentzVector FourMom = (*i)->GetMomentum();
215 FourMom.boost(theFragment.GetMomentum().boostVector());
216 (*i)->SetMomentum(FourMom);
217#ifdef PRECOMPOUND_TEST
218 (*i)->SetCreatorModel(G4String("G4StatMF"));
219#endif
220 }
221
222 // garbage collection
223 delete theMicrocanonicalEnsemble;
224 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
225 delete theChannel;
226
227 return theResult;
228}
229
230
231G4bool G4StatMF::FindTemperatureOfBreakingChannel(const G4Fragment & theFragment,
232 const G4StatMFChannel * aChannel,
233 G4double & Temperature)
234 // This finds temperature of breaking channel.
235{
236 G4double A = theFragment.GetA();
237 G4double Z = theFragment.GetZ();
238 G4double U = theFragment.GetExcitationEnergy();
239
240 G4double T = std::max(Temperature,0.0012*MeV);
241
242 G4double Ta = T;
243 G4double Tb = T;
244
245
246 G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
247
248 G4double Da = (U - TotalEnergy)/U;
249 G4double Db = 0.0;
250
251 // bracketing the solution
252 if (Da == 0.0) {
253 Temperature = T;
254 return true;
255 } else if (Da < 0.0) {
256 do {
257 Tb -= 0.5 * std::abs(Tb);
258 T = Tb;
259 if (Tb < 0.001*MeV) return false;
260
261 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
262
263 Db = (U - TotalEnergy)/U;
264 } while (Db < 0.0);
265
266 } else {
267 do {
268 Tb += 0.5 * std::abs(Tb);
269 T = Tb;
270
271 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
272
273 Db = (U - TotalEnergy)/U;
274 } while (Db > 0.0);
275 }
276
277 G4double eps = 1.0e-14 * std::abs(Tb-Ta);
278 //G4double eps = 1.0e-3 ;
279
280 // Start the bisection method
281 for (G4int j = 0; j < 1000; j++) {
282 G4double Tc = (Ta+Tb)/2.0;
283 if (std::abs(Ta-Tc) <= eps) {
284 Temperature = Tc;
285 return true;
286 }
287
288 T = Tc;
289
290 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
291
292 G4double Dc = (U - TotalEnergy)/U;
293
294 if (Dc == 0.0) {
295 Temperature = Tc;
296 return true;
297 }
298
299 if (Da*Dc < 0.0) {
300 Tb = Tc;
301 Db = Dc;
302 } else {
303 Ta = Tc;
304 Da = Dc;
305 }
306 }
307
308 Temperature = (Ta+Tb)/2.0;
309 return false;
310}
311
312
313
314G4double G4StatMF::CalcEnergy(const G4double A, const G4double Z, const G4StatMFChannel * aChannel,
315 const G4double T)
316{
317 G4double MassExcess0 = G4NucleiProperties::GetMassExcess(static_cast<G4int>(A),static_cast<G4int>(Z));
318
319 G4double Coulomb = (3./5.)*(elm_coupling*Z*Z)*std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
320 (G4StatMFParameters::Getr0()*std::pow(A,1./3.));
321
322 G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T);
323
324 return -MassExcess0 + Coulomb + ChannelEnergy;
325
326}
327
328
329
Note: See TracBrowser for help on using the repository browser.