source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFChannel.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: 15.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: G4StatMFChannel.cc,v 1.10 2008/11/19 14:33:31 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// Hadronic Process: Nuclear De-excitations
31// by V. Lara
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 semi-infinite loop
37
38
39#include "G4StatMFChannel.hh"
40#include "G4HadronicException.hh"
41#include <numeric>
42
43class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double>
44{
45public:
46 SumCoulombEnergy() : total(0.0) {}
47 G4double operator() (G4double& , G4StatMFFragment*& frag)
48 {
49 total += frag->GetCoulombEnergy();
50 return total;
51 }
52
53 G4double GetTotal() { return total; }
54public:
55 G4double total;
56};
57
58
59
60
61
62// Copy constructor
63G4StatMFChannel::G4StatMFChannel(const G4StatMFChannel & )
64{
65 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::copy_constructor meant to not be accessable");
66}
67
68// Operators
69
70G4StatMFChannel & G4StatMFChannel::
71operator=(const G4StatMFChannel & )
72{
73 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator= meant to not be accessable");
74 return *this;
75}
76
77
78G4bool G4StatMFChannel::operator==(const G4StatMFChannel & ) const
79{
80 // throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator== meant to not be accessable");
81 return false;
82}
83
84
85G4bool G4StatMFChannel::operator!=(const G4StatMFChannel & ) const
86{
87 // throw G4HadronicException(__FILE__, __LINE__, "G4StatMFChannel::operator!= meant to not be accessable");
88 return true;
89}
90
91
92G4bool G4StatMFChannel::CheckFragments(void)
93{
94 std::deque<G4StatMFFragment*>::iterator i;
95 for (i = _theFragments.begin();
96 i != _theFragments.end(); ++i)
97 {
98 G4int A = static_cast<G4int>((*i)->GetA());
99 G4int Z = static_cast<G4int>((*i)->GetZ());
100 if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
101 }
102
103 return true;
104}
105
106
107
108
109void G4StatMFChannel::CreateFragment(const G4double A, const G4double Z)
110 // Create a new fragment.
111 // Fragments are automatically sorted: first charged fragments,
112 // then neutral ones.
113{
114 if (Z <= 0.5) {
115 _theFragments.push_back(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
116 _NumOfNeutralFragments++;
117 } else {
118 _theFragments.push_front(new G4StatMFFragment(static_cast<G4int>(A),static_cast<G4int>(Z)));
119 _NumOfChargedFragments++;
120 }
121
122 return;
123}
124
125
126G4double G4StatMFChannel::GetFragmentsCoulombEnergy(void)
127{
128 G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
129 0.0,SumCoulombEnergy());
130// G4double Coulomb = 0.0;
131// for (unsigned int i = 0;i < _theFragments.size(); i++)
132// Coulomb += _theFragments[i]->GetCoulombEnergy();
133 return Coulomb;
134}
135
136
137
138G4double G4StatMFChannel::GetFragmentsEnergy(const G4double T) const
139{
140 G4double Energy = 0.0;
141
142 G4double TranslationalEnergy = (3./2.)*T*static_cast<G4double>(_theFragments.size());
143
144 std::deque<G4StatMFFragment*>::const_iterator i;
145 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
146 {
147 Energy += (*i)->GetEnergy(T);
148 }
149 return Energy + TranslationalEnergy;
150}
151
152G4FragmentVector * G4StatMFChannel::GetFragments(const G4double anA,
153 const G4double anZ,
154 const G4double T)
155 //
156{
157 // calculate momenta of charged fragments
158 CoulombImpulse(anA,anZ,T);
159
160 // calculate momenta of neutral fragments
161 FragmentsMomenta(_NumOfNeutralFragments, _NumOfChargedFragments, T);
162
163
164 G4FragmentVector * theResult = new G4FragmentVector;
165 std::deque<G4StatMFFragment*>::iterator i;
166 for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
167 theResult->push_back((*i)->GetFragment(T));
168
169 return theResult;
170
171}
172
173
174
175void G4StatMFChannel::CoulombImpulse(const G4double anA, const G4double anZ, const G4double T)
176 // Aafter breakup, fragments fly away under Coulomb field.
177 // This method calculates asymptotic fragments momenta.
178{
179 // First, we have to place the fragments inside of the original nucleus volume
180 PlaceFragments(anA);
181
182 // Second, we sample initial charged fragments momenta. There are
183 // _NumOfChargedFragments charged fragments and they start at the begining
184 // of the vector _theFragments (i.e. 0)
185 FragmentsMomenta(_NumOfChargedFragments, 0, T);
186
187 // Third, we have to figure out the asymptotic momenta of charged fragments
188 // For taht we have to solve equations of motion for fragments
189 SolveEqOfMotion(anA,anZ,T);
190
191 return;
192}
193
194
195
196void G4StatMFChannel::PlaceFragments(const G4double anA)
197 // This gives the position of fragments at the breakup instant.
198 // Fragments positions are sampled inside prolongated ellipsoid.
199{
200 const G4double R0 = G4StatMFParameters::Getr0();
201 const G4double Rsys = 2.0*R0*std::pow(anA,1./3.);
202
203 G4bool TooMuchIterations;
204 do
205 {
206 TooMuchIterations = false;
207
208 // Sample the position of the first fragment
209 G4double R = (Rsys - R0*std::pow(_theFragments[0]->GetA(),1./3.))*
210 std::pow(G4UniformRand(),1./3.);
211 _theFragments[0]->SetPosition(IsotropicVector(R));
212
213
214 // Sample the position of the remaining fragments
215 G4bool ThereAreOverlaps = false;
216 std::deque<G4StatMFFragment*>::iterator i;
217 for (i = _theFragments.begin()+1; i != _theFragments.end(); ++i)
218 {
219 G4int counter = 0;
220 do
221 {
222 R = (Rsys - R0*std::pow((*i)->GetA(),1./3.))*std::pow(G4UniformRand(),1./3.);
223 (*i)->SetPosition(IsotropicVector(R));
224
225 // Check that there are not overlapping fragments
226 std::deque<G4StatMFFragment*>::iterator j;
227 for (j = _theFragments.begin(); j != i; ++j)
228 {
229 G4ThreeVector FragToFragVector = (*i)->GetPosition() - (*j)->GetPosition();
230 G4double Rmin = R0*(std::pow((*i)->GetA(),1./3.) +
231 std::pow((*j)->GetA(),1./3));
232 if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)) ) break;
233 }
234 counter++;
235 } while (ThereAreOverlaps && counter < 1000);
236
237 if (counter >= 1000)
238 {
239 TooMuchIterations = true;
240 break;
241 }
242 }
243 } while (TooMuchIterations);
244
245 return;
246}
247
248
249void G4StatMFChannel::FragmentsMomenta(const G4int NF, const G4int idx,
250 const G4double T)
251 // Calculate fragments momenta at the breakup instant
252 // Fragment kinetic energies are calculated according to the
253 // Boltzmann distribution at given temperature.
254 // NF is number of fragments
255 // idx is index of first fragment
256{
257 G4double KinE = (3./2.)*T*static_cast<G4double>(NF);
258
259 G4ThreeVector p;
260
261 if (NF <= 0) return;
262 else if (NF == 1)
263 {
264 // We have only one fragment to deal with
265 p = IsotropicVector(std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE));
266 _theFragments[idx]->SetMomentum(p);
267 }
268 else if (NF == 2)
269 {
270 // We have only two fragment to deal with
271 G4double M1 = _theFragments[idx]->GetNuclearMass();
272 G4double M2 = _theFragments[idx+1]->GetNuclearMass();
273 p = IsotropicVector(std::sqrt(2.0*KinE*(M1*M2)/(M1+M2)));
274 _theFragments[idx]->SetMomentum(p);
275 _theFragments[idx+1]->SetMomentum(-p);
276 }
277 else
278 {
279 // We have more than two fragments
280 G4double AvailableE;
281 G4int i1,i2;
282 G4double SummedE;
283 G4ThreeVector SummedP;
284 do
285 {
286 // Fisrt sample momenta of NF-2 fragments
287 // according to Boltzmann distribution
288 AvailableE = 0.0;
289 SummedE = 0.0;
290 SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
291 for (G4int i = idx; i < idx+NF-2; i++)
292 {
293 G4double E;
294 G4double RandE;
295 G4double Boltzmann;
296 do
297 {
298 E = 9.0*T*G4UniformRand();
299 Boltzmann = std::sqrt(E)*std::exp(-E/T);
300 RandE = std::sqrt(T/2.)*std::exp(-0.5)*G4UniformRand();
301 }
302 while (RandE > Boltzmann);
303 p = IsotropicVector(std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass()));
304 _theFragments[i]->SetMomentum(p);
305 SummedE += E;
306 SummedP += p;
307 }
308 // Calculate momenta of last two fragments in such a way
309 // that constraints are satisfied
310 i1 = idx+NF-2; // before last fragment index
311 i2 = idx+NF-1; // last fragment index
312 p = -SummedP;
313 AvailableE = KinE - SummedE;
314 // Available Kinetic Energy should be shared between two last fragments
315 }
316 while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
317 _theFragments[i2]->GetNuclearMass())));
318
319 G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()/_theFragments[i1]->GetNuclearMass();
320 G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()*AvailableE/p.mag2());
321 G4double CosTheta1;
322 G4double Sign;
323
324 if (CTM12 > 0.9999) {CosTheta1 = 1.;}
325 else {
326 do
327 {
328 do
329 {
330 CosTheta1 = 1.0 - 2.0*G4UniformRand();
331 }
332 while (CosTheta1*CosTheta1 < CTM12);
333 }
334 while (CTM12 >= 0.0 && CosTheta1 < 0.0);
335 }
336
337 if (CTM12 < 0.0) Sign = 1.0;
338 else if (G4UniformRand() <= 0.5) Sign = -1.0;
339 else Sign = 1.0;
340
341
342 G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()*(CosTheta1*CosTheta1-CTM12)))/H;
343 G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
344 G4double Phi = twopi*G4UniformRand();
345 G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
346 G4double CosPhi1 = std::cos(Phi);
347 G4double SinPhi1 = std::sin(Phi);
348 G4double CosPhi2 = -CosPhi1;
349 G4double SinPhi2 = -SinPhi1;
350 G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
351 G4double SinTheta2 = 0.0;
352 if (CosTheta2 > -1.0 && CosTheta2 < 1.0) SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
353
354 G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
355 G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
356 G4ThreeVector b(1.0,0.0,0.0);
357
358 p1 = RotateMomentum(p,b,p1);
359 p2 = RotateMomentum(p,b,p2);
360
361 SummedP += p1 + p2;
362 SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) +
363 p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());
364
365 _theFragments[i1]->SetMomentum(p1);
366 _theFragments[i2]->SetMomentum(p2);
367
368 }
369
370 return;
371}
372
373
374void G4StatMFChannel::SolveEqOfMotion(const G4double anA, const G4double anZ, const G4double T)
375 // This method will find a solution of Newton's equation of motion
376 // for fragments in the self-consistent time-dependent Coulomb field
377{
378 G4double CoulombEnergy = (3./5.)*(elm_coupling*anZ*anZ)*
379 std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
380 (G4StatMFParameters::Getr0()*std::pow(anA,1./3.))
381 - GetFragmentsCoulombEnergy();
382 if (CoulombEnergy <= 0.0) return;
383
384 G4int Iterations = 0;
385 G4double TimeN = 0.0;
386 G4double TimeS = 0.0;
387 G4double DeltaTime = 10.0;
388
389 G4ThreeVector * Pos = new G4ThreeVector[_NumOfChargedFragments];
390 G4ThreeVector * Vel = new G4ThreeVector[_NumOfChargedFragments];
391 G4ThreeVector * Accel = new G4ThreeVector[_NumOfChargedFragments];
392
393 G4int i;
394 for (i = 0; i < _NumOfChargedFragments; i++)
395 {
396 Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
397 _theFragments[i]->GetMomentum();
398 Pos[i] = _theFragments[i]->GetPosition();
399 }
400
401 do
402 {
403
404 G4ThreeVector distance;
405 G4ThreeVector force;
406
407 for (i = 0; i < _NumOfChargedFragments; i++)
408 {
409 force.setX(0.0); force.setY(0.0); force.setZ(0.0);
410 for (G4int j = 0; j < _NumOfChargedFragments; j++)
411 {
412 if (i != j)
413 {
414 distance = Pos[i] - Pos[j];
415 force += (elm_coupling*(_theFragments[i]->GetZ()*_theFragments[j]->GetZ())/
416 (distance.mag2()*distance.mag()))*distance;
417 }
418 }
419 Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
420 }
421
422 TimeN = TimeS + DeltaTime;
423
424 G4ThreeVector SavedVel;
425 for ( i = 0; i < _NumOfChargedFragments; i++)
426 {
427 SavedVel = Vel[i];
428 Vel[i] += Accel[i]*(TimeN-TimeS);
429 Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
430 }
431
432 // if (Iterations >= 50 && Iterations < 75) DeltaTime = 4.;
433 // else if (Iterations >= 75) DeltaTime = 10.;
434
435 TimeS = TimeN;
436
437 }
438 while (Iterations++ < 100);
439
440 // Summed fragment kinetic energy
441 G4double TotalKineticEnergy = 0.0;
442 for (i = 0; i < _NumOfChargedFragments; i++)
443 {
444 TotalKineticEnergy += _theFragments[i]->GetNuclearMass()*
445 0.5*Vel[i].mag2();
446 }
447 // Scaling of fragment velocities
448 G4double KineticEnergy = (3./2.)*static_cast<G4double>(_theFragments.size())*T;
449 G4double Eta = ( CoulombEnergy + KineticEnergy ) / TotalKineticEnergy;
450 for (i = 0; i < _NumOfChargedFragments; i++)
451 {
452 Vel[i] *= Eta;
453 }
454
455 // Finally calculate fragments momenta
456 for (i = 0; i < _NumOfChargedFragments; i++)
457 {
458 _theFragments[i]->SetMomentum(_theFragments[i]->GetNuclearMass()*Vel[i]);
459 }
460
461 // garbage collection
462 delete [] Pos;
463 delete [] Vel;
464 delete [] Accel;
465
466 return;
467}
468
469
470
471G4ThreeVector G4StatMFChannel::RotateMomentum(G4ThreeVector Pa,
472 G4ThreeVector V, G4ThreeVector P)
473 // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0
474{
475 G4ThreeVector U = Pa.unit();
476
477 G4double Alpha1 = U * V;
478
479 G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
480
481 G4ThreeVector N = (1./Alpha2)*U.cross(V);
482
483 G4ThreeVector RotatedMomentum(
484 ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
485 ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
486 ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
487 );
488 return RotatedMomentum;
489}
490
491
492
493
494
495G4ThreeVector G4StatMFChannel::IsotropicVector(const G4double Magnitude)
496 // Samples a isotropic random vector with a magnitud given by Magnitude.
497 // By default Magnitude = 1
498{
499 G4double CosTheta = 1.0 - 2.0*G4UniformRand();
500 G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
501 G4double Phi = twopi*G4UniformRand();
502 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
503 Magnitude*std::cos(Phi)*CosTheta,
504 Magnitude*std::sin(Phi));
505 return Vector;
506}
Note: See TracBrowser for help on using the repository browser.