76 G4int IterationsLimit = 100000;
86 G4double theMeanMult = theMicrocanonicalEnsemble->GetMeanMultiplicity();
87 if (theMeanMult <= MaxAverageMultiplicity) {
90 theChannel = theMicrocanonicalEnsemble->ChooseAandZ(theFragment);
91 _theEnsemble = theMicrocanonicalEnsemble;
99 _theEnsemble = theMacrocanonicalEnsemble;
105 theChannel = theMacrocanonicalEnsemble->
ChooseAandZ(theFragment);
109 if (!ChannelOk)
delete theChannel;
111 }
while (!ChannelOk);
116 theResult->push_back(
new G4Fragment(theFragment));
117 delete theMicrocanonicalEnsemble;
118 if (theMacrocanonicalEnsemble != 0)
delete theMacrocanonicalEnsemble;
130 if (FindTemperatureOfBreakingChannel(theFragment,theChannel,Temperature))
break;
140 }
while (Iterations++ < IterationsLimit );
145 if (Iterations >= IterationsLimit)
146 throw G4HadronicException(__FILE__, __LINE__,
"G4StatMF::BreakItUp: Was not possible to solve for temperature of breaking channel");
157 InitialMomentum.
boost(-InitialMomentum.boostVector());
162 G4FragmentVector::iterator j;
163 for (j = theResult->begin(); j != theResult->end(); j++)
164 FragmentsEnergy += (*j)->GetMomentum().
e();
165 SavedScaleFactor = ScaleFactor;
166 ScaleFactor = InitialMomentum.e()/FragmentsEnergy;
168 for (j = theResult->begin(); j != theResult->end(); j++) {
169 ScaledMomentum = ScaleFactor * (*j)->GetMomentum().vect();
170 G4double Mass = (*j)->GetMomentum().m();
172 NewMomentum.
setVect(ScaledMomentum);
173 NewMomentum.
setE(std::sqrt(ScaledMomentum.
mag2()+Mass*Mass));
174 (*j)->SetMomentum(NewMomentum);
176 }
while (ScaleFactor > 1.0+1.e-5 && std::abs(ScaleFactor-SavedScaleFactor)/ScaleFactor > 1.e-10);
180 G4FragmentVector::iterator i;
181 for (i = theResult->begin(); i != theResult->end(); i++) {
184 (*i)->SetMomentum(FourMom);
185 #ifdef PRECOMPOUND_TEST
186 (*i)->SetCreatorModel(
G4String(
"G4StatMF"));
191 delete theMicrocanonicalEnsemble;
192 if (theMacrocanonicalEnsemble != 0)
delete theMacrocanonicalEnsemble;
199 G4bool G4StatMF::FindTemperatureOfBreakingChannel(
const G4Fragment & theFragment,
214 G4double TotalEnergy = CalcEnergy(A,Z,aChannel,T);
223 }
else if (Da < 0.0) {
225 Tb -= 0.5 * std::fabs(Tb);
227 if (Tb < 0.001*
MeV)
return false;
229 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
231 Db = (U - TotalEnergy)/U;
236 Tb += 0.5 * std::fabs(Tb);
239 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
241 Db = (U - TotalEnergy)/U;
245 G4double eps = 1.0e-14 * std::abs(Tb-Ta);
249 for (
G4int j = 0; j < 1000; j++) {
251 if (std::abs(Ta-Tc) <= eps) {
258 TotalEnergy = CalcEnergy(A,Z,aChannel,T);
276 Temperature = (Ta+Tb)/2.0;
291 return -MassExcess0 + Coulomb + ChannelEnergy;
static G4Pow * GetInstance()
Hep3Vector boostVector() const
static G4double GetKappaCoulomb()
G4double GetFragmentsEnergy(G4double T) const
G4double GetMeanTemperature(void) const
size_t GetMultiplicity(void)
G4StatMFChannel * ChooseAandZ(const G4Fragment &theFragment)
static G4double GetMaxAverageMultiplicity(G4int A)
static G4double GetMassExcess(const G4int A, const G4int Z)
G4double Z13(G4int Z) const
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
const G4LorentzVector & GetMomentum() const
HepLorentzVector & boost(double, double, double)
std::vector< G4Fragment * > G4FragmentVector
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void setVect(const Hep3Vector &)
G4bool CheckFragments(void)
G4double GetExcitationEnergy() const