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