Changeset 1315 for trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:

• trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc (modified) (7 diffs)

trunk/source/processes/hadronic/models/de_excitation/multifragmentation/src/G4StatMFMacroTemperature.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4StatMFMacroTemperature.cc,v 1.7 2008/11/19 14:33:31 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03 $
+// $Id: G4StatMFMacroTemperature.cc,v 1.8 2010/04/16 17:04:08 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
 //
 // Hadronic Process: Nuclear De-excitations
@@ -36 +36 @@
 //          Moscow, pshenich@fias.uni-frankfurt.de) make algorithm closer to
 //          original MF model
+// 16.04.10 V.Ivanchenko improved logic of solving equation for tempetature
+//          to protect code from rare unwanted exception; moved constructor 
+//          and destructor to source  
 
 #include "G4StatMFMacroTemperature.hh"
+
+G4StatMFMacroTemperature::G4StatMFMacroTemperature(const G4double anA, const G4double aZ, 
+  const G4double ExEnergy, const G4double FreeE0, const G4double kappa, 
+  std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
+  theA(anA),
+  theZ(aZ),
+  _ExEnergy(ExEnergy),
+  _FreeInternalE0(FreeE0),
+  _Kappa(kappa),
+  _MeanMultiplicity(0.0),
+  _MeanTemperature(0.0),
+  _ChemPotentialMu(0.0),
+  _ChemPotentialNu(0.0),
+  _theClusters(ClusterVector) 
+{}
+
+G4StatMFMacroTemperature::G4StatMFMacroTemperature() 
+{}
+        
+G4StatMFMacroTemperature::~G4StatMFMacroTemperature() 
+{}
 
 // operators definitions
@@ -82 +106 @@
 
     G4int iterations = 0;  
-    while (fTa < 0.0 && iterations++ < 10) {
+    while (fTa < 0.0 && ++iterations < 10) {
         Ta -= 0.5*Ta;
         fTa = this->operator()(Ta);
@@ -89 +113 @@
     iterations = 0;  
     while (fTa*fTb > 0.0 && iterations++ < 10) {
-        Tb += 2.*std::abs(Tb-Ta);
+        Tb += 2.*std::fabs(Tb-Ta);
         fTb = this->operator()(Tb);
     }
@@ -96 +120 @@
       G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
       G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
-        throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
+      throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
     }
 
@@ -102 +126 @@
     theSolver->SetIntervalLimits(Ta,Tb);
     if (!theSolver->Crenshaw(*this)){ 
-      G4cerr <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
-      G4cerr <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
+      G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
+      G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
     }
     _MeanTemperature = theSolver->GetRoot();
@@ -111 +135 @@
     // Verify if the root is found and it is indeed within the physical domain, 
     // say, between 1 and 50 MeV, otherwise try Brent method:
-    if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
-    G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
-    G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
-    theSolverBrent->SetIntervalLimits(Ta,Tb);
-    if (!theSolverBrent->Brent(*this)){
-      G4cerr <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
-      G4cerr <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 
+    if (std::fabs(FunctionValureAtRoot) > 5.e-2) {
+      if (_MeanTemperature < 1. || _MeanTemperature > 50.) {
+        G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot 
+               << " solution? = " << _MeanTemperature << " MeV " << G4endl;
+        G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
+        theSolverBrent->SetIntervalLimits(Ta,Tb);
+        if (!theSolverBrent->Brent(*this)){
+          G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
+          G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl; 
+          throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
+        }
+
+        _MeanTemperature = theSolverBrent->GetRoot();
+        FunctionValureAtRoot =  this->operator()(_MeanTemperature);
+        delete theSolverBrent;
+      }
+      if (std::abs(FunctionValureAtRoot) > 5.e-2) {
+        //if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
+        G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
         throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
-    }
-
-    _MeanTemperature = theSolverBrent->GetRoot();
-    FunctionValureAtRoot =  this->operator()(_MeanTemperature);
-    delete theSolverBrent;
-
-     if (_MeanTemperature < 1. || _MeanTemperature > 50. || std::abs(FunctionValureAtRoot) > 5.e-2) {
-    G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
-        throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
-     }
-    }
-
+      }
+    }
+    //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = " << FunctionValureAtRoot 
+    //     << " T(MeV)= " << _MeanTemperature << G4endl;
     return _MeanTemperature;
 }