ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMFMacroCanonical.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4StatMFMacroCanonical.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 //
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 infinite loop for
35 // a fagment with Z=A; fixed memory leak
36 
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4Pow.hh"
41 
42 // constructor
44 {
45 
46  // Get memory for clusters
47  _theClusters.push_back(new G4StatMFMacroNucleon); // Size 1
48  _theClusters.push_back(new G4StatMFMacroBiNucleon); // Size 2
49  _theClusters.push_back(new G4StatMFMacroTriNucleon); // Size 3
50  _theClusters.push_back(new G4StatMFMacroTetraNucleon); // Size 4
51  for (G4int i = 4; i < theFragment.GetA_asInt(); i++)
52  _theClusters.push_back(new G4StatMFMacroMultiNucleon(i+1)); // Size 5 ... A
53 
54  // Perform class initialization
55  Initialize(theFragment);
56 
57 }
58 
59 // destructor
61 {
62  // garbage collection
63  if (!_theClusters.empty())
64  {
65  std::for_each(_theClusters.begin(),_theClusters.end(),DeleteFragment());
66  }
67 }
68 
69 // Initialization method
71 {
72 
73  G4int A = theFragment.GetA_asInt();
74  G4int Z = theFragment.GetZ_asInt();
75  G4double x = 1.0 - 2.0*Z/G4double(A);
76  G4Pow* g4calc = G4Pow::GetInstance();
77 
78  // Free Internal energy at T = 0
79  __FreeInternalE0 = A*( -G4StatMFParameters::GetE0() + // Volume term (for T = 0)
80  G4StatMFParameters::GetGamma0()*x*x) // Symmetry term
81  + G4StatMFParameters::GetBeta0()*g4calc->Z23(A) + // Surface term (for T = 0)
82  0.6*elm_coupling*Z*Z/(G4StatMFParameters::Getr0()* // Coulomb term
83  g4calc->Z13(A));
84 
85  CalculateTemperature(theFragment);
86  return;
87 }
88 
90 {
91  // Excitation Energy
92  G4double U = theFragment.GetExcitationEnergy();
93 
94  G4int A = theFragment.GetA_asInt();
95  G4int Z = theFragment.GetZ_asInt();
96 
97  // Fragment Multiplicity
98  G4double FragMult = std::max((1.0+(2.31/MeV)*(U/A - 3.5*MeV))*A/100.0, 2.0);
99 
100  // Parameter Kappa
101  G4Pow* g4calc = G4Pow::GetInstance();
102  _Kappa = (1.0+elm_coupling*(g4calc->A13(FragMult)-1)/
103  (G4StatMFParameters::Getr0()*g4calc->Z13(A)));
104  _Kappa = _Kappa*_Kappa*_Kappa - 1.0;
105 
106  G4StatMFMacroTemperature * theTemp = new
108 
109  __MeanTemperature = theTemp->CalcTemperature();
113  __MeanEntropy = theTemp->GetEntropy();
114 
115  delete theTemp;
116 
117  return;
118 }
119 
120 // --------------------------------------------------------------------------
121 
123 // Calculate total fragments multiplicity, fragment atomic numbers and charges
124 {
125  G4int A = theFragment.GetA_asInt();
126  G4int Z = theFragment.GetZ_asInt();
127 
128  std::vector<G4int> ANumbers(A);
129 
130  G4double Multiplicity = ChooseA(A,ANumbers);
131 
132  std::vector<G4int> FragmentsA;
133 
134  G4int i = 0;
135  for (i = 0; i < A; i++)
136  {
137  for (G4int j = 0; j < ANumbers[i]; j++) FragmentsA.push_back(i+1);
138  }
139 
140  // Sort fragments in decreasing order
141  G4int im = 0;
142  for (G4int j = 0; j < Multiplicity; j++)
143  {
144  G4int FragmentsAMax = 0;
145  im = j;
146  for (i = j; i < Multiplicity; i++)
147  {
148  if (FragmentsA[i] <= FragmentsAMax) { continue; }
149  else
150  {
151  im = i;
152  FragmentsAMax = FragmentsA[im];
153  }
154  }
155  if (im != j)
156  {
157  FragmentsA[im] = FragmentsA[j];
158  FragmentsA[j] = FragmentsAMax;
159  }
160  }
161  return ChooseZ(Z,FragmentsA);
162 }
163 
164 G4double G4StatMFMacroCanonical::ChooseA(G4int A, std::vector<G4int> & ANumbers)
165  // Determines fragments multiplicities and compute total fragment multiplicity
166 {
167  G4double multiplicity = 0.0;
168  G4int i;
169 
170  std::vector<G4double> AcumMultiplicity;
171  AcumMultiplicity.reserve(A);
172 
173  AcumMultiplicity.push_back((*(_theClusters.begin()))->GetMeanMultiplicity());
174  for (std::vector<G4VStatMFMacroCluster*>::iterator it = _theClusters.begin()+1;
175  it != _theClusters.end(); ++it)
176  {
177  AcumMultiplicity.push_back((*it)->GetMeanMultiplicity()+AcumMultiplicity.back());
178  }
179 
180  G4int CheckA;
181  do {
182  CheckA = -1;
183  G4int SumA = 0;
184  G4int ThisOne = 0;
185  multiplicity = 0.0;
186  for (i = 0; i < A; i++) ANumbers[i] = 0;
187  do {
189  for (i = 0; i < A; i++) {
190  if (RandNumber < AcumMultiplicity[i]) {
191  ThisOne = i;
192  break;
193  }
194  }
195  multiplicity++;
196  ANumbers[ThisOne] = ANumbers[ThisOne]+1;
197  SumA += ThisOne+1;
198  CheckA = A - SumA;
199 
200  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
201  } while (CheckA > 0);
202 
203  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
204  } while (CheckA < 0 || std::abs(__MeanMultiplicity - multiplicity) > std::sqrt(__MeanMultiplicity) + 0.5);
205 
206  return multiplicity;
207 }
208 
210  std::vector<G4int> & FragmentsA)
211  //
212 {
213  G4Pow* g4calc = G4Pow::GetInstance();
214  std::vector<G4int> FragmentsZ;
215 
216  G4int DeltaZ = 0;
218  G4int multiplicity = FragmentsA.size();
219 
220  do {
221  FragmentsZ.clear();
222  G4int SumZ = 0;
223  for (G4int i = 0; i < multiplicity; i++)
224  {
225  G4int A = FragmentsA[i];
226  if (A <= 1)
227  {
228  G4double RandNumber = G4UniformRand();
229  if (RandNumber < (*_theClusters.begin())->GetZARatio())
230  {
231  FragmentsZ.push_back(1);
232  SumZ += FragmentsZ[i];
233  }
234  else FragmentsZ.push_back(0);
235  }
236  else
237  {
238  G4double RandZ;
240  + 2*CP*g4calc->Z23(FragmentsA[i]);
241  G4double ZMean;
242  if (FragmentsA[i] > 1 && FragmentsA[i] < 5) { ZMean = 0.5*FragmentsA[i]; }
243  else {
244  ZMean = FragmentsA[i]*(4.0*G4StatMFParameters::GetGamma0()
245  + _ChemPotentialNu)/CC;
246  }
247  G4double ZDispersion = std::sqrt(FragmentsA[i]*__MeanTemperature/CC);
248  G4int z;
249  do
250  {
251  RandZ = G4RandGauss::shoot(ZMean,ZDispersion);
252  z = G4lrint(RandZ+0.5);
253  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
254  } while (z < 0 || z > A);
255  FragmentsZ.push_back(z);
256  SumZ += z;
257  }
258  }
259  DeltaZ = Z - SumZ;
260  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
261  } while (std::abs(DeltaZ) > 1);
262 
263  // DeltaZ can be 0, 1 or -1
264  G4int idx = 0;
265  if (DeltaZ < 0.0)
266  {
267  while (FragmentsZ[idx] < 1) { ++idx; }
268  }
269  FragmentsZ[idx] += DeltaZ;
270 
271  G4StatMFChannel * theChannel = new G4StatMFChannel;
272  for (G4int i = multiplicity-1; i >= 0; i--)
273  {
274  theChannel->CreateFragment(FragmentsA[i],FragmentsZ[i]);
275  }
276 
277  return theChannel;
278 }