Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4StatMFMicroCanonical.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$
28 //
29 // Hadronic Process: Nuclear De-excitations
30 // by V. Lara
31 
32 #include <numeric>
33 
35 #include "G4PhysicalConstants.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4HadronicException.hh"
38 #include "G4Pow.hh"
39 
40 // constructor
42 {
43  // Perform class initialization
44  Initialize(theFragment);
45 
46 }
47 
48 
49 // destructor
51 {
52  // garbage collection
53  if (!_ThePartitionManagerVector.empty()) {
54  std::for_each(_ThePartitionManagerVector.begin(),
55  _ThePartitionManagerVector.end(),
56  DeleteFragment());
57  }
58 }
59 
60 
61 
62 // Initialization method
63 
64 void G4StatMFMicroCanonical::Initialize(const G4Fragment & theFragment)
65 {
66 
67  std::vector<G4StatMFMicroManager*>::iterator it;
68 
69  // Excitation Energy
70  G4double U = theFragment.GetExcitationEnergy();
71 
72  G4int A = theFragment.GetA_asInt();
73  G4int Z = theFragment.GetZ_asInt();
74  G4double x = 1.0 - 2.0*Z/G4double(A);
75  G4Pow* g4pow = G4Pow::GetInstance();
76 
77  // Configuration temperature
78  G4double TConfiguration = std::sqrt(8.0*U/G4double(A));
79 
80  // Free internal energy at Temperature T = 0
81  __FreeInternalE0 = A*(
82  // Volume term (for T = 0)
84  // Symmetry term
86  ) +
87  // Surface term (for T = 0)
88  G4StatMFParameters::GetBeta0()*g4pow->Z23(A) +
89  // Coulomb term
90  elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*g4pow->Z13(A));
91 
92  // Statistical weight
93  G4double W = 0.0;
94 
95 
96  // Mean breakup multiplicity
97  __MeanMultiplicity = 0.0;
98 
99  // Mean channel temperature
100  __MeanTemperature = 0.0;
101 
102  // Mean channel entropy
103  __MeanEntropy = 0.0;
104 
105  // Calculate entropy of compound nucleus
106  G4double SCompoundNucleus = CalcEntropyOfCompoundNucleus(theFragment,TConfiguration);
107 
108  // Statistical weight of compound nucleus
109  _WCompoundNucleus = 1.0; // std::exp(SCompoundNucleus - SCompoundNucleus);
110 
111  W += _WCompoundNucleus;
112 
113 
114 
115  // Maximal fragment multiplicity allowed in direct simulation
117  if (A > 110) MaxMult -= 1;
118 
119 
120 
121  for (G4int im = 2; im <= MaxMult; im++) {
122  G4StatMFMicroManager * aMicroManager =
123  new G4StatMFMicroManager(theFragment,im,__FreeInternalE0,SCompoundNucleus);
124  _ThePartitionManagerVector.push_back(aMicroManager);
125  }
126 
127  // W is the total probability
128  W = std::accumulate(_ThePartitionManagerVector.begin(),
129  _ThePartitionManagerVector.end(),
130  W,SumProbabilities());
131 
132  // Normalization of statistical weights
133  for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
134  {
135  (*it)->Normalize(W);
136  }
137 
138  _WCompoundNucleus /= W;
139 
140  __MeanMultiplicity += 1.0 * _WCompoundNucleus;
141  __MeanTemperature += TConfiguration * _WCompoundNucleus;
142  __MeanEntropy += SCompoundNucleus * _WCompoundNucleus;
143 
144 
145  for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it)
146  {
147  __MeanMultiplicity += (*it)->GetMeanMultiplicity();
148  __MeanTemperature += (*it)->GetMeanTemperature();
149  __MeanEntropy += (*it)->GetMeanEntropy();
150  }
151 
152  return;
153 }
154 
155 
156 G4double G4StatMFMicroCanonical::CalcFreeInternalEnergy(const G4Fragment & theFragment,
157  G4double T)
158 {
159  G4int A = theFragment.GetA_asInt();
160  G4int Z = theFragment.GetZ_asInt();
162 
163  G4double InvLevelDensityPar = G4StatMFParameters::GetEpsilon0()*(1.0 + 3.0/G4double(A-1));
164 
165  G4double VolumeTerm = (-G4StatMFParameters::GetE0()+T*T/InvLevelDensityPar)*A;
166 
167  G4double SymmetryTerm = G4StatMFParameters::GetGamma0()*(A - 2*Z)*(A - 2*Z)/G4double(A);
168 
169  G4double SurfaceTerm = (G4StatMFParameters::Beta(T)-T*G4StatMFParameters::DBetaDT(T))*A13*A13;
170 
171  G4double CoulombTerm = elm_coupling*(3.0/5.0)*Z*Z/(G4StatMFParameters::Getr0()*A13);
172 
173  return VolumeTerm + SymmetryTerm + SurfaceTerm + CoulombTerm;
174 }
175 
176 G4double
177 G4StatMFMicroCanonical::CalcEntropyOfCompoundNucleus(const G4Fragment & theFragment,
178  G4double & TConf)
179  // Calculates Temperature and Entropy of compound nucleus
180 {
181  G4int A = theFragment.GetA_asInt();
182  // const G4double Z = theFragment.GetZ();
183  G4double U = theFragment.GetExcitationEnergy();
184  G4double A13 = G4Pow::GetInstance()->Z13(A);
185 
186  G4double Ta = std::max(std::sqrt(U/(0.125*A)),0.0012*MeV);
187  G4double Tb = Ta;
188 
189  G4double ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Ta);
190  G4double Da = (U+__FreeInternalE0-ECompoundNucleus)/U;
191  G4double Db = 0.0;
192 
193  G4double InvLevelDensity = CalcInvLevelDensity(static_cast<G4int>(A));
194 
195  // bracketing the solution
196  if (Da == 0.0) {
197  TConf = Ta;
198  return 2*Ta*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Ta)*A13*A13;
199  } else if (Da < 0.0) {
200  do {
201  Tb -= 0.5*Tb;
202  ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
203  Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
204  } while (Db < 0.0);
205  } else {
206  do {
207  Tb += 0.5*Tb;
208  ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tb);
209  Db = (U+__FreeInternalE0-ECompoundNucleus)/U;
210  } while (Db > 0.0);
211  }
212 
213  G4double eps = 1.0e-14 * std::fabs(Tb-Ta);
214 
215  for (G4int i = 0; i < 1000; i++) {
216  G4double Tc = (Ta+Tb)/2.0;
217  if (std::abs(Ta-Tb) <= eps) {
218  TConf = Tc;
219  return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
220  }
221  ECompoundNucleus = CalcFreeInternalEnergy(theFragment,Tc);
222  G4double Dc = (U+__FreeInternalE0-ECompoundNucleus)/U;
223 
224  if (Dc == 0.0) {
225  TConf = Tc;
226  return 2*Tc*A/InvLevelDensity - G4StatMFParameters::DBetaDT(Tc)*A13*A13;
227  }
228 
229  if (Da*Dc < 0.0) {
230  Tb = Tc;
231  Db = Dc;
232  } else {
233  Ta = Tc;
234  Da = Dc;
235  }
236  }
237 
238  G4cerr <<
239  "G4StatMFMicrocanoncal::CalcEntropyOfCompoundNucleus: I can't calculate the temperature"
240  << G4endl;
241 
242  return 0.0;
243 }
244 
246  // Choice of fragment atomic numbers and charges
247 {
248  // We choose a multiplicity (1,2,3,...) and then a channel
249  G4double RandNumber = G4UniformRand();
250 
251  if (RandNumber < _WCompoundNucleus) {
252 
253  G4StatMFChannel * aChannel = new G4StatMFChannel;
254  aChannel->CreateFragment(theFragment.GetA_asInt(),theFragment.GetZ_asInt());
255  return aChannel;
256 
257  } else {
258 
259  G4double AccumWeight = _WCompoundNucleus;
260  std::vector<G4StatMFMicroManager*>::iterator it;
261  for (it = _ThePartitionManagerVector.begin(); it != _ThePartitionManagerVector.end(); ++it) {
262  AccumWeight += (*it)->GetProbability();
263  if (RandNumber < AccumWeight) {
264  return (*it)->ChooseChannel(theFragment.GetA(),theFragment.GetZ(),__MeanTemperature);
265  }
266  }
267  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroCanonical::ChooseAandZ: wrong normalization!");
268  }
269 
270  return 0;
271 }
272 
273 G4double G4StatMFMicroCanonical::CalcInvLevelDensity(G4int anA)
274 {
275  // Calculate Inverse Density Level
276  // Epsilon0*(1 + 3 /(Af - 1))
277  if (anA == 1) return 0.0;
278  else return
279  G4StatMFParameters::GetEpsilon0()*(1.0+3.0/(anA - 1.0));
280 }