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