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 #include "G4StatMFMacroMultiplicity.hh"
00039 #include "G4PhysicalConstants.hh"
00040
00041
00042 G4StatMFMacroMultiplicity &
00043 G4StatMFMacroMultiplicity::operator=(const G4StatMFMacroMultiplicity & )
00044 {
00045 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator= meant to not be accessable");
00046 return *this;
00047 }
00048
00049 G4bool G4StatMFMacroMultiplicity::operator==(const G4StatMFMacroMultiplicity & ) const
00050 {
00051 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator== meant to not be accessable");
00052 return false;
00053 }
00054
00055
00056 G4bool G4StatMFMacroMultiplicity::operator!=(const G4StatMFMacroMultiplicity & ) const
00057 {
00058 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::operator!= meant to not be accessable");
00059 return true;
00060 }
00061
00062
00063
00064
00065 G4double G4StatMFMacroMultiplicity::CalcChemicalPotentialMu(void)
00066
00067
00068 {
00069 G4double CP = ((3./5.)*elm_coupling/G4StatMFParameters::Getr0())*
00070 (1.0-1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1.0/3.0));
00071
00072
00073
00074 G4double ZA5 = _theClusters->operator[](4)->GetZARatio();
00075 G4double ILD5 = _theClusters->operator[](4)->GetInvLevelDensity();
00076 _ChemPotentialMu = -G4StatMFParameters::GetE0()-
00077 _MeanTemperature*_MeanTemperature/ILD5 -
00078 _ChemPotentialNu*ZA5 +
00079 G4StatMFParameters::GetGamma0()*(1.0-2.0*ZA5)*(1.0-2.0*ZA5) +
00080 (2.0/3.0)*G4StatMFParameters::Beta(_MeanTemperature)/std::pow(5.,1./3.) +
00081 (5.0/3.0)*CP*ZA5*ZA5*std::pow(5.,2./3.) -
00082 1.5*_MeanTemperature/5.0;
00083
00084
00085
00086 G4double ChemPa = _ChemPotentialMu;
00087 if (ChemPa/_MeanTemperature > 10.0) ChemPa = 10.0*_MeanTemperature;
00088 G4double ChemPb = ChemPa - 0.5*std::abs(ChemPa);
00089
00090
00091 G4double fChemPa = this->operator()(ChemPa);
00092 G4double fChemPb = this->operator()(ChemPb);
00093
00094
00095
00096
00097 G4double intervalWidth = 1.e-4;
00098
00099
00100 G4int iterations = 0;
00101 while (fChemPa*fChemPb > 0.0 && iterations < 100)
00102 {
00103 if (std::abs(fChemPa) <= std::abs(fChemPb))
00104 {
00105 ChemPa += 0.6*(ChemPa-ChemPb);
00106 fChemPa = this->operator()(ChemPa);
00107 iterations++;
00108 }
00109 else
00110 {
00111 ChemPb += 0.6*(ChemPb-ChemPa);
00112 fChemPb = this->operator()(ChemPb);
00113 iterations++;
00114 }
00115 }
00116
00117 if (fChemPa*fChemPb > 0.0)
00118 {
00119 G4cerr <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
00120 G4cerr <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
00121 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't bracket the root.");
00122 }
00123 else if (fChemPa*fChemPb < 0.0 && std::abs(ChemPa-ChemPb) > intervalWidth)
00124 {
00125 G4Solver<G4StatMFMacroMultiplicity> * theSolver = new G4Solver<G4StatMFMacroMultiplicity>(100,intervalWidth);
00126 theSolver->SetIntervalLimits(ChemPa,ChemPb);
00127
00128 if (!theSolver->Brent(*this))
00129 {
00130 G4cerr <<"G4StatMFMacroMultiplicity:"<<" ChemPa="<<ChemPa<<" ChemPb="<<ChemPb<< G4endl;
00131 G4cerr <<"G4StatMFMacroMultiplicity:"<<" fChemPa="<<fChemPa<<" fChemPb="<<fChemPb<< G4endl;
00132 throw G4HadronicException(__FILE__, __LINE__, "G4StatMFMacroMultiplicity::CalcChemicalPotentialMu: I couldn't find the root.");
00133 }
00134 _ChemPotentialMu = theSolver->GetRoot();
00135 delete theSolver;
00136 }
00137 else
00138 {
00139 _ChemPotentialMu = ChemPa;
00140 }
00141
00142 return _ChemPotentialMu;
00143 }
00144
00145
00146
00147 G4double G4StatMFMacroMultiplicity::CalcMeanA(const G4double mu)
00148 {
00149 G4double r03 = G4StatMFParameters::Getr0(); r03 *= r03*r03;
00150 G4double V0 = (4.0/3.0)*pi*theA*r03;
00151
00152 G4double MeanA = 0.0;
00153
00154 _MeanMultiplicity = 0.0;
00155
00156
00157 G4int n = 1;
00158 for (std::vector<G4VStatMFMacroCluster*>::iterator i = _theClusters->begin();
00159 i != _theClusters->end(); ++i)
00160 {
00161 G4double multip = (*i)->CalcMeanMultiplicity(V0*_Kappa,mu,_ChemPotentialNu,_MeanTemperature);
00162 MeanA += multip*static_cast<G4double>(n++);
00163 _MeanMultiplicity += multip;
00164 }
00165
00166 return MeanA;
00167 }