Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroPartition.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroPartition.cc (Version 11.3.0) and /processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMicroPartition.cc (Version 10.7.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //                                                 27 //
 28 // by V. Lara                                      28 // by V. Lara
 29 // -------------------------------------------     29 // --------------------------------------------------------------------
 30                                                    30 
 31 #include "G4StatMFMicroPartition.hh"               31 #include "G4StatMFMicroPartition.hh"
 32 #include "G4PhysicalConstants.hh"                  32 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"                      33 #include "G4SystemOfUnits.hh"
 34 #include "G4HadronicException.hh"                  34 #include "G4HadronicException.hh"
 35 #include "Randomize.hh"                            35 #include "Randomize.hh"
 36 #include "G4Log.hh"                                36 #include "G4Log.hh"
 37 #include "G4Exp.hh"                                37 #include "G4Exp.hh"
 38 #include "G4Pow.hh"                                38 #include "G4Pow.hh"
 39                                                    39 
 40 // Copy constructor                                40 // Copy constructor
 41 G4StatMFMicroPartition::G4StatMFMicroPartition     41 G4StatMFMicroPartition::G4StatMFMicroPartition(const G4StatMFMicroPartition & )
 42 {                                                  42 {
 43   throw G4HadronicException(__FILE__, __LINE__     43   throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::copy_constructor meant to not be accessible");
 44 }                                                  44 }
 45                                                    45 
 46 // Operators                                       46 // Operators
 47                                                    47 
 48 G4StatMFMicroPartition & G4StatMFMicroPartitio     48 G4StatMFMicroPartition & G4StatMFMicroPartition::
 49 operator=(const G4StatMFMicroPartition & )         49 operator=(const G4StatMFMicroPartition & )
 50 {                                                  50 {
 51   throw G4HadronicException(__FILE__, __LINE__     51   throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator= meant to not be accessible");
 52   return *this;                                    52   return *this;
 53 }                                                  53 }
 54                                                    54 
 55                                                    55 
 56 G4bool G4StatMFMicroPartition::operator==(cons     56 G4bool G4StatMFMicroPartition::operator==(const G4StatMFMicroPartition & ) const
 57 {                                                  57 {
 58   //throw G4HadronicException(__FILE__, __LINE     58   //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator== meant to not be accessible");
 59   return false;                                    59   return false;
 60 }                                                  60 }
 61                                                    61  
 62                                                    62 
 63 G4bool G4StatMFMicroPartition::operator!=(cons     63 G4bool G4StatMFMicroPartition::operator!=(const G4StatMFMicroPartition & ) const
 64 {                                                  64 {
 65   //throw G4HadronicException(__FILE__, __LINE     65   //throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMicroPartition::operator!= meant to not be accessible");
 66   return true;                                     66   return true;
 67 }                                                  67 }
 68                                                    68 
 69 void G4StatMFMicroPartition::CoulombFreeEnergy     69 void G4StatMFMicroPartition::CoulombFreeEnergy(G4int anA)
 70 {                                                  70 {
 71   // This Z independent factor in the Coulomb      71   // This Z independent factor in the Coulomb free energy 
 72   G4double  CoulombConstFactor = G4StatMFParam     72   G4double  CoulombConstFactor = G4StatMFParameters::GetCoulomb();
 73                                                    73 
 74   // We use the aproximation Z_f ~ Z/A * A_f       74   // We use the aproximation Z_f ~ Z/A * A_f
 75                                                    75 
 76   G4double ZA = G4double(theZ)/G4double(theA);     76   G4double ZA = G4double(theZ)/G4double(theA);
 77                                                    77                     
 78   if (anA == 0 || anA == 1)                        78   if (anA == 0 || anA == 1) 
 79     {                                              79     {
 80       _theCoulombFreeEnergy.push_back(CoulombC     80       _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA);
 81     }                                              81     } 
 82   else if (anA == 2 || anA == 3 || anA == 4)       82   else if (anA == 2 || anA == 3 || anA == 4) 
 83     {                                              83     {
 84       // Z/A ~ 1/2                                 84       // Z/A ~ 1/2
 85       _theCoulombFreeEnergy.push_back(CoulombC     85       _theCoulombFreeEnergy.push_back(CoulombConstFactor*0.5
 86               *anA*G4Pow::GetInstance()->Z23(a     86               *anA*G4Pow::GetInstance()->Z23(anA));
 87     }                                              87     } 
 88   else  // anA > 4                                 88   else  // anA > 4
 89     {                                              89     {
 90       _theCoulombFreeEnergy.push_back(CoulombC     90       _theCoulombFreeEnergy.push_back(CoulombConstFactor*ZA*ZA  
 91               *anA*G4Pow::GetInstance()->Z23(a     91               *anA*G4Pow::GetInstance()->Z23(anA));
 92     }                                              92     }
 93 }                                                  93 }
 94                                                    94 
 95 G4double G4StatMFMicroPartition::GetCoulombEne     95 G4double G4StatMFMicroPartition::GetCoulombEnergy(void)
 96 {                                                  96 {
 97   G4Pow* g4calc = G4Pow::GetInstance();            97   G4Pow* g4calc = G4Pow::GetInstance();
 98   G4double  CoulombFactor = 1.0/g4calc->A13(1.     98   G4double  CoulombFactor = 1.0/g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb()); 
 99                                                    99       
100   G4double CoulombEnergy = elm_coupling*0.6*th    100   G4double CoulombEnergy = elm_coupling*0.6*theZ*theZ*CoulombFactor/
101     (G4StatMFParameters::Getr0()*g4calc->Z13(t    101     (G4StatMFParameters::Getr0()*g4calc->Z13(theA));
102                                                   102   
103   G4double ZA = G4double(theZ)/G4double(theA);    103   G4double ZA = G4double(theZ)/G4double(theA);
104   for (unsigned int i = 0; i < _thePartition.s    104   for (unsigned int i = 0; i < _thePartition.size(); i++) 
105     CoulombEnergy += _theCoulombFreeEnergy[i]     105     CoulombEnergy += _theCoulombFreeEnergy[i] - elm_coupling*0.6*
106       ZA*ZA*_thePartition[i]*g4calc->Z23(_theP    106       ZA*ZA*_thePartition[i]*g4calc->Z23(_thePartition[i])/
107       G4StatMFParameters::Getr0();                107       G4StatMFParameters::Getr0();
108                                                   108     
109   return CoulombEnergy;                           109   return CoulombEnergy;
110 }                                                 110 }
111                                                   111 
112 G4double G4StatMFMicroPartition::GetPartitionE    112 G4double G4StatMFMicroPartition::GetPartitionEnergy(G4double T)
113 {                                                 113 {
114   G4Pow* g4calc = G4Pow::GetInstance();           114   G4Pow* g4calc = G4Pow::GetInstance();
115   G4double  CoulombFactor = 1.0/g4calc->A13(1.    115   G4double  CoulombFactor = 1.0/g4calc->A13(1.0+G4StatMFParameters::GetKappaCoulomb()); 
116                                                   116   
117   G4double PartitionEnergy = 0.0;                 117   G4double PartitionEnergy = 0.0;
118                                                   118   
119   // We use the aprox that Z_f ~ Z/A * A_f        119   // We use the aprox that Z_f ~ Z/A * A_f
120   for (unsigned int i = 0; i < _thePartition.s    120   for (unsigned int i = 0; i < _thePartition.size(); i++) 
121     {                                             121     {
122       if (_thePartition[i] == 0 || _thePartiti    122       if (_thePartition[i] == 0 || _thePartition[i] == 1) 
123         {                                         123         { 
124           PartitionEnergy += _theCoulombFreeEn    124           PartitionEnergy += _theCoulombFreeEnergy[i];
125         }                                         125         }
126       else if (_thePartition[i] == 2)             126       else if (_thePartition[i] == 2) 
127         {                                         127         {   
128           PartitionEnergy +=                      128           PartitionEnergy +=  
129             -2.796 // Binding Energy of deuter    129             -2.796 // Binding Energy of deuteron ??????
130             + _theCoulombFreeEnergy[i];           130             + _theCoulombFreeEnergy[i];   
131   }                                               131   }
132       else if (_thePartition[i] == 3)             132       else if (_thePartition[i] == 3) 
133         {                                         133         { 
134           PartitionEnergy +=                      134           PartitionEnergy +=  
135             -9.224 // Binding Energy of trtion    135             -9.224 // Binding Energy of trtion/He3 ??????
136             + _theCoulombFreeEnergy[i];           136             + _theCoulombFreeEnergy[i];   
137   }                                               137   } 
138       else if (_thePartition[i] == 4)             138       else if (_thePartition[i] == 4) 
139         {                                         139         { 
140           PartitionEnergy +=                      140           PartitionEnergy +=
141             -30.11 // Binding Energy of ALPHA     141             -30.11 // Binding Energy of ALPHA ??????
142             + _theCoulombFreeEnergy[i]            142             + _theCoulombFreeEnergy[i] 
143             + 4.*T*T/InvLevelDensity(4.);         143             + 4.*T*T/InvLevelDensity(4.);
144   }                                               144   } 
145       else                                        145       else 
146         {                                         146         {
147           PartitionEnergy +=                      147           PartitionEnergy +=
148             //Volume term                         148             //Volume term           
149             (- G4StatMFParameters::GetE0() +      149             (- G4StatMFParameters::GetE0() + 
150              T*T/InvLevelDensity(_thePartition    150              T*T/InvLevelDensity(_thePartition[i]))
151             *_thePartition[i] +                   151             *_thePartition[i] + 
152                                                   152             
153             // Symmetry term                      153             // Symmetry term
154             G4StatMFParameters::GetGamma0()*      154             G4StatMFParameters::GetGamma0()*
155             (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/    155             (1.0-2.0*theZ/theA)*(1.0-2.0*theZ/theA)*_thePartition[i] +  
156                                                   156             
157             // Surface term                       157             // Surface term
158             (G4StatMFParameters::Beta(T) - T*G    158             (G4StatMFParameters::Beta(T) - T*G4StatMFParameters::DBetaDT(T))*
159             g4calc->Z23(_thePartition[i]) +       159             g4calc->Z23(_thePartition[i]) +
160                                                   160             
161             // Coulomb term                       161             // Coulomb term 
162             _theCoulombFreeEnergy[i];             162             _theCoulombFreeEnergy[i];
163   }                                               163   }
164     }                                             164     }
165                                                   165   
166   PartitionEnergy += elm_coupling*0.6*theZ*the    166   PartitionEnergy += elm_coupling*0.6*theZ*theZ*CoulombFactor/
167     (G4StatMFParameters::Getr0()*g4calc->Z13(t    167     (G4StatMFParameters::Getr0()*g4calc->Z13(theA))
168     + 1.5*T*(_thePartition.size()-1);             168     + 1.5*T*(_thePartition.size()-1);
169                                                   169   
170   return PartitionEnergy;                         170   return PartitionEnergy;
171 }                                                 171 }
172                                                   172 
173 G4double G4StatMFMicroPartition::CalcPartition    173 G4double G4StatMFMicroPartition::CalcPartitionTemperature(G4double U,
174                 G4double FreeInternalE0)          174                 G4double FreeInternalE0)
175 {                                                 175 {
176   G4double PartitionEnergy = GetPartitionEnerg    176   G4double PartitionEnergy = GetPartitionEnergy(0.0);
177                                                   177   
178   // If this happens, T = 0 MeV, which means t    178   // If this happens, T = 0 MeV, which means that probability for this
179   // partition will be 0                          179   // partition will be 0
180   if (std::fabs(U + FreeInternalE0 - Partition    180   if (std::fabs(U + FreeInternalE0 - PartitionEnergy) < 0.003) return -1.0;
181                                                   181     
182   // Calculate temperature by midpoint method     182   // Calculate temperature by midpoint method
183                                                   183   
184   // Bracketing the solution                      184   // Bracketing the solution
185   G4double Ta = 0.001;                            185   G4double Ta = 0.001;
186   G4double Tb = std::max(std::sqrt(8.0*U/theA)    186   G4double Tb = std::max(std::sqrt(8.0*U/theA),0.0012*MeV);
187   G4double Tmid = 0.0;                            187   G4double Tmid = 0.0;
188                                                   188   
189   G4double Da = (U + FreeInternalE0 - GetParti    189   G4double Da = (U + FreeInternalE0 - GetPartitionEnergy(Ta))/U;
190   G4double Db = (U + FreeInternalE0 - GetParti    190   G4double Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
191                                                   191   
192   G4int maxit = 0;                                192   G4int maxit = 0;
193   // Loop checking, 05-Aug-2015, Vladimir Ivan    193   // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
194   while (Da*Db > 0.0 && maxit < 1000)             194   while (Da*Db > 0.0 && maxit < 1000) 
195     {                                             195     {
196       ++maxit;                                    196       ++maxit;
197       Tb += 0.5*Tb;                               197       Tb += 0.5*Tb;   
198       Db = (U + FreeInternalE0 - GetPartitionE    198       Db = (U + FreeInternalE0 - GetPartitionEnergy(Tb))/U;
199     }                                             199     }
200                                                   200   
201   G4double eps = 1.0e-14*std::abs(Ta-Tb);         201   G4double eps = 1.0e-14*std::abs(Ta-Tb);
202                                                   202   
203   for (G4int i = 0; i < 1000; i++)                203   for (G4int i = 0; i < 1000; i++) 
204     {                                             204     {
205       Tmid = (Ta+Tb)/2.0;                         205       Tmid = (Ta+Tb)/2.0;
206       if (std::fabs(Ta-Tb) <= eps) return Tmid    206       if (std::fabs(Ta-Tb) <= eps) return Tmid;
207       G4double Dmid = (U + FreeInternalE0 - Ge    207       G4double Dmid = (U + FreeInternalE0 - GetPartitionEnergy(Tmid))/U;
208       if (std::fabs(Dmid) < 0.003) return Tmid    208       if (std::fabs(Dmid) < 0.003) return Tmid;
209       if (Da*Dmid < 0.0)                          209       if (Da*Dmid < 0.0) 
210         {                                         210         {
211           Tb = Tmid;                              211           Tb = Tmid;
212           Db = Dmid;                              212           Db = Dmid;
213         }                                         213         } 
214       else                                        214       else 
215         {                                         215         {
216           Ta = Tmid;                              216           Ta = Tmid;
217           Da = Dmid;                              217           Da = Dmid;
218         }                                         218         } 
219     }                                             219     }
220   // if we arrive here the temperature could n    220   // if we arrive here the temperature could not be calculated
221   G4cout << "G4StatMFMicroPartition::CalcParti    221   G4cout << "G4StatMFMicroPartition::CalcPartitionTemperature: I can't calculate the temperature"  
222          << G4endl;                               222          << G4endl;
223   // and set probability to 0 returning T < 0     223   // and set probability to 0 returning T < 0
224   return -1.0;                                    224   return -1.0;
225                                                   225   
226 }                                                 226 }
227                                                   227 
228 G4double G4StatMFMicroPartition::CalcPartition    228 G4double G4StatMFMicroPartition::CalcPartitionProbability(G4double U,
229                 G4double FreeInternalE0,          229                 G4double FreeInternalE0,
230                 G4double SCompound)               230                 G4double SCompound)
231 {                                                 231 { 
232   G4double T = CalcPartitionTemperature(U,Free    232   G4double T = CalcPartitionTemperature(U,FreeInternalE0);
233   if ( T <= 0.0) return _Probability = 0.0;       233   if ( T <= 0.0) return _Probability = 0.0;
234   _Temperature = T;                               234   _Temperature = T;
235                                                   235   
236   G4Pow* g4calc = G4Pow::GetInstance();           236   G4Pow* g4calc = G4Pow::GetInstance();
237                                                   237   
238   // Factorial of fragment multiplicity           238   // Factorial of fragment multiplicity
239   G4double Fact = 1.0;                            239   G4double Fact = 1.0;
240   unsigned int i;                                 240   unsigned int i;
241   for (i = 0; i < _thePartition.size() - 1; i+    241   for (i = 0; i < _thePartition.size() - 1; i++) 
242     {                                             242     {
243       G4double f = 1.0;                           243       G4double f = 1.0;
244       for (unsigned int ii = i+1; i< _theParti    244       for (unsigned int ii = i+1; i< _thePartition.size(); i++) 
245         {                                         245         {
246           if (_thePartition[i] == _thePartitio    246           if (_thePartition[i] == _thePartition[ii]) f++;
247         }                                         247         }
248       Fact *= f;                                  248       Fact *= f;
249   }                                               249   }
250                                                   250   
251   G4double ProbDegeneracy = 1.0;                  251   G4double ProbDegeneracy = 1.0;
252   G4double ProbA32 = 1.0;                         252   G4double ProbA32 = 1.0; 
253                                                   253   
254   for (i = 0; i < _thePartition.size(); i++)      254   for (i = 0; i < _thePartition.size(); i++) 
255     {                                             255     {
256       ProbDegeneracy *= GetDegeneracyFactor(_t    256       ProbDegeneracy *= GetDegeneracyFactor(_thePartition[i]);
257       ProbA32 *= _thePartition[i]*std::sqrt((G    257       ProbA32 *= _thePartition[i]*std::sqrt((G4double)_thePartition[i]);
258     }                                             258     }
259                                                   259   
260   // Compute entropy                              260   // Compute entropy
261   G4double PartitionEntropy = 0.0;                261   G4double PartitionEntropy = 0.0;
262   for (i = 0; i < _thePartition.size(); i++)      262   for (i = 0; i < _thePartition.size(); i++) 
263     {                                             263     {
264       // interaction entropy for alpha            264       // interaction entropy for alpha
265       if (_thePartition[i] == 4)                  265       if (_thePartition[i] == 4) 
266         {                                         266         {
267           PartitionEntropy +=                     267           PartitionEntropy += 
268             2.0*T*_thePartition[i]/InvLevelDen    268             2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i]);
269         }                                         269         }
270       // interaction entropy for Af > 4           270       // interaction entropy for Af > 4
271       else if (_thePartition[i] > 4)              271       else if (_thePartition[i] > 4) 
272         {                                         272         {
273           PartitionEntropy +=                     273           PartitionEntropy += 
274             2.0*T*_thePartition[i]/InvLevelDen    274             2.0*T*_thePartition[i]/InvLevelDensity(_thePartition[i])
275             - G4StatMFParameters::DBetaDT(T) *    275             - G4StatMFParameters::DBetaDT(T) * g4calc->Z23(_thePartition[i]);
276         }                                         276         } 
277     }                                             277     }
278                                                   278   
279   // Thermal Wave Lenght = std::sqrt(2 pi hbar    279   // Thermal Wave Lenght = std::sqrt(2 pi hbar^2 / nucleon_mass T)
280   G4double ThermalWaveLenght3 = 16.15*fermi/st    280   G4double ThermalWaveLenght3 = 16.15*fermi/std::sqrt(T);
281   ThermalWaveLenght3 = ThermalWaveLenght3*Ther    281   ThermalWaveLenght3 = ThermalWaveLenght3*ThermalWaveLenght3*ThermalWaveLenght3;
282                                                   282   
283   // Translational Entropy                        283   // Translational Entropy
284   G4double kappa = 1. + elm_coupling*(g4calc-> << 284   G4double kappa = 1. + elm_coupling*(g4calc->Z13(_thePartition.size())-1.0)
285                     /(G4StatMFParameters::Getr    285                     /(G4StatMFParameters::Getr0()*g4calc->Z13(theA));
286   kappa = kappa*kappa*kappa;                      286   kappa = kappa*kappa*kappa;
287   kappa -= 1.;                                    287   kappa -= 1.;
288   G4double V0 = (4./3.)*pi*theA*G4StatMFParame    288   G4double V0 = (4./3.)*pi*theA*G4StatMFParameters::Getr0()*G4StatMFParameters::Getr0()*
289     G4StatMFParameters::Getr0();                  289     G4StatMFParameters::Getr0();
290   G4double FreeVolume = kappa*V0;                 290   G4double FreeVolume = kappa*V0;
291   G4double TranslationalS = std::max(0.0, G4Lo    291   G4double TranslationalS = std::max(0.0, G4Log(ProbA32/Fact) +
292       (_thePartition.size()-1.0)*G4Log(FreeVol    292       (_thePartition.size()-1.0)*G4Log(FreeVolume/ThermalWaveLenght3) +
293       1.5*(_thePartition.size()-1.0) - 1.5*g4c    293       1.5*(_thePartition.size()-1.0) - 1.5*g4calc->logZ(theA));
294                                                   294   
295   PartitionEntropy += G4Log(ProbDegeneracy) +     295   PartitionEntropy += G4Log(ProbDegeneracy) + TranslationalS;
296   _Entropy = PartitionEntropy;                    296   _Entropy = PartitionEntropy;
297                                                   297   
298   // And finally compute probability of fragme    298   // And finally compute probability of fragment configuration
299   G4double exponent = PartitionEntropy-SCompou    299   G4double exponent = PartitionEntropy-SCompound;
300   if (exponent > 300.0) exponent = 300.0;         300   if (exponent > 300.0) exponent = 300.0;
301   return _Probability = G4Exp(exponent);          301   return _Probability = G4Exp(exponent);
302 }                                                 302 }
303                                                   303 
304 G4double G4StatMFMicroPartition::GetDegeneracy    304 G4double G4StatMFMicroPartition::GetDegeneracyFactor(G4int A)
305 {                                                 305 {
306   // Degeneracy factors are statistical factor    306   // Degeneracy factors are statistical factors
307   // DegeneracyFactor for nucleon is (2S_n + 1    307   // DegeneracyFactor for nucleon is (2S_n + 1)(2I_n + 1) = 4
308   G4double DegFactor = 0;                         308   G4double DegFactor = 0;
309   if (A > 4) DegFactor = 1.0;                     309   if (A > 4) DegFactor = 1.0;
310   else if (A == 1) DegFactor = 4.0;     // nuc    310   else if (A == 1) DegFactor = 4.0;     // nucleon
311   else if (A == 2) DegFactor = 3.0;     // Deu    311   else if (A == 2) DegFactor = 3.0;     // Deuteron
312   else if (A == 3) DegFactor = 4.0;     // Tri    312   else if (A == 3) DegFactor = 4.0;     // Triton + He3
313   else if (A == 4) DegFactor = 1.0;     // alp    313   else if (A == 4) DegFactor = 1.0;     // alpha
314   return DegFactor;                               314   return DegFactor;
315 }                                                 315 }
316                                                   316 
317 G4StatMFChannel * G4StatMFMicroPartition::Choo    317 G4StatMFChannel * G4StatMFMicroPartition::ChooseZ(G4int A0, G4int Z0, G4double MeanT)
318 // Gives fragments charges                        318 // Gives fragments charges
319 {                                                 319 {
320   std::vector<G4int> FragmentsZ;                  320   std::vector<G4int> FragmentsZ;
321                                                   321   
322   G4int ZBalance = 0;                             322   G4int ZBalance = 0;
323   do                                              323   do 
324     {                                             324     {
325       G4double CC = G4StatMFParameters::GetGam    325       G4double CC = G4StatMFParameters::GetGamma0()*8.0;
326       G4int SumZ = 0;                             326       G4int SumZ = 0;
327       for (unsigned int i = 0; i < _thePartiti    327       for (unsigned int i = 0; i < _thePartition.size(); i++) 
328         {                                         328         {
329           G4double ZMean;                         329           G4double ZMean;
330           G4double Af = _thePartition[i];         330           G4double Af = _thePartition[i];
331           if (Af > 1.5 && Af < 4.5) ZMean = 0.    331           if (Af > 1.5 && Af < 4.5) ZMean = 0.5*Af;
332           else ZMean = Af*Z0/A0;                  332           else ZMean = Af*Z0/A0;
333           G4double ZDispersion = std::sqrt(Af     333           G4double ZDispersion = std::sqrt(Af * MeanT/CC);
334           G4int Zf;                               334           G4int Zf;
335           do                                      335           do 
336             {                                     336             {
337               Zf = static_cast<G4int>(G4RandGa    337               Zf = static_cast<G4int>(G4RandGauss::shoot(ZMean,ZDispersion));
338       }                                           338       } 
339           // Loop checking, 05-Aug-2015, Vladi    339           // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
340           while (Zf < 0 || Zf > Af);              340           while (Zf < 0 || Zf > Af);
341           FragmentsZ.push_back(Zf);               341           FragmentsZ.push_back(Zf);
342           SumZ += Zf;                             342           SumZ += Zf;
343   }                                               343   }
344       ZBalance = Z0 - SumZ;                       344       ZBalance = Z0 - SumZ;
345     }                                             345     } 
346   // Loop checking, 05-Aug-2015, Vladimir Ivan    346   // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
347   while (std::abs(ZBalance) > 1);                 347   while (std::abs(ZBalance) > 1);
348   FragmentsZ[0] += ZBalance;                      348   FragmentsZ[0] += ZBalance;
349                                                   349   
350   G4StatMFChannel * theChannel = new G4StatMFC    350   G4StatMFChannel * theChannel = new G4StatMFChannel;
351   for (unsigned int i = 0; i < _thePartition.s    351   for (unsigned int i = 0; i < _thePartition.size(); i++)
352     {                                             352     {
353       theChannel->CreateFragment(_thePartition    353       theChannel->CreateFragment(_thePartition[i],FragmentsZ[i]);
354     }                                             354     }
355                                                   355 
356   return theChannel;                              356   return theChannel;
357 }                                                 357 }
358                                                   358