Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMF.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/G4StatMF.cc (Version 11.3.0) and /processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMF.cc (Version 9.1)


  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 // $Id: G4StatMF.cc,v 1.5 2006/06/29 20:24:43 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-01 $
 27 //                                                 29 //
 28 // Hadronic Process: Nuclear De-excitations        30 // Hadronic Process: Nuclear De-excitations
 29 // by V. Lara                                      31 // by V. Lara
 30                                                    32 
 31 #include "G4StatMF.hh"                             33 #include "G4StatMF.hh"
 32 #include "G4PhysicalConstants.hh"              << 
 33 #include "G4SystemOfUnits.hh"                  << 
 34 #include "G4Pow.hh"                            << 
 35 #include "G4PhysicsModelCatalog.hh"            << 
 36                                                    34 
 37 G4StatMF::G4StatMF()                           <<  35 
                                                   >>  36 
                                                   >>  37 // Default constructor
                                                   >>  38 G4StatMF::G4StatMF() : _theEnsemble(0) {}
                                                   >>  39 
                                                   >>  40 
                                                   >>  41 // Destructor
                                                   >>  42 G4StatMF::~G4StatMF() {} //{if (_theEnsemble != 0) delete _theEnsemble;}
                                                   >>  43 
                                                   >>  44 
                                                   >>  45 // Copy constructor
                                                   >>  46 G4StatMF::G4StatMF(const G4StatMF & ) : G4VMultiFragmentation()
 38 {                                                  47 {
 39   _secID = G4PhysicsModelCatalog::GetModelID(" <<  48     throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::copy_constructor meant to not be accessable");
 40 }                                                  49 }
 41                                                    50 
 42 G4StatMF::~G4StatMF() {}                       << 
 43                                                    51 
 44 G4FragmentVector* G4StatMF::BreakItUp(const G4 <<  52 // Operators
                                                   >>  53 
                                                   >>  54 G4StatMF & G4StatMF::operator=(const G4StatMF & )
                                                   >>  55 {
                                                   >>  56     throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator= meant to not be accessable");
                                                   >>  57     return *this;
                                                   >>  58 }
                                                   >>  59 
                                                   >>  60 
                                                   >>  61 G4bool G4StatMF::operator==(const G4StatMF & )
                                                   >>  62 {
                                                   >>  63     throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator== meant to not be accessable");
                                                   >>  64     return false;
                                                   >>  65 }
                                                   >>  66  
                                                   >>  67 
                                                   >>  68 G4bool G4StatMF::operator!=(const G4StatMF & )
 45 {                                                  69 {
                                                   >>  70     throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::operator!= meant to not be accessable");
                                                   >>  71     return true;
                                                   >>  72 }
                                                   >>  73 
                                                   >>  74 
                                                   >>  75 
                                                   >>  76 
                                                   >>  77 
                                                   >>  78 G4FragmentVector * G4StatMF::BreakItUp(const G4Fragment &theFragment)
                                                   >>  79 {
                                                   >>  80   //  G4FragmentVector * theResult = new G4FragmentVector;
                                                   >>  81 
 46   if (theFragment.GetExcitationEnergy() <= 0.0     82   if (theFragment.GetExcitationEnergy() <= 0.0) {
 47     return nullptr;                            <<  83     G4FragmentVector * theResult = new G4FragmentVector;
                                                   >>  84     theResult->push_back(new G4Fragment(theFragment));
                                                   >>  85     return 0;
 48   }                                                86   }
 49                                                    87 
                                                   >>  88 
 50   // Maximun average multiplicity: M_0 = 2.6 f     89   // Maximun average multiplicity: M_0 = 2.6 for A ~ 200 
 51   // and M_0 = 3.3 for A <= 110                    90   // and M_0 = 3.3 for A <= 110
 52   G4double MaxAverageMultiplicity =                91   G4double MaxAverageMultiplicity = 
 53     G4StatMFParameters::GetMaxAverageMultiplic <<  92     G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA()));
 54                                                    93 
 55                                                    94   
 56     // We'll use two kinds of ensembles            95     // We'll use two kinds of ensembles
 57   G4StatMFMicroCanonical * theMicrocanonicalEn     96   G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
 58   G4StatMFMacroCanonical * theMacrocanonicalEn     97   G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
                                                   >>  98 
 59                                                    99   
 60   //------------------------------------------ << 100   //-------------------------------------------------------
 61   // Direct simulation part (Microcanonical en << 101   // Direct simulation part (Microcanonical ensemble)
 62   //------------------------------------------ << 102   //-------------------------------------------------------
 63                                                   103   
 64   // Microcanonical ensemble initialization    << 104   // Microcanonical ensemble initialization 
 65   theMicrocanonicalEnsemble = new G4StatMFMicr    105   theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
 66                                                   106 
 67   G4int Iterations = 0;                           107   G4int Iterations = 0;
 68   G4int IterationsLimit = 100000;              << 
 69   G4double Temperature = 0.0;                     108   G4double Temperature = 0.0;
 70                                                   109   
 71   G4bool FirstTime = true;                        110   G4bool FirstTime = true;
 72   G4StatMFChannel * theChannel = 0;               111   G4StatMFChannel * theChannel = 0;
 73                                                << 112 
 74   G4bool ChannelOk;                               113   G4bool ChannelOk;
 75   do {  // Try to de-excite as much as Iterati << 114   do {  // Try to de-excite as much as 10 times
 76     do {                                          115     do {
 77                                                   116       
 78       G4double theMeanMult = theMicrocanonical    117       G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
 79       if (theMeanMult <= MaxAverageMultiplicit    118       if (theMeanMult <= MaxAverageMultiplicity) {
 80   // G4cout << "MICROCANONICAL" << G4endl;        119   // G4cout << "MICROCANONICAL" << G4endl;
 81   // Choose fragments atomic numbers and charg    120   // Choose fragments atomic numbers and charges from direct simulation
 82   theChannel = theMicrocanonicalEnsemble->Choo    121   theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
 83   _theEnsemble = theMicrocanonicalEnsemble;       122   _theEnsemble = theMicrocanonicalEnsemble;
 84       } else {                                    123       } else {
 85   //------------------------------------------    124   //-----------------------------------------------------
 86   // Non direct simulation part (Macrocanonica    125   // Non direct simulation part (Macrocanonical Ensemble)
 87   //------------------------------------------    126   //-----------------------------------------------------
 88   if (FirstTime) {                                127   if (FirstTime) {
 89     // Macrocanonical ensemble initialization     128     // Macrocanonical ensemble initialization 
 90     theMacrocanonicalEnsemble = new G4StatMFMa    129     theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
 91     _theEnsemble = theMacrocanonicalEnsemble;     130     _theEnsemble = theMacrocanonicalEnsemble;
 92     FirstTime = false;                            131     FirstTime = false;
 93   }                                               132   }
 94   // G4cout << "MACROCANONICAL" << G4endl;        133   // G4cout << "MACROCANONICAL" << G4endl;
 95   // Select calculated fragment total multipli    134   // Select calculated fragment total multiplicity, 
 96   // fragment atomic numbers and fragment char    135   // fragment atomic numbers and fragment charges.
 97   theChannel = theMacrocanonicalEnsemble->Choo    136   theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
 98       }                                           137       }
 99                                                   138       
100       ChannelOk = theChannel->CheckFragments() << 139       if (!(ChannelOk = theChannel->CheckFragments())) delete theChannel; 
101       if (!ChannelOk) delete theChannel;       << 
102                                                   140       
103       // Loop checking, 05-Aug-2015, Vladimir  << 
104     } while (!ChannelOk);                         141     } while (!ChannelOk);
105                                                   142     
106                                                   143     
107     if (theChannel->GetMultiplicity() <= 1) {     144     if (theChannel->GetMultiplicity() <= 1) {
108       G4FragmentVector * theResult = new G4Fra    145       G4FragmentVector * theResult = new G4FragmentVector;
109       theResult->push_back(new G4Fragment(theF    146       theResult->push_back(new G4Fragment(theFragment));
110       delete theMicrocanonicalEnsemble;           147       delete theMicrocanonicalEnsemble;
111       if (theMacrocanonicalEnsemble != 0) dele    148       if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
112       delete theChannel;                          149       delete theChannel;
113       return theResult;                           150       return theResult;
114     }                                             151     }
115                                                   152     
116     //--------------------------------------      153     //--------------------------------------
117     // Second part of simulation procedure.       154     // Second part of simulation procedure.
118     //--------------------------------------      155     //--------------------------------------
119                                                   156     
120     // Find temperature of breaking channel.      157     // Find temperature of breaking channel.
121     Temperature = _theEnsemble->GetMeanTempera << 158     Temperature = _theEnsemble->GetMeanTemperature(); // Initial value for Temperature 
122                                                << 159     
123     if (FindTemperatureOfBreakingChannel(theFr    160     if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
124                                                << 161     
125     // Do not forget to delete this unusable c << 162   } while (Iterations++ < 10);
126     // otherwise for very proton-reach nuclei  << 163   
127     // number of iterations. N.B. "theChannel" << 164   
128                                                << 165   // If Iterations >= 10 means that we couldn't solve for temperature
129     // G4cout << " Iteration # " << Iterations << 166   if (Iterations >= 10) 
130                                                << 
131     delete theChannel;                         << 
132                                                << 
133     // Loop checking, 05-Aug-2015, Vladimir Iv << 
134   } while (Iterations++ < IterationsLimit );   << 
135                                                << 
136   // If Iterations >= IterationsLimit means th << 
137   if (Iterations >= IterationsLimit)           << 
138     throw G4HadronicException(__FILE__, __LINE    167     throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
139                                                << 168   
                                                   >> 169   
140   G4FragmentVector * theResult = theChannel->     170   G4FragmentVector * theResult = theChannel->
141     GetFragments(theFragment.GetA_asInt(),theF << 171     GetFragments(theFragment.GetA(),theFragment.GetZ(),Temperature);
142                                                << 172   
                                                   >> 173   
                                                   >> 174   
143   // ~~~~~~ Energy conservation Patch !!!!!!!!    175   // ~~~~~~ Energy conservation Patch !!!!!!!!!!!!!!!!!!!!!!
144   // Original nucleus 4-momentum in CM system     176   // Original nucleus 4-momentum in CM system
145   G4LorentzVector InitialMomentum(theFragment.    177   G4LorentzVector InitialMomentum(theFragment.GetMomentum());
146   InitialMomentum.boost(-InitialMomentum.boost    178   InitialMomentum.boost(-InitialMomentum.boostVector());
147   G4double ScaleFactor = 0.0;                     179   G4double ScaleFactor = 0.0;
148   G4double SavedScaleFactor = 0.0;                180   G4double SavedScaleFactor = 0.0;
149   do {                                            181   do {
150     G4double FragmentsEnergy = 0.0;               182     G4double FragmentsEnergy = 0.0;
151     for (auto const & ptr : *theResult) {      << 183     G4FragmentVector::iterator j;
152       FragmentsEnergy += ptr->GetMomentum().e( << 184     for (j = theResult->begin(); j != theResult->end(); j++) 
153     }                                          << 185       FragmentsEnergy += (*j)->GetMomentum().e();
154     if (0.0 == FragmentsEnergy) { break; }     << 
155     SavedScaleFactor = ScaleFactor;               186     SavedScaleFactor = ScaleFactor;
156     ScaleFactor = InitialMomentum.e()/Fragment    187     ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
157     G4ThreeVector ScaledMomentum(0.0,0.0,0.0);    188     G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
158     for (auto const & ptr : *theResult) {      << 189     for (j = theResult->begin(); j != theResult->end(); j++) {
159       ScaledMomentum = ScaleFactor * ptr->GetM << 190       ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
160       G4double Mass = ptr->GetMomentum().mag() << 191       G4double Mass = (*j)->GetMomentum().m();
161       G4LorentzVector NewMomentum;                192       G4LorentzVector NewMomentum;
162       NewMomentum.setVect(ScaledMomentum);        193       NewMomentum.setVect(ScaledMomentum);
163       NewMomentum.setE(std::sqrt(ScaledMomentu    194       NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
164       ptr->SetMomentum(NewMomentum);           << 195       (*j)->SetMomentum(NewMomentum);   
165     }                                             196     }
166     // Loop checking, 05-Aug-2015, Vladimir Iv << 
167   } while (ScaleFactor > 1.0+1.e-5 && std::abs    197   } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
168   // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!    198   // ~~~~~~ End of patch !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
169                                                   199   
170   // Perform Lorentz boost                        200   // Perform Lorentz boost
171   G4FragmentVector::iterator i;                   201   G4FragmentVector::iterator i;
172   for (i = theResult->begin(); i != theResult-    202   for (i = theResult->begin(); i != theResult->end(); i++) {
173     G4LorentzVector FourMom = (*i)->GetMomentu    203     G4LorentzVector FourMom = (*i)->GetMomentum();
174     FourMom.boost(theFragment.GetMomentum().bo    204     FourMom.boost(theFragment.GetMomentum().boostVector());
175     (*i)->SetMomentum(FourMom);                   205     (*i)->SetMomentum(FourMom);
176     (*i)->SetCreatorModelID(_secID);           << 206 #ifdef PRECOMPOUND_TEST
                                                   >> 207     (*i)->SetCreatorModel(G4String("G4StatMF"));
                                                   >> 208 #endif
177   }                                               209   }
178                                                   210   
179   // garbage collection                           211   // garbage collection
180   delete theMicrocanonicalEnsemble;               212   delete theMicrocanonicalEnsemble;
181   if (theMacrocanonicalEnsemble != 0) delete t    213   if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
182   delete theChannel;                              214   delete theChannel;
183                                                   215   
184   return theResult;                               216   return theResult;
185 }                                                 217 }
186                                                   218 
187                                                   219 
188 G4bool G4StatMF::FindTemperatureOfBreakingChan    220 G4bool G4StatMF::FindTemperatureOfBreakingChannel(const G4Fragment & theFragment,
189               const G4StatMFChannel * aChannel    221               const G4StatMFChannel * aChannel,
190               G4double & Temperature)             222               G4double & Temperature)
191   // This finds temperature of breaking channe    223   // This finds temperature of breaking channel.
192 {                                                 224 {
193   G4int A = theFragment.GetA_asInt();          << 225   G4double A = theFragment.GetA();
194   G4int Z = theFragment.GetZ_asInt();          << 226   G4double Z = theFragment.GetZ();
195   G4double U = theFragment.GetExcitationEnergy    227   G4double U = theFragment.GetExcitationEnergy();
196                                                   228   
197   G4double T = std::max(Temperature,0.0012*MeV << 229   G4double T = std::max(Temperature,0.0012*MeV);
                                                   >> 230   
198   G4double Ta = T;                                231   G4double Ta = T;
                                                   >> 232   G4double Tb = T;
                                                   >> 233   
                                                   >> 234   
199   G4double TotalEnergy = CalcEnergy(A,Z,aChann    235   G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
200                                                   236   
201   G4double Da = (U - TotalEnergy)/U;              237   G4double Da = (U - TotalEnergy)/U;
202   G4double Db = 0.0;                              238   G4double Db = 0.0;
203                                                   239   
204   // bracketing the solution                      240   // bracketing the solution
205   if (Da == 0.0) {                                241   if (Da == 0.0) {
206     Temperature = T;                              242     Temperature = T;
207     return true;                                  243     return true;
208   } else if (Da < 0.0) {                          244   } else if (Da < 0.0) {
209     do {                                          245     do {
210       T *= 0.5;                                << 246       Tb -= 0.5 * std::abs(Tb);
211       if (T < 0.001*MeV) return false;         << 247       T = Tb;
                                                   >> 248       if (Tb < 0.001*MeV) return false;
212                                                   249       
213       TotalEnergy = CalcEnergy(A,Z,aChannel,T)    250       TotalEnergy = CalcEnergy(A,Z,aChannel,T);
214                                                   251       
215       Db = (U - TotalEnergy)/U;                   252       Db = (U - TotalEnergy)/U;
216       // Loop checking, 05-Aug-2015, Vladimir  << 
217     } while (Db < 0.0);                           253     } while (Db < 0.0);
218                                                   254     
219   } else {                                        255   } else {
220     do {                                          256     do {
221       T *= 1.5;                                << 257       Tb += 0.5 * std::abs(Tb);
                                                   >> 258       T = Tb;
222                                                   259       
223       TotalEnergy = CalcEnergy(A,Z,aChannel,T)    260       TotalEnergy = CalcEnergy(A,Z,aChannel,T);
224                                                   261       
225       Db = (U - TotalEnergy)/U;                   262       Db = (U - TotalEnergy)/U;
226       // Loop checking, 05-Aug-2015, Vladimir  << 
227     } while (Db > 0.0);                           263     } while (Db > 0.0); 
228   }                                               264   }
229                                                   265   
230   G4double eps = 1.0e-14 * std::abs(T-Ta);     << 266   G4double eps = 1.0e-14 * std::abs(Tb-Ta);
231   //G4double eps = 1.0e-3 ;                       267   //G4double eps = 1.0e-3 ;
232                                                   268   
233   // Start the bisection method                   269   // Start the bisection method
234   for (G4int j = 0; j < 1000; j++) {              270   for (G4int j = 0; j < 1000; j++) {
235     G4double Tc =  (Ta+T)*0.5;                 << 271     G4double Tc =  (Ta+Tb)/2.0;
236     if (std::abs(Ta-Tc) <= eps) {                 272     if (std::abs(Ta-Tc) <= eps) {
237       Temperature = Tc;                           273       Temperature = Tc;
238       return true;                                274       return true;
239     }                                             275     }
240                                                   276     
241     T = Tc;                                    << 277     T = Tc;
                                                   >> 278     
242     TotalEnergy = CalcEnergy(A,Z,aChannel,T);     279     TotalEnergy = CalcEnergy(A,Z,aChannel,T);
                                                   >> 280     
243     G4double Dc = (U - TotalEnergy)/U;            281     G4double Dc = (U - TotalEnergy)/U; 
244                                                   282     
245     if (Dc == 0.0) {                              283     if (Dc == 0.0) {
246       Temperature  = Tc;                          284       Temperature  = Tc;
247       return true;                                285       return true;
248     }                                             286     }
                                                   >> 287     
249     if (Da*Dc < 0.0) {                            288     if (Da*Dc < 0.0) {
250       T  = Tc;                                 << 289       Tb = Tc;
251       Db = Dc;                                    290       Db = Dc;
252     } else {                                      291     } else {
253       Ta = Tc;                                    292       Ta = Tc;
254       Da = Dc;                                    293       Da = Dc;
255     }                                             294     }
256   }                                               295   }
257                                                   296   
258   Temperature  = (Ta+T)*0.5;                   << 297   Temperature  = (Ta+Tb)/2.0;
259   return false;                                   298   return false;
260 }                                                 299 }
261                                                   300 
262 G4double G4StatMF::CalcEnergy(G4int A, G4int Z << 301 
263             G4double T)                        << 302 
                                                   >> 303 G4double G4StatMF::CalcEnergy(const G4double A, const G4double Z, const G4StatMFChannel * aChannel,
                                                   >> 304             const G4double T)
264 {                                                 305 {
265   G4double MassExcess0 = G4NucleiProperties::G << 306     G4double MassExcess0 = G4NucleiProperties::GetMassExcess(static_cast<G4int>(A),static_cast<G4int>(Z));
266   G4double ChannelEnergy = aChannel->GetFragme << 307   
267   return -MassExcess0 + G4StatMFParameters::Ge << 308     G4double Coulomb = (3./5.)*(elm_coupling*Z*Z)*std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
                                                   >> 309   (G4StatMFParameters::Getr0()*std::pow(A,1./3.));
                                                   >> 310 
                                                   >> 311     G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T);
                                                   >> 312   
                                                   >> 313     return -MassExcess0 + Coulomb + ChannelEnergy;
                                                   >> 314 
268 }                                                 315 }
269                                                   316 
270                                                   317 
271                                                   318 
272                                                   319