source: trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroCanonical.cc @ 1340

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

update ti head
File size: 9.1 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: G4StatMFMacroCanonical.cc,v 1.8 2008/11/19 14:33:31 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-03-ref-09 $
29//
30// by V. Lara
31// --------------------------------------------------------------------
32//
33// Modified:
34// 25.07.08 I.Pshenichnov (in collaboration with Alexander Botvina and Igor
35//          Mishustin (FIAS, Frankfurt, INR, Moscow and Kurchatov Institute,
36//          Moscow, pshenich@fias.uni-frankfurt.de) fixed infinite loop for
37//          a fagment with Z=A; fixed memory leak
38
39#include "G4StatMFMacroCanonical.hh"
40
41
42// constructor
43G4StatMFMacroCanonical::G4StatMFMacroCanonical(const G4Fragment & theFragment) 
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(); i++)   
52    _theClusters.push_back(new G4StatMFMacroMultiNucleon(i+1)); // Size 5 ... A
53 
54  // Perform class initialization
55  Initialize(theFragment);
56   
57}
58
59
60// destructor
61G4StatMFMacroCanonical::~G4StatMFMacroCanonical() 
62{
63  // garbage collection
64  if (!_theClusters.empty()) 
65    {
66      std::for_each(_theClusters.begin(),_theClusters.end(),DeleteFragment());
67    }
68}
69
70// operators definitions
71G4StatMFMacroCanonical & 
72G4StatMFMacroCanonical::operator=(const G4StatMFMacroCanonical & ) 
73{
74  throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroCanonical::operator= meant to not be accessable");
75  return *this;
76}
77
78G4bool G4StatMFMacroCanonical::operator==(const G4StatMFMacroCanonical & ) const 
79{
80  return false;
81}
82
83
84G4bool G4StatMFMacroCanonical::operator!=(const G4StatMFMacroCanonical & ) const 
85{
86  return true;
87}
88
89
90// Initialization method
91
92
93void G4StatMFMacroCanonical::Initialize(const G4Fragment & theFragment) 
94{
95 
96  G4double A = theFragment.GetA();
97  G4double Z = theFragment.GetZ();
98 
99  // Free Internal energy at T = 0
100  __FreeInternalE0 = A*( -G4StatMFParameters::GetE0() +          // Volume term (for T = 0)
101                         G4StatMFParameters::GetGamma0()*        // Symmetry term
102                         (1.0-2.0*Z/A)*(1.0-2.0*Z/A) ) +
103    G4StatMFParameters::GetBeta0()*std::pow(A,2.0/3.0) +              // Surface term (for T = 0)
104    (3.0/5.0)*elm_coupling*Z*Z/(G4StatMFParameters::Getr0()*     // Coulomb term
105                                std::pow(A,1.0/3.0));
106 
107 
108 
109  CalculateTemperature(theFragment);
110 
111  return;
112}
113
114
115
116
117
118void G4StatMFMacroCanonical::CalculateTemperature(const G4Fragment & theFragment)
119{
120  // Excitation Energy
121  G4double U = theFragment.GetExcitationEnergy();
122 
123  G4double A = theFragment.GetA();
124  G4double Z = theFragment.GetZ();
125 
126  // Fragment Multiplicity
127  G4double FragMult = std::max((1.0+(2.31/MeV)*(U/A - 3.5*MeV))*A/100.0, 2.0);
128
129
130  // Parameter Kappa
131  _Kappa = (1.0+elm_coupling*(std::pow(FragMult,1./3.)-1)/
132            (G4StatMFParameters::Getr0()*std::pow(A,1./3.)));
133  _Kappa = _Kappa*_Kappa*_Kappa - 1.0;
134
135       
136  G4StatMFMacroTemperature * theTemp = new     
137    G4StatMFMacroTemperature(A,Z,U,__FreeInternalE0,_Kappa,&_theClusters);
138 
139  __MeanTemperature = theTemp->CalcTemperature();
140  _ChemPotentialNu = theTemp->GetChemicalPotentialNu();
141  _ChemPotentialMu = theTemp->GetChemicalPotentialMu();
142  __MeanMultiplicity = theTemp->GetMeanMultiplicity();
143  __MeanEntropy = theTemp->GetEntropy();
144       
145  delete theTemp;                       
146 
147  return;
148}
149
150
151// --------------------------------------------------------------------------
152
153G4StatMFChannel * G4StatMFMacroCanonical::ChooseAandZ(const G4Fragment &theFragment)
154    // Calculate total fragments multiplicity, fragment atomic numbers and charges
155{
156  G4double A = theFragment.GetA();
157  G4double Z = theFragment.GetZ();
158 
159  std::vector<G4double> ANumbers(static_cast<G4int>(A));
160 
161  G4double Multiplicity = ChooseA(A,ANumbers);
162 
163 
164  std::vector<G4double> FragmentsA;
165 
166  G4int i = 0; 
167    for (i = 0; i < A; i++) 
168      {
169        for (G4int j = 0; j < ANumbers[i]; j++) FragmentsA.push_back(i+1);
170      }
171
172   
173    // Sort fragments in decreasing order
174    G4int im = 0;
175    for (G4int j = 0; j < Multiplicity; j++) 
176      {
177        G4double FragmentsAMax = 0.0;
178        im = j;
179        for (i = j; i < Multiplicity; i++) 
180          {
181            if (FragmentsA[i] <= FragmentsAMax) continue;
182            else 
183              {
184                im = i;
185                FragmentsAMax = FragmentsA[im];
186              }
187          }
188       
189        if (im != j) 
190          {
191            FragmentsA[im] = FragmentsA[j];
192            FragmentsA[j] = FragmentsAMax;
193          }
194      }
195   
196    return ChooseZ(static_cast<G4int>(Z),FragmentsA);
197}
198
199
200
201G4double G4StatMFMacroCanonical::ChooseA(const G4double A, std::vector<G4double> & ANumbers)
202  // Determines fragments multiplicities and compute total fragment multiplicity
203{
204  G4double multiplicity = 0.0;
205  G4int i;
206 
207 
208  std::vector<G4double> AcumMultiplicity;
209  AcumMultiplicity.reserve(static_cast<G4int>(A));
210
211  AcumMultiplicity.push_back((*(_theClusters.begin()))->GetMeanMultiplicity());
212  for (std::vector<G4VStatMFMacroCluster*>::iterator it = _theClusters.begin()+1;
213       it != _theClusters.end(); ++it)
214    {
215      AcumMultiplicity.push_back((*it)->GetMeanMultiplicity()+AcumMultiplicity.back());
216    }
217 
218  G4int CheckA;
219  do {
220    CheckA = -1;
221    G4int SumA = 0;
222    G4int ThisOne = 0;
223    multiplicity = 0.0;
224    for (i = 0; i < A; i++) ANumbers[i] = 0.0;
225    do {
226      G4double RandNumber = G4UniformRand()*__MeanMultiplicity;
227      for (i = 0; i < A; i++) {
228        if (RandNumber < AcumMultiplicity[i]) {
229          ThisOne = i;
230          break;
231        }
232      }
233      multiplicity++;
234      ANumbers[ThisOne] = ANumbers[ThisOne]+1;
235      SumA += ThisOne+1;
236      CheckA = static_cast<G4int>(A) - SumA;
237     
238    } while (CheckA > 0);
239   
240  } while (CheckA < 0 || std::abs(__MeanMultiplicity - multiplicity) > std::sqrt(__MeanMultiplicity) + 1./2.);
241 
242  return multiplicity;
243}
244
245
246G4StatMFChannel * G4StatMFMacroCanonical::ChooseZ(const G4int & Z, 
247                                                  std::vector<G4double> & FragmentsA)
248    //
249{
250  std::vector<G4double> FragmentsZ;
251 
252  G4double DeltaZ = 0.0;
253  G4double CP = (3./5.)*(elm_coupling/G4StatMFParameters::Getr0())*
254    (1.0 - 1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.));
255 
256  G4int multiplicity = FragmentsA.size();
257 
258  do 
259    {
260      FragmentsZ.clear();
261      G4int SumZ = 0;
262      for (G4int i = 0; i < multiplicity; i++) 
263        {
264          G4double A = FragmentsA[i];
265          if (A <= 1.0) 
266            {
267              G4double RandNumber = G4UniformRand();
268              if (RandNumber < (*_theClusters.begin())->GetZARatio()) 
269                {
270                  FragmentsZ.push_back(1.0);
271                  SumZ += static_cast<G4int>(FragmentsZ[i]);
272                } 
273              else FragmentsZ.push_back(0.0);
274            } 
275          else 
276            {
277              G4double RandZ;
278              G4double CC = 8.0*G4StatMFParameters::GetGamma0()+2.0*CP*std::pow(FragmentsA[i],2./3.);
279              G4double ZMean;
280              if (FragmentsA[i] > 1.5 && FragmentsA[i] < 4.5) ZMean = 0.5*FragmentsA[i];
281              else ZMean = FragmentsA[i]*(4.0*G4StatMFParameters::GetGamma0()+_ChemPotentialNu)/CC;
282              G4double ZDispersion = std::sqrt(FragmentsA[i]*__MeanTemperature/CC);
283              G4int z;
284              do 
285                {
286                  RandZ = G4RandGauss::shoot(ZMean,ZDispersion);
287                  z = static_cast<G4int>(RandZ+0.5);
288                } while (z < 0 || z > A);
289              FragmentsZ.push_back(z);
290              SumZ += z;
291            }
292        }
293      DeltaZ = Z - SumZ;
294    }
295  while (std::abs(DeltaZ) > 1.1);
296   
297  // DeltaZ can be 0, 1 or -1
298  G4int idx = 0;
299  if (DeltaZ < 0.0)
300    {
301      while (FragmentsZ[idx] < 0.5) ++idx;
302    }
303  FragmentsZ[idx] += DeltaZ;
304 
305  G4StatMFChannel * theChannel = new G4StatMFChannel;
306  for (G4int i = multiplicity-1; i >= 0; i--) 
307    {
308      theChannel->CreateFragment(FragmentsA[i],FragmentsZ[i]);
309    }
310 
311   
312    return theChannel;
313}
Note: See TracBrowser for help on using the repository browser.