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 "G4StatMFMacroCanonical.hh"
00039 #include "G4PhysicalConstants.hh"
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4Pow.hh"
00042
00043
00044 G4StatMFMacroCanonical::G4StatMFMacroCanonical(const G4Fragment & theFragment)
00045 {
00046
00047
00048 _theClusters.push_back(new G4StatMFMacroNucleon);
00049 _theClusters.push_back(new G4StatMFMacroBiNucleon);
00050 _theClusters.push_back(new G4StatMFMacroTriNucleon);
00051 _theClusters.push_back(new G4StatMFMacroTetraNucleon);
00052 for (G4int i = 4; i < theFragment.GetA_asInt(); i++)
00053 _theClusters.push_back(new G4StatMFMacroMultiNucleon(i+1));
00054
00055
00056 Initialize(theFragment);
00057
00058 }
00059
00060
00061 G4StatMFMacroCanonical::~G4StatMFMacroCanonical()
00062 {
00063
00064 if (!_theClusters.empty())
00065 {
00066 std::for_each(_theClusters.begin(),_theClusters.end(),DeleteFragment());
00067 }
00068 }
00069
00070
00071 void G4StatMFMacroCanonical::Initialize(const G4Fragment & theFragment)
00072 {
00073
00074 G4int A = theFragment.GetA_asInt();
00075 G4int Z = theFragment.GetZ_asInt();
00076 G4double x = 1.0 - 2.0*Z/G4double(A);
00077 G4Pow* g4pow = G4Pow::GetInstance();
00078
00079
00080 __FreeInternalE0 = A*( -G4StatMFParameters::GetE0() +
00081 G4StatMFParameters::GetGamma0()*x*x)
00082 +
00083 G4StatMFParameters::GetBeta0()*g4pow->Z23(A) +
00084 (3.0/5.0)*elm_coupling*Z*Z/(G4StatMFParameters::Getr0()*
00085 g4pow->Z13(A));
00086
00087 CalculateTemperature(theFragment);
00088
00089 return;
00090 }
00091
00092 void G4StatMFMacroCanonical::CalculateTemperature(const G4Fragment & theFragment)
00093 {
00094
00095 G4double U = theFragment.GetExcitationEnergy();
00096
00097 G4int A = theFragment.GetA_asInt();
00098 G4int Z = theFragment.GetZ_asInt();
00099
00100
00101 G4double FragMult = std::max((1.0+(2.31/MeV)*(U/A - 3.5*MeV))*A/100.0, 2.0);
00102
00103
00104
00105 _Kappa = (1.0+elm_coupling*(std::pow(FragMult,1./3.)-1)/
00106 (G4StatMFParameters::Getr0()*G4Pow::GetInstance()->Z13(A)));
00107 _Kappa = _Kappa*_Kappa*_Kappa - 1.0;
00108
00109
00110 G4StatMFMacroTemperature * theTemp = new
00111 G4StatMFMacroTemperature(A,Z,U,__FreeInternalE0,_Kappa,&_theClusters);
00112
00113 __MeanTemperature = theTemp->CalcTemperature();
00114 _ChemPotentialNu = theTemp->GetChemicalPotentialNu();
00115 _ChemPotentialMu = theTemp->GetChemicalPotentialMu();
00116 __MeanMultiplicity = theTemp->GetMeanMultiplicity();
00117 __MeanEntropy = theTemp->GetEntropy();
00118
00119 delete theTemp;
00120
00121 return;
00122 }
00123
00124
00125
00126
00127 G4StatMFChannel * G4StatMFMacroCanonical::ChooseAandZ(const G4Fragment &theFragment)
00128
00129 {
00130 G4int A = theFragment.GetA_asInt();
00131 G4int Z = theFragment.GetZ_asInt();
00132
00133 std::vector<G4int> ANumbers(A);
00134
00135 G4double Multiplicity = ChooseA(A,ANumbers);
00136
00137 std::vector<G4int> FragmentsA;
00138
00139 G4int i = 0;
00140 for (i = 0; i < A; i++)
00141 {
00142 for (G4int j = 0; j < ANumbers[i]; j++) FragmentsA.push_back(i+1);
00143 }
00144
00145
00146 G4int im = 0;
00147 for (G4int j = 0; j < Multiplicity; j++)
00148 {
00149 G4int FragmentsAMax = 0;
00150 im = j;
00151 for (i = j; i < Multiplicity; i++)
00152 {
00153 if (FragmentsA[i] <= FragmentsAMax) { continue; }
00154 else
00155 {
00156 im = i;
00157 FragmentsAMax = FragmentsA[im];
00158 }
00159 }
00160
00161 if (im != j)
00162 {
00163 FragmentsA[im] = FragmentsA[j];
00164 FragmentsA[j] = FragmentsAMax;
00165 }
00166 }
00167
00168 return ChooseZ(Z,FragmentsA);
00169 }
00170
00171
00172 G4double G4StatMFMacroCanonical::ChooseA(G4int A, std::vector<G4int> & ANumbers)
00173
00174 {
00175 G4double multiplicity = 0.0;
00176 G4int i;
00177
00178 std::vector<G4double> AcumMultiplicity;
00179 AcumMultiplicity.reserve(A);
00180
00181 AcumMultiplicity.push_back((*(_theClusters.begin()))->GetMeanMultiplicity());
00182 for (std::vector<G4VStatMFMacroCluster*>::iterator it = _theClusters.begin()+1;
00183 it != _theClusters.end(); ++it)
00184 {
00185 AcumMultiplicity.push_back((*it)->GetMeanMultiplicity()+AcumMultiplicity.back());
00186 }
00187
00188 G4int CheckA;
00189 do {
00190 CheckA = -1;
00191 G4int SumA = 0;
00192 G4int ThisOne = 0;
00193 multiplicity = 0.0;
00194 for (i = 0; i < A; i++) ANumbers[i] = 0;
00195 do {
00196 G4double RandNumber = G4UniformRand()*__MeanMultiplicity;
00197 for (i = 0; i < A; i++) {
00198 if (RandNumber < AcumMultiplicity[i]) {
00199 ThisOne = i;
00200 break;
00201 }
00202 }
00203 multiplicity++;
00204 ANumbers[ThisOne] = ANumbers[ThisOne]+1;
00205 SumA += ThisOne+1;
00206 CheckA = A - SumA;
00207
00208 } while (CheckA > 0);
00209
00210 } while (CheckA < 0 || std::abs(__MeanMultiplicity - multiplicity) > std::sqrt(__MeanMultiplicity) + 1./2.);
00211
00212 return multiplicity;
00213 }
00214
00215
00216 G4StatMFChannel * G4StatMFMacroCanonical::ChooseZ(G4int & Z,
00217 std::vector<G4int> & FragmentsA)
00218
00219 {
00220 G4Pow* g4pow = G4Pow::GetInstance();
00221 std::vector<G4int> FragmentsZ;
00222
00223 G4int DeltaZ = 0;
00224 G4double CP = (3./5.)*(elm_coupling/G4StatMFParameters::Getr0())*
00225 (1.0 - 1.0/std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.));
00226
00227 G4int multiplicity = FragmentsA.size();
00228
00229 do
00230 {
00231 FragmentsZ.clear();
00232 G4int SumZ = 0;
00233 for (G4int i = 0; i < multiplicity; i++)
00234 {
00235 G4int A = FragmentsA[i];
00236 if (A <= 1)
00237 {
00238 G4double RandNumber = G4UniformRand();
00239 if (RandNumber < (*_theClusters.begin())->GetZARatio())
00240 {
00241 FragmentsZ.push_back(1);
00242 SumZ += FragmentsZ[i];
00243 }
00244 else FragmentsZ.push_back(0);
00245 }
00246 else
00247 {
00248 G4double RandZ;
00249 G4double CC = 8.0*G4StatMFParameters::GetGamma0()+2.0*CP*g4pow->Z23(FragmentsA[i]);
00250 G4double ZMean;
00251 if (FragmentsA[i] > 1 && FragmentsA[i] < 5) { ZMean = 0.5*FragmentsA[i]; }
00252 else ZMean = FragmentsA[i]*(4.0*G4StatMFParameters::GetGamma0()+_ChemPotentialNu)/CC;
00253 G4double ZDispersion = std::sqrt(FragmentsA[i]*__MeanTemperature/CC);
00254 G4int z;
00255 do
00256 {
00257 RandZ = G4RandGauss::shoot(ZMean,ZDispersion);
00258 z = static_cast<G4int>(RandZ+0.5);
00259 } while (z < 0 || z > A);
00260 FragmentsZ.push_back(z);
00261 SumZ += z;
00262 }
00263 }
00264 DeltaZ = Z - SumZ;
00265 }
00266 while (std::abs(DeltaZ) > 1);
00267
00268
00269 G4int idx = 0;
00270 if (DeltaZ < 0.0)
00271 {
00272 while (FragmentsZ[idx] < 1) { ++idx; }
00273 }
00274 FragmentsZ[idx] += DeltaZ;
00275
00276 G4StatMFChannel * theChannel = new G4StatMFChannel;
00277 for (G4int i = multiplicity-1; i >= 0; i--)
00278 {
00279 theChannel->CreateFragment(FragmentsA[i],FragmentsZ[i]);
00280 }
00281
00282 return theChannel;
00283 }