174 G4cout <<
" >>> G4EquilibriumEvaporator::deExcite" <<
G4endl;
193 const G4double prob_cut_off = 1.0e-15;
194 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
195 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
196 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
197 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
199 const G4double fisssion_cut = 1000.0;
200 const G4double cut_off_energy = 0.1;
203 const G4int itry_max = 1000;
204 const G4int itry_global_max = 1000;
206 const G4int itry_gam_max = 100;
221 toTheNucleiSystemRestFrame.
setBullet(dummy);
226 if (explosion(
A,
Z,
EEXS)) {
228 theBigBanger.
deExcite(target, output);
235 if (
EEXS < cut_off_energy) {
244 G4double coul_coeff = (
A >= 100.0) ? 1.4 : 1.2;
249 G4bool fission_open =
true;
250 G4int itry_global = 0;
253 while (try_again && itry_global < itry_global_max) {
262 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
266 if (explosion(
A,
Z,
EEXS)) {
268 G4cout <<
" big bang in eql step " << itry_global <<
G4endl;
276 if (
EEXS < cut_off_energy) {
278 G4cout <<
" no energy for evaporation in eql step " << itry_global
291 const std::vector<G4double>& AK = parms.first;
292 const std::vector<G4double>& CPA = parms.second;
297 for (i = 0; i < 6; i++) {
300 u[i] = parlev * A1[i];
304 if (goodRemnant(A1[i], Z1[i])) {
306 V[i] = coul_coeff *
Z * Q[i] * AK[i] / (1.0 +
EEXS / E0) /
308 TM[i] =
EEXS - QB - V[i] * A / A1[i];
310 G4cout <<
" i=" << i <<
" : QB " << QB <<
" u " << u[i]
311 <<
" V " << V[i] <<
" TM " << TM[i] <<
G4endl;
320 if (TM[0] > cut_off_energy) {
323 W[0] = BE * A13*A13 * G[0] * AL;
324 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
326 if (TM1 > huge_num) TM1 = huge_num;
327 else if (TM1 < small) TM1 = small;
333 for (i = 1; i < 6; i++) {
335 if (TM[i] > cut_off_energy) {
337 W[i] = BE * A13*A13 * G[i] * (1.0 + CPA[i]);
338 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
340 if (TM1 > huge_num) TM1 = huge_num;
341 else if (TM1 < small) TM1 = small;
350 if (A >= 100.0 && fission_open) {
353 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
358 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
360 if (TM1 > huge_num) TM1 = huge_num;
361 else if (TM1 < small) TM1 = small;
363 W[6] = BF *
G4Exp(TM1);
364 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
372 G4cout <<
" Evaporation probabilities: sum = " << prob_sum
373 <<
"\n\t n " << W[0] <<
" p " << W[1] <<
" d " << W[2]
374 <<
" He3 " << W[3] <<
"\n\t t " << W[4] <<
" alpha " << W[5]
375 <<
" fission " << W[6] <<
G4endl;
380 if (prob_sum < prob_cut_off) {
382 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
388 while (
EEXS > cut_off_energy && try_again) {
395 FMAX = (T04*T04*T04*T04) *
G4Exp((
EEXS - T04) / T00);
401 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
404 while (itry < itry_max) {
412 if (itry == itry_max) {
430 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
431 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
432 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
447 if (itry_gam == itry_gam_max) try_again =
false;
456 for (i = 0; i < 7; i++) {
464 if (icase < 0)
continue;
468 G4cout <<
" particle/light-ion escape ("
469 << (icase==0 ?
"neutron" : icase==1 ?
"proton" :
470 icase==2 ?
"deuteron" : icase==3 ?
"He3" :
471 icase==4 ?
"triton" : icase==5 ?
"alpha" :
"ERROR!")
475 G4double Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
479 while (itry1 < itry_max && bad) {
486 while (itry < itry_max) {
491 Ptest = (X/Xmax)*
G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
498 if (S > V[icase] && S <
EEXS) {
501 G4cout <<
" escape itry1 " << itry1 <<
" icase "
502 << icase <<
" S (MeV) " << S <<
G4endl;
507 G4int ptype = 2 - icase;
515 G4double pmod = std::sqrt((2.0 * mass + S) * S);
523 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
524 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
525 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
533 if (EEXS_new < 0.0)
continue;
552 G4cout <<
" nucleus A " << AN[icase] <<
" Z " << Q[icase]
553 <<
" escape icase " << icase <<
G4endl;
560 G4double pmod = std::sqrt((2.0 * mass + S) * S);
568 <<
" pz " <<
PEX.
pz() <<
" E " <<
PEX.
e() << G4endl
569 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
570 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
578 if (EEXS_new < 0.0)
continue;
588 AN[icase], Q[icase], 0.*GeV,
600 if (itry1 == itry_max || bad) try_again =
false;
603 G4cout <<
" fission: A " << A <<
" Z " <<
Z <<
" eexs " <<
EEXS
604 <<
" Wn " << W[0] <<
" Wf " << W[6] <<
G4endl;
608 fission_output.
reset();
624 fission_open =
false;
632 if (itry_global == itry_global_max) {
634 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void getTargetData(const G4Fragment &target)
G4GLOB_DLL std::ostream G4cout
static constexpr double m
G4double getNucleiMass() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4cbrt(G4double x)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
static const G4double * SL[nLA]
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
static constexpr double GeV
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
G4int numberOfFragments() const
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
void addRecoilFragment(const G4Fragment *aFragment)
const G4Fragment & getRecoilFragment(G4int index=0) const
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void toTheTargetRestFrame()
void getParams(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4double bindingEnergy(G4int A, G4int Z)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void setTarget(const G4InuclParticle *target)