00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "G4StatMFMacroTemperature.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045
00046 G4StatMFMacroTemperature::G4StatMFMacroTemperature(const G4double anA, const G4double aZ,
00047 const G4double ExEnergy, const G4double FreeE0, const G4double kappa,
00048 std::vector<G4VStatMFMacroCluster*> * ClusterVector) :
00049 theA(anA),
00050 theZ(aZ),
00051 _ExEnergy(ExEnergy),
00052 _FreeInternalE0(FreeE0),
00053 _Kappa(kappa),
00054 _MeanMultiplicity(0.0),
00055 _MeanTemperature(0.0),
00056 _ChemPotentialMu(0.0),
00057 _ChemPotentialNu(0.0),
00058 _MeanEntropy(0.0),
00059 _theClusters(ClusterVector)
00060 {}
00061
00062 G4StatMFMacroTemperature::~G4StatMFMacroTemperature()
00063 {}
00064
00065
00066 G4double G4StatMFMacroTemperature::CalcTemperature(void)
00067 {
00068
00069 G4double Ta = 0.5;
00070 G4double Tb = std::max(std::sqrt(_ExEnergy/(theA*0.12)),0.01*MeV);
00071
00072 G4double fTa = this->operator()(Ta);
00073 G4double fTb = this->operator()(Tb);
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084 G4int iterations = 0;
00085 while (fTa < 0.0 && ++iterations < 10) {
00086 Ta -= 0.5*Ta;
00087 fTa = this->operator()(Ta);
00088 }
00089
00090 iterations = 0;
00091 while (fTa*fTb > 0.0 && iterations++ < 10) {
00092 Tb += 2.*std::fabs(Tb-Ta);
00093 fTb = this->operator()(Tb);
00094 }
00095
00096 if (fTa*fTb > 0.0) {
00097 G4cerr <<"G4StatMFMacroTemperature:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
00098 G4cerr <<"G4StatMFMacroTemperature:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
00099 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't bracket the solution.");
00100 }
00101
00102 G4Solver<G4StatMFMacroTemperature> * theSolver = new G4Solver<G4StatMFMacroTemperature>(100,1.e-4);
00103 theSolver->SetIntervalLimits(Ta,Tb);
00104 if (!theSolver->Crenshaw(*this)){
00105 G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
00106 G4cout <<"G4StatMFMacroTemperature, Crenshaw method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
00107 }
00108 _MeanTemperature = theSolver->GetRoot();
00109 G4double FunctionValureAtRoot = this->operator()(_MeanTemperature);
00110 delete theSolver;
00111
00112
00113
00114 if (std::fabs(FunctionValureAtRoot) > 5.e-2) {
00115 if (_MeanTemperature < 1. || _MeanTemperature > 50.) {
00116 G4cout << "Crenshaw method failed; function = " << FunctionValureAtRoot
00117 << " solution? = " << _MeanTemperature << " MeV " << G4endl;
00118 G4Solver<G4StatMFMacroTemperature> * theSolverBrent = new G4Solver<G4StatMFMacroTemperature>(200,1.e-3);
00119 theSolverBrent->SetIntervalLimits(Ta,Tb);
00120 if (!theSolverBrent->Brent(*this)){
00121 G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" Ta="<<Ta<<" Tb="<<Tb<< G4endl;
00122 G4cout <<"G4StatMFMacroTemperature, Brent method failed:"<<" fTa="<<fTa<<" fTb="<<fTb<< G4endl;
00123 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
00124 }
00125
00126 _MeanTemperature = theSolverBrent->GetRoot();
00127 FunctionValureAtRoot = this->operator()(_MeanTemperature);
00128 delete theSolverBrent;
00129 }
00130 if (std::abs(FunctionValureAtRoot) > 5.e-2) {
00131
00132 G4cout << "Brent method failed; function = " << FunctionValureAtRoot << " solution? = " << _MeanTemperature << " MeV " << G4endl;
00133 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroTemperature::CalcTemperature: I couldn't find the root with any method.");
00134 }
00135 }
00136
00137
00138 return _MeanTemperature;
00139 }
00140
00141
00142
00143 G4double G4StatMFMacroTemperature::FragsExcitEnergy(const G4double T)
00144
00145 {
00146
00147
00148 G4double R0 = G4StatMFParameters::Getr0()*std::pow(theA,1./3.);
00149 G4double R = R0*std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(), 1./3.);
00150 G4double FreeVol = _Kappa*(4.*pi/3.)*R0*R0*R0;
00151
00152
00153
00154 CalcChemicalPotentialNu(T);
00155
00156
00157
00158 G4double AverageEnergy = 0.0;
00159 std::vector<G4VStatMFMacroCluster*>::iterator i;
00160 for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
00161 {
00162 AverageEnergy += (*i)->GetMeanMultiplicity() * (*i)->CalcEnergy(T);
00163 }
00164
00165
00166 AverageEnergy += (3./5.)*elm_coupling*theZ*theZ/R;
00167
00168
00169 _MeanEntropy = 0.0;
00170 for (i = _theClusters->begin(); i != _theClusters->end(); ++i)
00171 {
00172 _MeanEntropy += (*i)->CalcEntropy(T,FreeVol);
00173 }
00174
00175
00176 return AverageEnergy - _FreeInternalE0;
00177
00178 }
00179
00180
00181 void G4StatMFMacroTemperature::CalcChemicalPotentialNu(const G4double T)
00182
00183
00184 {
00185 G4StatMFMacroChemicalPotential * theChemPot = new
00186 G4StatMFMacroChemicalPotential(theA,theZ,_Kappa,T,_theClusters);
00187
00188
00189 _ChemPotentialNu = theChemPot->CalcChemicalPotentialNu();
00190 _ChemPotentialMu = theChemPot->GetChemicalPotentialMu();
00191 _MeanMultiplicity = theChemPot->GetMeanMultiplicity();
00192
00193 delete theChemPot;
00194
00195 return;
00196
00197 }
00198
00199