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 #include "G4StatMF.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4Pow.hh"
00036
00037
00038
00039 G4StatMF::G4StatMF() : _theEnsemble(0) {}
00040
00041
00042
00043 G4StatMF::~G4StatMF() {}
00044
00045
00046 G4FragmentVector * G4StatMF::BreakItUp(const G4Fragment &theFragment)
00047 {
00048
00049
00050 if (theFragment.GetExcitationEnergy() <= 0.0) {
00051
00052
00053 return 0;
00054 }
00055
00056
00057
00058
00059 G4double MaxAverageMultiplicity =
00060 G4StatMFParameters::GetMaxAverageMultiplicity(static_cast<G4int>(theFragment.GetA()));
00061
00062
00063
00064 G4StatMFMicroCanonical * theMicrocanonicalEnsemble = 0;
00065 G4StatMFMacroCanonical * theMacrocanonicalEnsemble = 0;
00066
00067
00068
00069
00070
00071
00072
00073 theMicrocanonicalEnsemble = new G4StatMFMicroCanonical(theFragment);
00074
00075 G4int Iterations = 0;
00076 G4int IterationsLimit = 100000;
00077 G4double Temperature = 0.0;
00078
00079 G4bool FirstTime = true;
00080 G4StatMFChannel * theChannel = 0;
00081
00082 G4bool ChannelOk;
00083 do {
00084 do {
00085
00086 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
00087 if (theMeanMult <= MaxAverageMultiplicity) {
00088
00089
00090 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
00091 _theEnsemble = theMicrocanonicalEnsemble;
00092 } else {
00093
00094
00095
00096 if (FirstTime) {
00097
00098 theMacrocanonicalEnsemble = new G4StatMFMacroCanonical(theFragment);
00099 _theEnsemble = theMacrocanonicalEnsemble;
00100 FirstTime = false;
00101 }
00102
00103
00104
00105 theChannel = theMacrocanonicalEnsemble->ChooseAandZ(theFragment);
00106 }
00107
00108 ChannelOk = theChannel->CheckFragments();
00109 if (!ChannelOk) delete theChannel;
00110
00111 } while (!ChannelOk);
00112
00113
00114 if (theChannel->GetMultiplicity() <= 1) {
00115 G4FragmentVector * theResult = new G4FragmentVector;
00116 theResult->push_back(new G4Fragment(theFragment));
00117 delete theMicrocanonicalEnsemble;
00118 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
00119 delete theChannel;
00120 return theResult;
00121 }
00122
00123
00124
00125
00126
00127
00128 Temperature = _theEnsemble->GetMeanTemperature();
00129
00130 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature)) break;
00131
00132
00133
00134
00135
00136
00137
00138 delete theChannel;
00139
00140 } while (Iterations++ < IterationsLimit );
00141
00142
00143
00144
00145 if (Iterations >= IterationsLimit)
00146 throw G4HadronicException(__FILE__, __LINE__, "G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
00147
00148
00149 G4FragmentVector * theResult = theChannel->
00150 GetFragments(theFragment.GetA_asInt(),theFragment.GetZ_asInt(),Temperature);
00151
00152
00153
00154
00155
00156 G4LorentzVector InitialMomentum(theFragment.GetMomentum());
00157 InitialMomentum.boost(-InitialMomentum.boostVector());
00158 G4double ScaleFactor = 0.0;
00159 G4double SavedScaleFactor = 0.0;
00160 do {
00161 G4double FragmentsEnergy = 0.0;
00162 G4FragmentVector::iterator j;
00163 for (j = theResult->begin(); j != theResult->end(); j++)
00164 FragmentsEnergy += (*j)->GetMomentum().e();
00165 SavedScaleFactor = ScaleFactor;
00166 ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
00167 G4ThreeVector ScaledMomentum(0.0,0.0,0.0);
00168 for (j = theResult->begin(); j != theResult->end(); j++) {
00169 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
00170 G4double Mass = (*j)->GetMomentum().m();
00171 G4LorentzVector NewMomentum;
00172 NewMomentum.setVect(ScaledMomentum);
00173 NewMomentum.setE(std::sqrt(ScaledMomentum.mag2()+Mass*Mass));
00174 (*j)->SetMomentum(NewMomentum);
00175 }
00176 } while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
00177
00178
00179
00180 G4FragmentVector::iterator i;
00181 for (i = theResult->begin(); i != theResult->end(); i++) {
00182 G4LorentzVector FourMom = (*i)->GetMomentum();
00183 FourMom.boost(theFragment.GetMomentum().boostVector());
00184 (*i)->SetMomentum(FourMom);
00185 #ifdef PRECOMPOUND_TEST
00186 (*i)->SetCreatorModel(G4String("G4StatMF"));
00187 #endif
00188 }
00189
00190
00191 delete theMicrocanonicalEnsemble;
00192 if (theMacrocanonicalEnsemble != 0) delete theMacrocanonicalEnsemble;
00193 delete theChannel;
00194
00195 return theResult;
00196 }
00197
00198
00199 G4bool G4StatMF::FindTemperatureOfBreakingChannel(const G4Fragment & theFragment,
00200 const G4StatMFChannel * aChannel,
00201 G4double & Temperature)
00202
00203 {
00204 G4int A = theFragment.GetA_asInt();
00205 G4int Z = theFragment.GetZ_asInt();
00206 G4double U = theFragment.GetExcitationEnergy();
00207
00208 G4double T = std::max(Temperature,0.0012*MeV);
00209
00210 G4double Ta = T;
00211 G4double Tb = T;
00212
00213
00214 G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
00215
00216 G4double Da = (U - TotalEnergy)/U;
00217 G4double Db = 0.0;
00218
00219
00220 if (Da == 0.0) {
00221 Temperature = T;
00222 return true;
00223 } else if (Da < 0.0) {
00224 do {
00225 Tb -= 0.5 * std::fabs(Tb);
00226 T = Tb;
00227 if (Tb < 0.001*MeV) return false;
00228
00229 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
00230
00231 Db = (U - TotalEnergy)/U;
00232 } while (Db < 0.0);
00233
00234 } else {
00235 do {
00236 Tb += 0.5 * std::fabs(Tb);
00237 T = Tb;
00238
00239 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
00240
00241 Db = (U - TotalEnergy)/U;
00242 } while (Db > 0.0);
00243 }
00244
00245 G4double eps = 1.0e-14 * std::abs(Tb-Ta);
00246
00247
00248
00249 for (G4int j = 0; j < 1000; j++) {
00250 G4double Tc = (Ta+Tb)/2.0;
00251 if (std::abs(Ta-Tc) <= eps) {
00252 Temperature = Tc;
00253 return true;
00254 }
00255
00256 T = Tc;
00257
00258 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
00259
00260 G4double Dc = (U - TotalEnergy)/U;
00261
00262 if (Dc == 0.0) {
00263 Temperature = Tc;
00264 return true;
00265 }
00266
00267 if (Da*Dc < 0.0) {
00268 Tb = Tc;
00269 Db = Dc;
00270 } else {
00271 Ta = Tc;
00272 Da = Dc;
00273 }
00274 }
00275
00276 Temperature = (Ta+Tb)/2.0;
00277 return false;
00278 }
00279
00280 G4double G4StatMF::CalcEnergy(G4int A, G4int Z, const G4StatMFChannel * aChannel,
00281 G4double T)
00282 {
00283 G4double MassExcess0 = G4NucleiProperties::GetMassExcess(A,Z);
00284
00285 G4double Coulomb = (3./5.)*(elm_coupling*Z*Z)
00286 *std::pow(1.0+G4StatMFParameters::GetKappaCoulomb(),1./3.)/
00287 (G4StatMFParameters::Getr0()*G4Pow::GetInstance()->Z13(A));
00288
00289 G4double ChannelEnergy = aChannel->GetFragmentsEnergy(T);
00290
00291 return -MassExcess0 + Coulomb + ChannelEnergy;
00292
00293 }
00294
00295
00296