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