100 22.5, 22.0, 21.0, 21.0, 20.0,
103 20.6, 20.6, 18.6, 15.8, 13.5, 6.5,
106 6.65, 6.22, 6.27, 6.5, 6.7, 6.2, 6.25, 5.9, 6.1, 5.75,
109 6.46, 5.7, 6.28, 5.8, 6.15, 5.6, 5.8, 5.2, 5.8,
112 6.2 , 5.9 , 5.9, 6.0, 5.8, 5.7, 5.4, 5.4,
115 5.6, 6.1, 5.57, 6.3, 5.5, 5.8, 4.7, 6.2, 6.4, 6.2,
118 6.5, 6.2, 6.5, 5.3, 6.4, 5.7, 5.7, 6.2, 5.7,
121 6.3, 5.8, 6.7, 5.8, 6.6, 6.1, 4.3,
124 6.2, 3.8, 5.6, 4.0, 4.0, 4.2, 4.2, 3.5 };
128 0.6761, 0.677, 0.6788, 0.6803, 0.685,
130 0.6889, 0.6914, 0.6991, 0.7068, 0.725, 0.7391,
132 0.74, 0.741, 0.742, 0.743, 0.744, 0.7509, 0.752, 0.7531, 0.7543, 0.7548,
134 0.7557, 0.7566, 0.7576,
136 0.7587, 0.7597, 0.7608, 0.762, 0.7632, 0.7644, 0.7675, 0.7686, 0.7697, 0.7709,
138 0.7714, 0.7721, 0.7723, 0.7733, 0.7743, 0.7753, 0.7764,
140 0.7775, 0.7786, 0.7801, 0.781, 0.7821, 0.7831, 0.7842, 0.7852,
142 0.7864, 0.7875, 0.7880, 0.7887, 0.7889, 0.7899, 0.7909, 0.7919, 0.7930,
144 0.7941, 0.7953, 0.7965, 0.7977, 0.7987, 0.7989,
146 0.7997, 0.8075, 0.8097, 0.8119, 0.8143, 0.8164, 0.8174, 0.8274 };
154 theParaMaker(verboseLevel), QFinterp(XREP) {
155 parms.first.resize(6,0.);
156 parms.second.resize(6,0.);
173 G4cout <<
" >>> G4EquilibriumEvaporator::deExcite" <<
G4endl;
192 const G4double prob_cut_off = 1.0e-15;
193 const G4double Q1[6] = { 0.0, 0.0, 2.23, 8.49, 7.72, 28.3 };
194 const G4int AN[6] = { 1, 1, 2, 3, 3, 4 };
195 const G4int Q[6] = { 0, 1, 1, 1, 2, 2 };
196 const G4double G[6] = { 2.0, 2.0, 6.0, 6.0, 6.0, 4.0 };
198 const G4double fisssion_cut = 1000.0;
199 const G4double cut_off_energy = 0.1;
202 const G4int itry_max = 1000;
203 const G4int itry_global_max = 1000;
205 const G4int itry_gam_max = 100;
220 toTheNucleiSystemRestFrame.
setBullet(dummy);
225 if (explosion(
A,
Z,
EEXS)) {
227 theBigBanger.
deExcite(target, output);
234 if (
EEXS < cut_off_energy) {
243 G4double coul_coeff = (
A >= 100.0) ? 1.4 : 1.2;
248 G4bool fission_open =
true;
249 G4int itry_global = 0;
252 while (try_again && itry_global < itry_global_max) {
261 G4cout <<
" A " <<
A <<
" Z " <<
Z <<
" mass " << nuc_mass
265 if (explosion(
A,
Z,
EEXS)) {
267 G4cout <<
" big bang in eql step " << itry_global <<
G4endl;
275 if (
EEXS < cut_off_energy) {
277 G4cout <<
" no energy for evaporation in eql step " << itry_global
290 const std::vector<G4double>& AK = parms.first;
291 const std::vector<G4double>& CPA = parms.second;
296 for (i = 0; i < 6; i++) {
299 u[i] = parlev * A1[i];
303 if (goodRemnant(A1[i], Z1[i])) {
305 V[i] = coul_coeff *
Z * Q[i] * AK[i] / (1.0 +
EEXS / E0) /
307 TM[i] =
EEXS - QB - V[i] *
A / A1[i];
309 G4cout <<
" i=" << i <<
" : QB " << QB <<
" u " << u[i]
310 <<
" V " << V[i] <<
" TM " << TM[i] <<
G4endl;
319 if (TM[0] > cut_off_energy) {
322 W[0] = BE *
A13*
A13 * G[0] * AL;
323 G4double TM1 = 2.0 * std::sqrt(u[0] * TM[0]) - ue;
325 if (TM1 > huge_num) TM1 = huge_num;
326 else if (TM1 < small) TM1 = small;
332 for (i = 1; i < 6; i++) {
334 if (TM[i] > cut_off_energy) {
336 W[i] = BE *
A13*
A13 * G[i] * (1.0 + CPA[i]);
337 G4double TM1 = 2.0 * std::sqrt(u[i] * TM[i]) - ue;
339 if (TM1 > huge_num) TM1 = huge_num;
340 else if (TM1 < small) TM1 = small;
349 if (
A >= 100.0 && fission_open) {
352 G4double X = 0.019316 * X2 / (1.0 - 1.79 * X1 * X1);
357 G4double TM1 = 2.0 * std::sqrt(AF * EF) - ue;
359 if (TM1 > huge_num) TM1 = huge_num;
360 else if (TM1 < small) TM1 = small;
362 W[6] = BF *
G4Exp(TM1);
363 if (W[6] > fisssion_cut*W[0]) W[6] = fisssion_cut*W[0];
371 G4cout <<
" Evaporation probabilities: sum = " << prob_sum
372 <<
"\n\t n " << W[0] <<
" p " << W[1] <<
" d " << W[2]
373 <<
" He3 " << W[3] <<
"\n\t t " << W[4] <<
" alpha " << W[5]
374 <<
" fission " << W[6] <<
G4endl;
379 if (prob_sum < prob_cut_off) {
381 G4double T00 = 1.0 / (std::sqrt(u1 / UCR0) - 1.25 / UCR0);
387 while (
EEXS > cut_off_energy && try_again) {
394 FMAX = (T04*T04*T04*T04) *
G4Exp((
EEXS - T04) / T00);
400 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
403 while (itry < itry_max) {
411 if (itry == itry_max) {
430 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
431 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
446 if (itry_gam == itry_gam_max) try_again =
false;
455 for (i = 0; i < 7; i++) {
463 if (icase < 0)
continue;
467 G4cout <<
" particle/light-ion escape ("
468 << (icase==0 ?
"neutron" : icase==1 ?
"proton" :
469 icase==2 ?
"deuteron" : icase==3 ?
"He3" :
470 icase==4 ?
"triton" : icase==5 ?
"alpha" :
"ERROR!")
474 G4double Xmax = (std::sqrt(u[icase]*TM[icase] + 0.25) - 0.5)/u[icase];
478 while (itry1 < itry_max && bad) {
485 while (itry < itry_max) {
490 Ptest = (X/Xmax)*
G4Exp(-2.*u[icase]*Xmax + 2.*std::sqrt(u[icase]*(TM[icase] - X)));
497 if (
S > V[icase] &&
S <
EEXS) {
500 G4cout <<
" escape itry1 " << itry1 <<
" icase "
501 << icase <<
" S (MeV) " <<
S <<
G4endl;
506 G4int ptype = 2 - icase;
514 G4double pmod = std::sqrt((2.0 * mass +
S) *
S);
523 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
524 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
531 G4double EEXS_new = ((
PEX-mom).m() - mass_new)*GeV;
532 if (EEXS_new < 0.0)
continue;
551 G4cout <<
" nucleus A " << AN[icase] <<
" Z " << Q[icase]
552 <<
" escape icase " << icase <<
G4endl;
559 G4double pmod = std::sqrt((2.0 * mass +
S) *
S);
568 <<
" evaporate px " << mom.
px() <<
" py " << mom.
py()
569 <<
" pz " << mom.
pz() <<
" E " << mom.
e() <<
G4endl;
576 G4double EEXS_new = ((
PEX-mom).m() - mass_new)*GeV;
577 if (EEXS_new < 0.0)
continue;
587 AN[icase], Q[icase], 0.*GeV,
599 if (itry1 == itry_max || bad) try_again =
false;
602 G4cout <<
" fission: A " <<
A <<
" Z " <<
Z <<
" eexs " <<
EEXS
603 <<
" Wn " << W[0] <<
" Wf " << W[6] <<
G4endl;
607 fission_output.
reset();
623 fission_open =
false;
631 if (itry_global == itry_global_max) {
633 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
657 G4cout <<
" >>> G4EquilibriumEvaporator::explosion? ";
663 G4bool bigb = (!(a >= 12 && z >= 0 && z < 3*(a-z)) &&
672G4bool G4EquilibriumEvaporator::goodRemnant(
G4int a,
675 G4cout <<
" >>> G4EquilibriumEvaporator::goodRemnant(" << a <<
"," << z
676 <<
")? " << (a>1 && z>0 && a>z) <<
G4endl;
679 return a > 1 && z > 0 && a > z;
688 G4cout <<
" >>> G4EquilibriumEvaporator::getQF ";
697 if (x < XMIN || x > XMAX) {
699 G4double FX = (0.73 + (3.33 * X1 - 0.66) * X1) * (X1*X1*X1);
706 if (QFF < 0.0) QFF = 0.0;
719 G4cout <<
" >>> G4EquilibriumEvaporator::getAF" <<
G4endl;
723 G4double AF = 1.285 * (1.0 - e / 1100.0);
725 if(AF < 1.06) AF = 1.06;
735 G4cout <<
" >>> G4EquilibriumEvaporator::getPARLEVDEN" <<
G4endl;
746 G4cout <<
" >>> G4EquilibriumEvaporator::getE0" <<
G4endl;
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4GLOB_DLL std::ostream G4cout
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
void getTargetData(const G4Fragment &target)
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
virtual void setVerboseLevel(G4int verbose=0)
virtual G4bool validateOutput(const G4Fragment &target, G4CollisionOutput &output)
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
void addRecoilFragment(const G4Fragment *aFragment)
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
void boostToLabFrame(const G4LorentzConvertor &convertor)
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
const G4Fragment & getRecoilFragment(G4int index=0) const
void addOutgoingNucleus(const G4InuclNuclei &nuclei)
G4int numberOfFragments() const
virtual ~G4EquilibriumEvaporator()
virtual void setVerboseLevel(G4int verbose)
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
G4EquilibriumEvaporator()
virtual void deExcite(const G4Fragment &target, G4CollisionOutput &output)
static G4double getParticleMass(G4int type)
G4double getNucleiMass() const
void getParams(G4double Z, std::pair< std::vector< G4double >, std::vector< G4double > > &parms)
void toTheTargetRestFrame()
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
void setTarget(const G4InuclParticle *target)
G4double bindingEnergy(G4int A, G4int Z)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
G4double G4cbrt(G4double x)