#include <G4StatMFMacroTemperature.hh>

Public Member Functions
	G4StatMFMacroTemperature (const G4double anA, const G4double aZ, const G4double ExEnergy, const G4double FreeE0, const G4double kappa, std::vector< G4VStatMFMacroCluster * > *ClusterVector)

	~G4StatMFMacroTemperature ()

G4double	operator() (const G4double T)

G4double	GetMeanMultiplicity (void) const

G4double	GetChemicalPotentialMu (void) const

G4double	GetChemicalPotentialNu (void) const

G4double	GetTemperature (void) const

G4double	GetEntropy (void) const

G4double	CalcTemperature (void)

Detailed Description

Definition at line 42 of file G4StatMFMacroTemperature.hh.

Constructor & Destructor Documentation

G4StatMFMacroTemperature::G4StatMFMacroTemperature	(	const G4double	anA,
		const G4double	aZ,
		const G4double	ExEnergy,
		const G4double	FreeE0,
		const G4double	kappa,
		std::vector< G4VStatMFMacroCluster * > *	ClusterVector
	)

Definition at line 46 of file G4StatMFMacroTemperature.cc.

References G4StatMFMacroTemperature().

Referenced by G4StatMFMacroTemperature().

                                                      :
   theA(anA),
   theZ(aZ),
   _ExEnergy(ExEnergy),
   _FreeInternalE0(FreeE0),
   _Kappa(kappa),
   _MeanMultiplicity(0.0),
   _MeanTemperature(0.0),
   _ChemPotentialMu(0.0),
   _ChemPotentialNu(0.0),
   _MeanEntropy(0.0),
   _theClusters(ClusterVector) 
 {}

G4StatMFMacroTemperature::~G4StatMFMacroTemperature ( )

Definition at line 62 of file G4StatMFMacroTemperature.cc.

63 {}

Member Function Documentation

G4double G4StatMFMacroTemperature::CalcTemperature ( void )

Definition at line 66 of file G4StatMFMacroTemperature.cc.

References G4Solver< Function >::Brent(), G4Solver< Function >::Crenshaw(), G4cerr, G4cout, G4endl, G4Solver< Function >::GetRoot(), G4INCL::Math::max(), python.hepunit::MeV, operator()(), and G4Solver< Function >::SetIntervalLimits().

 {
     // Inital guess for the interval of the ensemble temperature values
     G4double Ta = 0.5; 
     G4double Tb = std::max(std::sqrt(_ExEnergy/(theA*0.12)),0.01*MeV);
     
     G4double fTa = this->operator()(Ta); 
     G4double fTb = this->operator()(Tb); 
 
     // Bracketing the solution
     // T should be greater than 0.
     // The interval is [Ta,Tb]
     // We start with a value for Ta = 0.5 MeV
     // it should be enough to have fTa > 0 If it isn't 
     // the case, we decrease Ta. But carefully, because 
     // fTa growes very fast when Ta is near 0 and we could have
     // an overflow.
 
     G4int iterations = 0;  
     while (fTa < 0.0 && ++iterations < 10) {
     Ta -= 0.5*Ta;
     fTa = this->operator()(Ta);
     }
     // Usually, fTb will be less than 0, but if it is not the case: 
     iterations = 0;  
     while (fTa*fTb > 0.0 && iterations++ < 10) {
     Tb += 2.*std::fabs(Tb-Ta);
     fTb = this->operator()(Tb);
     }
     
     if (fTa*fTb > 0.0) {
       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.");
     }
 
     G4Solver<G4StatMFMacroTemperature> * theSolver = new G4Solver<G4StatMFMacroTemperature>(100,1.e-4);
     theSolver->SetIntervalLimits(Ta,Tb);
     if (!theSolver->Crenshaw(*this)){ 
       G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
       G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
     }
     _MeanTemperature = theSolver->GetRoot();
     G4double FunctionValureAtRoot =  this->operator()(_MeanTemperature);
     delete  theSolver;
 
     // 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 (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.");
       }
     }
     //G4cout << "G4StatMFMacroTemperature::CalcTemperature: function = " << FunctionValureAtRoot 
     //     << " T(MeV)= " << _MeanTemperature << G4endl;
     return _MeanTemperature;
 }