ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMFMicroPartition.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4StatMFMicroPartition.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 // by V. Lara
29 // --------------------------------------------------------------------
30 
32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4HadronicException.hh"
35 #include "Randomize.hh"
36 #include "G4Log.hh"
37 #include "G4Exp.hh"
38 #include "G4Pow.hh"
39 
40 // Copy constructor
42 {
43  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::copy_constructor meant to not be accessible");
44 }
45 
46 // Operators
47 
50 {
51  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator= meant to not be accessible");
52  return *this;
53 }
54 
55 
57 {
58  //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator== meant to not be accessible");
59  return false;
60 }
61 
62 
64 {
65  //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator!= meant to not be accessible");
66  return true;
67 }
68 
70 {
71  // This Z independent factor in the Coulomb free energy
72  G4double CoulombConstFactor = G4StatMFParameters::GetCoulomb();
73 
74  // We use the aproximation Z_f ~ Z/A * A_f
75 
77 
78  if (anA == 0 || anA == 1)
79  {
80  _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA);
81  }
82  else if (anA == 2 || anA == 3 || anA == 4)
83  {
84  // Z/A ~ 1/2
85  _theCoulombFreeEnergy.push_back(CoulombConstFactor*0.5
86  *anA*G4Pow::GetInstance()->Z23(anA));
87  }
88  else // anA > 4
89  {
90  _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA
91  *anA*G4Pow::GetInstance()->Z23(anA));
92  }
93 }
94 
96 {
97  G4Pow* g4calc = G4Pow::GetInstance();
98  G4double CoulombFactor = 1.0/g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb());
99 
100  G4double CoulombEnergy = elm_coupling*0.6*theZ*theZ*CoulombFactor/
101  (G4StatMFParameters::Getr0()*g4calc->Z13(theA));
102 
104  for (unsigned int i = 0; i < _thePartition.size(); i++)
105  CoulombEnergy += _theCoulombFreeEnergy[i] - elm_coupling*0.6*
106  ZA*ZA*_thePartition[i]*g4calc->Z23(_thePartition[i])/
108 
109  return CoulombEnergy;
110 }
111 
113 {
114  G4Pow* g4calc = G4Pow::GetInstance();
115  G4double CoulombFactor = 1.0/g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb());
116 
117  G4double PartitionEnergy = 0.0;
118 
119  // We use the aprox that Z_f ~ Z/A * A_f
120  for (unsigned int i = 0; i < _thePartition.size(); i++)
121  {
122  if (_thePartition[i] == 0 || _thePartition[i] == 1)
123  {
124  PartitionEnergy += _theCoulombFreeEnergy[i];
125  }
126  else if (_thePartition[i] == 2)
127  {
128  PartitionEnergy +=
129  -2.796 // Binding Energy of deuteron ??????
130  + _theCoulombFreeEnergy[i];
131  }
132  else if (_thePartition[i] == 3)
133  {
134  PartitionEnergy +=
135  -9.224 // Binding Energy of trtion/He3 ??????
136  + _theCoulombFreeEnergy[i];
137  }
138  else if (_thePartition[i] == 4)
139  {
140  PartitionEnergy +=
141  -30.11 // Binding Energy of ALPHA ??????
143  + 4.*T*T/InvLevelDensity(4.);
144  }
145  else
146  {
147  PartitionEnergy +=
148  //Volume term
151  *_thePartition[i] +
152 
153  // Symmetry term
155  (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/theA)*_thePartition[i] +
156 
157  // Surface term
159  g4calc->Z23(_thePartition[i]) +
160 
161  // Coulomb term
163  }
164  }
165 
166  PartitionEnergy += elm_coupling*0.6*theZ*theZ*CoulombFactor/
167  (G4StatMFParameters::Getr0()*g4calc->Z13(theA))
168  + 1.5*T*(_thePartition.size()-1);
169 
170  return PartitionEnergy;
171 }
172 
174  G4double FreeInternalE0)
175 {
176  G4double PartitionEnergy = GetPartitionEnergy(0.0);
177 
178  // If this happens, T = 0 MeV, which means that probability for this
179  // partition will be 0
180  if (std::fabs(U + FreeInternalE0 - PartitionEnergy) < 0.003) return -1.0;
181 
182  // Calculate temperature by midpoint method
183 
184  // Bracketing the solution
185  G4double Ta = 0.001;
186  G4double Tb = std::max(std::sqrt(8.0*U/theA),0.0012*MeV);
187  G4double Tmid = 0.0;
188 
189  G4double Da = (U + FreeInternalE0 - GetPartitionEnergy(Ta))/U;
190  G4double Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
191 
192  G4int maxit = 0;
193  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
194  while (Da*Db > 0.0 && maxit < 1000)
195  {
196  ++maxit;
197  Tb += 0.5*Tb;
198  Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
199  }
200 
201  G4double eps = 1.0e-14*std::abs(Ta-Tb);
202 
203  for (G4int i = 0; i < 1000; i++)
204  {
205  Tmid = (Ta+Tb)/2.0;
206  if (std::fabs(Ta-Tb) <= eps) return Tmid;
207  G4double Dmid = (U + FreeInternalE0 - GetPartitionEnergy(Tmid))/U;
208  if (std::fabs(Dmid) < 0.003) return Tmid;
209  if (Da*Dmid < 0.0)
210  {
211  Tb = Tmid;
212  Db = Dmid;
213  }
214  else
215  {
216  Ta = Tmid;
217  Da = Dmid;
218  }
219  }
220  // if we arrive here the temperature could not be calculated
221  G4cout << "G4StatMFMicroPartition::CalcPartitionTemperature: I can't calculate the temperature"
222  << G4endl;
223  // and set probability to 0 returning T < 0
224  return -1.0;
225 
226 }
227 
229  G4double FreeInternalE0,
230  G4double SCompound)
231 {
232  G4double T = CalcPartitionTemperature(U,FreeInternalE0);
233  if ( T <= 0.0) return _Probability = 0.0;
234  _Temperature = T;
235 
236  G4Pow* g4calc = G4Pow::GetInstance();
237 
238  // Factorial of fragment multiplicity
239  G4double Fact = 1.0;
240  unsigned int i;
241  for (i = 0; i < _thePartition.size() - 1; i++)
242  {
243  G4double f = 1.0;
244  for (unsigned int ii = i+1; i< _thePartition.size(); i++)
245  {
246  if (_thePartition[i] == _thePartition[ii]) f++;
247  }
248  Fact *= f;
249  }
250 
251  G4double ProbDegeneracy = 1.0;
252  G4double ProbA32 = 1.0;
253 
254  for (i = 0; i < _thePartition.size(); i++)
255  {
256  ProbDegeneracy *= GetDegeneracyFactor(_thePartition[i]);
257  ProbA32 *= _thePartition[i]*std::sqrt((G4double)_thePartition[i]);
258  }
259 
260  // Compute entropy
261  G4double PartitionEntropy = 0.0;
262  for (i = 0; i < _thePartition.size(); i++)
263  {
264  // interaction entropy for alpha
265  if (_thePartition[i] == 4)
266  {
267  PartitionEntropy +=
269  }
270  // interaction entropy for Af > 4
271  else if (_thePartition[i] > 4)
272  {
273  PartitionEntropy +=
275  - G4StatMFParameters::DBetaDT(T) * g4calc->Z23(_thePartition[i]);
276  }
277  }
278 
279  // Thermal Wave Lenght = std::sqrt(2 pi hbar^2 / nucleon_mass T)
280  G4double ThermalWaveLenght3 = 16.15*fermi/std::sqrt(T);
281  ThermalWaveLenght3 = ThermalWaveLenght3*ThermalWaveLenght3*ThermalWaveLenght3;
282 
283  // Translational Entropy
284  G4double kappa = 1. + elm_coupling*(g4calc->Z13(_thePartition.size())-1.0)
285  /(G4StatMFParameters::Getr0()*g4calc->Z13(theA));
286  kappa = kappa*kappa*kappa;
287  kappa -= 1.;
290  G4double FreeVolume = kappa*V0;
291  G4double TranslationalS = std::max(0.0, G4Log(ProbA32/Fact) +
292  (_thePartition.size()-1.0)*G4Log(FreeVolume/ThermalWaveLenght3) +
293  1.5*(_thePartition.size()-1.0) - 1.5*g4calc->logZ(theA));
294 
295  PartitionEntropy += G4Log(ProbDegeneracy) + TranslationalS;
296  _Entropy = PartitionEntropy;
297 
298  // And finally compute probability of fragment configuration
299  G4double exponent = PartitionEntropy-SCompound;
300  if (exponent > 300.0) exponent = 300.0;
301  return _Probability = G4Exp(exponent);
302 }
303 
305 {
306  // Degeneracy factors are statistical factors
307  // DegeneracyFactor for nucleon is (2S_n + 1)(2I_n + 1) = 4
308  G4double DegFactor = 0;
309  if (A > 4) DegFactor = 1.0;
310  else if (A == 1) DegFactor = 4.0; // nucleon
311  else if (A == 2) DegFactor = 3.0; // Deuteron
312  else if (A == 3) DegFactor = 4.0; // Triton + He3
313  else if (A == 4) DegFactor = 1.0; // alpha
314  return DegFactor;
315 }
316 
318 // Gives fragments charges
319 {
320  std::vector<G4int> FragmentsZ;
321 
322  G4int ZBalance = 0;
323  do
324  {
326  G4int SumZ = 0;
327  for (unsigned int i = 0; i < _thePartition.size(); i++)
328  {
329  G4double ZMean;
330  G4double Af = _thePartition[i];
331  if (Af > 1.5 && Af < 4.5) ZMean = 0.5*Af;
332  else ZMean = Af*Z0/A0;
333  G4double ZDispersion = std::sqrt(Af * MeanT/CC);
334  G4int Zf;
335  do
336  {
337  Zf = static_cast<G4int>(G4RandGauss::shoot(ZMean,ZDispersion));
338  }
339  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
340  while (Zf < 0 || Zf > Af);
341  FragmentsZ.push_back(Zf);
342  SumZ += Zf;
343  }
344  ZBalance = Z0 - SumZ;
345  }
346  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
347  while (std::abs(ZBalance) > 1);
348  FragmentsZ[0] += ZBalance;
349 
350  G4StatMFChannel * theChannel = new G4StatMFChannel;
351  for (unsigned int i = 0; i < _thePartition.size(); i++)
352  {
353  theChannel->CreateFragment(_thePartition[i],FragmentsZ[i]);
354  }
355 
356  return theChannel;
357 }