101 22.5, 22.0, 21.0, 21.0, 20.0,
104 20.6, 20.6, 18.6, 15.8, 13.5, 6.5,
107 6.65, 6.22, 6.27, 6.5, 6.7, 6.2, 6.25, 5.9, 6.1, 5.75,
110 6.46, 5.7, 6.28, 5.8, 6.15, 5.6, 5.8, 5.2, 5.8,
113 6.2 , 5.9 , 5.9, 6.0, 5.8, 5.7, 5.4, 5.4,
116 5.6, 6.1, 5.57, 6.3, 5.5, 5.8, 4.7, 6.2, 6.4, 6.2,
119 6.5, 6.2, 6.5, 5.3, 6.4, 5.7, 5.7, 6.2, 5.7,
122 6.3, 5.8, 6.7, 5.8, 6.6, 6.1, 4.3,
125 6.2, 3.8, 5.6, 4.0, 4.0, 4.2, 4.2, 3.5 };
129 0.6761, 0.677, 0.6788, 0.6803, 0.685,
131 0.6889, 0.6914, 0.6991, 0.7068, 0.725, 0.7391,
133 0.74, 0.741, 0.742, 0.743, 0.744, 0.7509, 0.752, 0.7531, 0.7543, 0.7548,
135 0.7557, 0.7566, 0.7576,
137 0.7587, 0.7597, 0.7608, 0.762, 0.7632, 0.7644, 0.7675, 0.7686, 0.7697, 0.7709,
139 0.7714, 0.7721, 0.7723, 0.7733, 0.7743, 0.7753, 0.7764,
141 0.7775, 0.7786, 0.7801, 0.781, 0.7821, 0.7831, 0.7842, 0.7852,
143 0.7864, 0.7875, 0.7880, 0.7887, 0.7889, 0.7899, 0.7909, 0.7919, 0.7930,
145 0.7941, 0.7953, 0.7965, 0.7977, 0.7987, 0.7989,
147 0.7997, 0.8075, 0.8097, 0.8119, 0.8143, 0.8164, 0.8174, 0.8274 };
155 theParaMaker(verboseLevel), QFinterp(XREP) {
156 parms.first.resize(6,0.);
157 parms.second.resize(6,0.);
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);
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
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];
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);
397 FMAX = EEXS*EEXS*EEXS*
EEXS;
401 G4cout <<
" T04 " << T04 <<
" FMAX (EEXS^4) " << FMAX <<
G4endl;
404 while (itry < itry_max) {
412 if (itry == itry_max) {
429 G4cout <<
" nucleus px " <<
PEX.px() <<
" py " <<
PEX.py()
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);
522 G4cout <<
" nucleus px " <<
PEX.px() <<
" py " <<
PEX.py()
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);
567 G4cout <<
" nucleus px " <<
PEX.px() <<
" py " <<
PEX.py()
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;
624 fission_open =
false;
632 if (itry_global == itry_global_max) {
634 G4cout <<
" ! itry_global " << itry_global_max <<
G4endl;
658 G4cout <<
" >>> G4EquilibriumEvaporator::explosion? ";
664 G4bool bigb = (!(a >= 12 && z >= 0 && z < 3*(a-z)) &&
676 G4cout <<
" >>> G4EquilibriumEvaporator::goodRemnant(" << a <<
"," << z
677 <<
")? " << (a>1 && z>0 && a>z) <<
G4endl;
680 return a > 1 && z > 0 && a > z;
689 G4cout <<
" >>> G4EquilibriumEvaporator::getQF ";
698 if (x < XMIN || x > XMAX) {
700 G4double FX = (0.73 + (3.33 * X1 - 0.66) * X1) * (X1*X1*X1);
702 QFF = G0 * FX * A13*
A13;
707 if (QFF < 0.0) QFF = 0.0;
720 G4cout <<
" >>> G4EquilibriumEvaporator::getAF" <<
G4endl;
724 G4double AF = 1.285 * (1.0 - e / 1100.0);
726 if(AF < 1.06) AF = 1.06;
736 G4cout <<
" >>> G4EquilibriumEvaporator::getPARLEVDEN" <<
G4endl;
747 G4cout <<
" >>> G4EquilibriumEvaporator::getE0" <<
G4endl;
virtual void setVerboseLevel(G4int verbose=0)
std::vector< ExP01TrackerHit * > a
G4double getQF(G4double x, G4double x2, G4int a, G4int z, G4double e) const
virtual ~G4EquilibriumEvaporator()
virtual G4bool explosion(G4int a, G4int z, G4double e) const
G4double getAF(G4double x, G4int a, G4int z, G4double e) const
G4CollisionOutput fission_output
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
static G4double getParticleMass(G4int type)
void setBullet(const G4InuclParticle *bullet)
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4CascadeInterpolator< 72 > QFinterp
G4bool goodRemnant(G4int a, G4int z) const
void getTargetData(const G4Fragment &target)
G4EquilibriumEvaporator()
G4GLOB_DLL std::ostream G4cout
static constexpr double m
G4double interpolate(const G4double x, const G4double(&yb)[nBins]) const
G4double getNucleiMass() const
std::pair< std::vector< G4double >, std::vector< G4double > > parms
void boostToLabFrame(const G4LorentzConvertor &convertor)
G4LorentzVector generateWithRandomAngles(G4double p, G4double mass=0.)
virtual void setVerboseLevel(G4int verbose)
G4double G4cbrt(G4double x)
G4double getE0(G4int A) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
const G4Fragment & makeFragment(G4LorentzVector mom, G4int A, G4int Z, G4double EX=0.)
G4double getPARLEVDEN(G4int A, G4int Z) const
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
static const G4double * SL[nLA]
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
static constexpr double GeV
G4InuclSpecialFunctions::paraMaker theParaMaker
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)
CLHEP::HepLorentzVector G4LorentzVector