Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4StatMFChannel.cc
이 파일의 문서화 페이지로 가기
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 107060 2017-11-01 15:00:04Z 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 #include "G4RandomDirection.hh"
46 
47 class SumCoulombEnergy : public std::binary_function<G4double,G4double,G4double>
48 {
49 public:
50  SumCoulombEnergy() : total(0.0) {}
52  {
53  total += frag->GetCoulombEnergy();
54  return total;
55  }
56 
57  G4double GetTotal() { return total; }
58 public:
60 };
61 
63  _NumOfNeutralFragments(0),
64  _NumOfChargedFragments(0)
65 {}
66 
68 {
69  if (!_theFragments.empty()) {
70  std::for_each(_theFragments.begin(),_theFragments.end(),
71  DeleteFragment());
72  }
73 }
74 
76 {
77  std::deque<G4StatMFFragment*>::iterator i;
78  for (i = _theFragments.begin();
79  i != _theFragments.end(); ++i)
80  {
81  G4int A = (*i)->GetA();
82  G4int Z = (*i)->GetZ();
83  if ( (A > 1 && (Z > A || Z <= 0)) || (A==1 && Z > A) || A <= 0 ) return false;
84  }
85  return true;
86 }
87 
89 // Create a new fragment.
90 // Fragments are automatically sorted: first charged fragments,
91 // then neutral ones.
92 {
93  if (Z <= 0.5) {
94  _theFragments.push_back(new G4StatMFFragment(A,Z));
96  } else {
97  _theFragments.push_front(new G4StatMFFragment(A,Z));
99  }
100 
101  return;
102 }
103 
105 {
106  G4double Coulomb = std::accumulate(_theFragments.begin(),_theFragments.end(),
107  0.0,SumCoulombEnergy());
108  // G4double Coulomb = 0.0;
109  // for (unsigned int i = 0;i < _theFragments.size(); i++)
110  // Coulomb += _theFragments[i]->GetCoulombEnergy();
111  return Coulomb;
112 }
113 
115 {
116  G4double Energy = 0.0;
117 
118  G4double TranslationalEnergy = 1.5*T*_theFragments.size();
119 
120  std::deque<G4StatMFFragment*>::const_iterator i;
121  for (i = _theFragments.begin(); i != _theFragments.end(); ++i)
122  {
123  Energy += (*i)->GetEnergy(T);
124  }
125  return Energy + TranslationalEnergy;
126 }
127 
129  G4int anZ,
130  G4double T)
131 {
132  // calculate momenta of charged fragments
133  CoulombImpulse(anA,anZ,T);
134 
135  // calculate momenta of neutral fragments
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 
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)
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 
166 // This gives the position of fragments at the breakup instant.
167 // Fragments positions are sampled inside prolongated ellipsoid.
168 {
169  G4Pow* g4calc = G4Pow::GetInstance();
170  const G4double R0 = G4StatMFParameters::Getr0();
171  G4double Rsys = 2.0*R0*g4calc->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*g4calc->Z13(_theFragments[0]->GetA()))*
180  g4calc->A13(G4UniformRand());
181  _theFragments[0]->SetPosition(R*G4RandomDirection());
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*g4calc->Z13((*i)->GetA()))*g4calc->A13(G4UniformRand());
193  (*i)->SetPosition(R*G4RandomDirection());
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 =
200  (*i)->GetPosition() - (*j)->GetPosition();
201  G4double Rmin = R0*(g4calc->Z13((*i)->GetA()) +
202  g4calc->Z13((*j)->GetA()));
203  if ( (ThereAreOverlaps = (FragToFragVector.mag2() < Rmin*Rmin)))
204  { break; }
205  }
206  counter++;
207  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
208  } while (ThereAreOverlaps && counter < 1000);
209 
210  if (counter >= 1000)
211  {
212  TooMuchIterations = true;
213  break;
214  }
215  }
216  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
217  } while (TooMuchIterations);
218  return;
219 }
220 
222  G4double T)
223 // Calculate fragments momenta at the breakup instant
224 // Fragment kinetic energies are calculated according to the
225 // Boltzmann distribution at given temperature.
226 // NF is number of fragments
227 // idx is index of first fragment
228 {
229  G4double KinE = 1.5*T*NF;
230  G4ThreeVector p(0.,0.,0.);
231 
232  if (NF <= 0) return;
233  else if (NF == 1)
234  {
235  // We have only one fragment to deal with
236  p = std::sqrt(2.0*_theFragments[idx]->GetNuclearMass()*KinE)
238  _theFragments[idx]->SetMomentum(p);
239  }
240  else if (NF == 2)
241  {
242  // We have only two fragment to deal with
243  G4double M1 = _theFragments[idx]->GetNuclearMass();
244  G4double M2 = _theFragments[idx+1]->GetNuclearMass();
245  p = std::sqrt(2.0*KinE*(M1*M2)/(M1+M2))*G4RandomDirection();
246  _theFragments[idx]->SetMomentum(p);
247  _theFragments[idx+1]->SetMomentum(-p);
248  }
249  else
250  {
251  // We have more than two fragments
252  G4double AvailableE;
253  G4int i1,i2;
254  G4double SummedE;
255  G4ThreeVector SummedP(0.,0.,0.);
256  do
257  {
258  // Fisrt sample momenta of NF-2 fragments
259  // according to Boltzmann distribution
260  AvailableE = 0.0;
261  SummedE = 0.0;
262  SummedP.setX(0.0);SummedP.setY(0.0);SummedP.setZ(0.0);
263  for (G4int i = idx; i < idx+NF-2; ++i)
264  {
265  G4double E;
266  G4double RandE;
267  do
268  {
269  E = 9.0*G4UniformRand();
270  RandE = std::sqrt(0.5/E)*G4Exp(E-0.5)*G4UniformRand();
271  }
272  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
273  while (RandE > 1.0);
274  E *= T;
275  p = std::sqrt(2.0*E*_theFragments[i]->GetNuclearMass())
277  _theFragments[i]->SetMomentum(p);
278  SummedE += E;
279  SummedP += p;
280  }
281  // Calculate momenta of last two fragments in such a way
282  // that constraints are satisfied
283  i1 = idx+NF-2; // before last fragment index
284  i2 = idx+NF-1; // last fragment index
285  p = -SummedP;
286  AvailableE = KinE - SummedE;
287  // Available Kinetic Energy should be shared between two last fragments
288  }
289  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
290  while (AvailableE <= p.mag2()/(2.0*(_theFragments[i1]->GetNuclearMass()+
291  _theFragments[i2]->GetNuclearMass())));
292  G4double H = 1.0 + _theFragments[i2]->GetNuclearMass()
293  /_theFragments[i1]->GetNuclearMass();
294  G4double CTM12 = H*(1.0 - 2.0*_theFragments[i2]->GetNuclearMass()
295  *AvailableE/p.mag2());
296  G4double CosTheta1;
297  G4double Sign;
298 
299  if (CTM12 > 1.) {CosTheta1 = 1.;}
300  else {
301  do
302  {
303  do
304  {
305  CosTheta1 = 1.0 - 2.0*G4UniformRand();
306  }
307  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
308  while (CosTheta1*CosTheta1 < CTM12);
309  }
310  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
311  while (CTM12 >= 0.0 && CosTheta1 < 0.0);
312  }
313 
314  if (CTM12 < 0.0) Sign = 1.0;
315  else if (G4UniformRand() <= 0.5) Sign = -1.0;
316  else Sign = 1.0;
317 
318  G4double P1 = (p.mag()*CosTheta1+Sign*std::sqrt(p.mag2()
319  *(CosTheta1*CosTheta1-CTM12)))/H;
320  G4double P2 = std::sqrt(P1*P1+p.mag2() - 2.0*P1*p.mag()*CosTheta1);
321  G4double Phi = twopi*G4UniformRand();
322  G4double SinTheta1 = std::sqrt(1.0 - CosTheta1*CosTheta1);
323  G4double CosPhi1 = std::cos(Phi);
324  G4double SinPhi1 = std::sin(Phi);
325  G4double CosPhi2 = -CosPhi1;
326  G4double SinPhi2 = -SinPhi1;
327  G4double CosTheta2 = (p.mag2() + P2*P2 - P1*P1)/(2.0*p.mag()*P2);
328  G4double SinTheta2 = 0.0;
329  if (CosTheta2 > -1.0 && CosTheta2 < 1.0) {
330  SinTheta2 = std::sqrt(1.0 - CosTheta2*CosTheta2);
331  }
332  G4ThreeVector p1(P1*SinTheta1*CosPhi1,P1*SinTheta1*SinPhi1,P1*CosTheta1);
333  G4ThreeVector p2(P2*SinTheta2*CosPhi2,P2*SinTheta2*SinPhi2,P2*CosTheta2);
334  G4ThreeVector b(1.0,0.0,0.0);
335 
336  p1 = RotateMomentum(p,b,p1);
337  p2 = RotateMomentum(p,b,p2);
338 
339  SummedP += p1 + p2;
340  SummedE += p1.mag2()/(2.0*_theFragments[i1]->GetNuclearMass()) +
341  p2.mag2()/(2.0*_theFragments[i2]->GetNuclearMass());
342 
343  _theFragments[i1]->SetMomentum(p1);
344  _theFragments[i2]->SetMomentum(p2);
345 
346  }
347  return;
348 }
349 
351 // This method will find a solution of Newton's equation of motion
352 // for fragments in the self-consistent time-dependent Coulomb field
353 {
354  G4Pow* g4calc = G4Pow::GetInstance();
355  G4double CoulombEnergy = 0.6*elm_coupling*anZ*anZ*
358  if (CoulombEnergy <= 0.0) return;
359 
360  G4int Iterations = 0;
361  G4double TimeN = 0.0;
362  G4double TimeS = 0.0;
363  G4double DeltaTime = 10.0;
364 
368 
369  G4int i;
370  for (i = 0; i < _NumOfChargedFragments; i++)
371  {
372  Vel[i] = (1.0/(_theFragments[i]->GetNuclearMass()))*
373  _theFragments[i]->GetMomentum();
374  Pos[i] = _theFragments[i]->GetPosition();
375  }
376 
377  G4ThreeVector distance(0.,0.,0.);
378  G4ThreeVector force(0.,0.,0.);
379  G4ThreeVector SavedVel(0.,0.,0.);
380  do {
381  for (i = 0; i < _NumOfChargedFragments; i++)
382  {
383  force.set(0.,0.,0.);
384  for (G4int j = 0; j < _NumOfChargedFragments; j++)
385  {
386  if (i != j)
387  {
388  distance = Pos[i] - Pos[j];
389  force += (elm_coupling*_theFragments[i]->GetZ()
390  *_theFragments[j]->GetZ()/
391  (distance.mag2()*distance.mag()))*distance;
392  }
393  }
394  Accel[i] = (1./(_theFragments[i]->GetNuclearMass()))*force;
395  }
396 
397  TimeN = TimeS + DeltaTime;
398 
399  for ( i = 0; i < _NumOfChargedFragments; i++)
400  {
401  SavedVel = Vel[i];
402  Vel[i] += Accel[i]*(TimeN-TimeS);
403  Pos[i] += (SavedVel+Vel[i])*(TimeN-TimeS)*0.5;
404  }
405  TimeS = TimeN;
406 
407  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
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 = 1.5*_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 
441  // Rotates a 3-vector P to close momentum triangle Pa + V + P = 0
442 {
443  G4ThreeVector U = Pa.unit();
444 
445  G4double Alpha1 = U * V;
446 
447  G4double Alpha2 = std::sqrt(V.mag2() - Alpha1*Alpha1);
448 
449  G4ThreeVector N = (1./Alpha2)*U.cross(V);
450 
451  G4ThreeVector RotatedMomentum(
452  ( (V.x() - Alpha1*U.x())/Alpha2 ) * P.x() + N.x() * P.y() + U.x() * P.z(),
453  ( (V.y() - Alpha1*U.y())/Alpha2 ) * P.x() + N.y() * P.y() + U.y() * P.z(),
454  ( (V.z() - Alpha1*U.z())/Alpha2 ) * P.x() + N.z() * P.y() + U.z() * P.z()
455  );
456  return RotatedMomentum;
457 }
458 
void set(double x, double y, double z)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetCoulombEnergy(void) const
G4double Z13(G4int Z) const
Definition: G4Pow.hh:132
G4double A13(G4double A) const
Definition: G4Pow.cc:138
G4bool CheckFragments(void)
static G4double GetKappaCoulomb()
ush Pos
Definition: deflate.h:92
const char * p
Definition: xmltok.h:285
static const G4double * P2[nN]
void CoulombImpulse(G4int anA, G4int anZ, G4double T)
double z() const
G4ThreeVector G4RandomDirection()
void setX(double)
static G4double Getr0()
Float_t Z
void CreateFragment(G4int A, G4int Z)
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4double GetFragmentsEnergy(G4double T) const
**D E S C R I P T I O N
G4ThreeVector RotateMomentum(G4ThreeVector Pa, G4ThreeVector V, G4ThreeVector P)
void FragmentsMomenta(G4int NF, G4int idx, G4double T)
void setZ(double)
static constexpr double elm_coupling
void PlaceFragments(G4int anA)
G4double GetFragmentsCoulombEnergy(void)
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
Double_t R
double A(double temperature)
std::deque< G4StatMFFragment * > _theFragments
G4double operator()(G4double &, G4StatMFFragment *&frag)
double mag2() const
G4FragmentVector * GetFragments(G4int anA, G4int anZ, G4double T)
void SolveEqOfMotion(G4int anA, G4int anZ, G4double T)
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector unit() const
static double P[]
double mag() const
int G4int
Definition: G4Types.hh:78
G4int _NumOfNeutralFragments
Definition: G4Pow.hh:56
static const G4double * P1[nN]
double x() const
double y() const
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
void setY(double)