source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroPartition.cc @ 1315

Last change on this file since 1315 was 1228, checked in by garnier, 15 years ago

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